How to compare two double precision real numbers in Fortran?

I wrote a Fortran subroutine to compute the time of flight(TOF) between two points on elliptical orbit. Obviously this TOF has to be a positive number. I tested my subroutine with different data but in some cases I get negative result despite of I coded a possible way to solve this problem.

Here is my subroutine:

!*************************************************************
    subroutine TOF_between_2TA(e,xn,theta1,theta2,delta_t)
!*************************************************************
!  Compute time of flight bewteen two points (point 1 and point 2) on elliptical orbits.
!*************************************************************
!  Input:
!  xn     = mean motion, any units
!  e      = eccentricity
!  theta1 = true anomaly of first point, in interval [0, twopi]
!  theta2 = true anomaly of second point, in interval [0, twopi]
!
!  Output:
!  delta_t = time of flight between two points (same units as given by xn)
!*************************************************************
 
       
        implicit none

        ! Arguments
        double precision, intent(in)  :: e,xn,theta1,theta2
        double precision, intent(out) :: delta_t

        ! Locals
        integer :: i
        double precision               :: xe,sxe,cxe,cth,sth,den,theta
        double precision, dimension(2) :: theta_vec,xm_vec
        
        !! To get positive time intervals, theta_vec must be sorted in ascending order  
        if (theta1 < theta2) then     
           theta_vec = [theta1, theta2]   ! Case theta1 < theta2
        elseif (theta2 < theta1) then 
           theta_vec = [theta2, theta1]  ! Case theta1 > theta2
        endif
        
        do i=1,2          
            theta = theta_vec(i)
            
            cth = cos(theta)
            sth = sin(theta)
            den = 1.0 + e*cth
            cxe = (e + cth)/den
            sxe = sqrt(1.0-e*e)*sth/den 
            xe = atan2(sxe,cxe)
        
            ! atan2 returns angles in interval -pi, +pi, we need angles in interval [0,2pi]
            xe = mod(xe,twopi)
            if (xe .lt. 0.0_pr) xe = xe + twopi 
        
            xm_vec(i) = xe - e*sxe         
            
        enddo
        
        delta_t= (xm_vec(2)-xm_vec(1)) / xn
        
        return
        
    end subroutine TOF_between_2TA

Can you suggest me a way to increase the robustness of my soubrutine and protect me against undesiderable result, i.e. negative numbers? I'm almost pretty sure that the problem is when I try to compare the vatiables theta1 and theta2, what is the smartest way to compare two real numbers?

Edit - After many tests, I report here the working routines:

!*************************************************************
    subroutine TOF_since_pericenter(e,xn,theta,delta_t_peri)
!*************************************************************
!  Time of Flight since pericenter: From true anomaly to time
!*************************************************************
!  Input:  
!  xn    = mean motion, any units
!  e     = eccentricity
!  theta = true anomaly, in interval [0, twopi]
!
!  Output:
!  delta_t_peri = time since pericenter (same units as given by xn)
!*************************************************************
 
        implicit none

        ! Arguments
        double precision, intent(in)  :: e,xn,theta
        double precision, intent(out) :: delta_t_peri

        ! Locals
        double precision :: xe,sxe,cxe,cth,sth,xm,den
        
        cth = cos(theta)
        sth = sin(theta)
        den = 1.d0 + e*cth
        cxe = (e + cth)/den
        sxe = sqrt(1.d0 - e*e)*sth/den 
        xe  = atan2(sxe,cxe)
        ! atan2 returns angles in interval [-pi, +pi], we need angles in interval [0,2pi]
        xe = mod(xe,twopi)
        if (xe .lt. 0.d0) xe = xe + twopi 
        
        xm = xe - e*sxe
        delta_t_peri = xm/xn
        
        return

    end subroutine ToF_since_pericenter
    
    
    
!*************************************************************
    subroutine TOF_between_2TA(e,xn,theta_1,theta_2,delta_t)
!*************************************************************
!  Compute time of flight required by point 1 to reach the point 2 on elliptical orbits.
!  The point 1 is the moving object and it rotates counterclockwise on its orbit 
!*************************************************************
!  Input:
!  xn        = mean motion, any units
!  e         = eccentricity
!  theta_1   = true anomaly of asteroid, in interval [0, twopi]
!  theta_2   = true anomaly of intersection point, in interval [0, twopi]
!
!  Output:
!  delta_t = time of flight between two points (same units as given by xn)
!*************************************************************
 
       
        implicit none

        ! Arguments
        double precision, intent(in)  :: e,xn,theta_1,theta_2
        double precision, intent(out) :: delta_t

        ! Locals
        integer :: i
        double precision               :: T,delta_t1_peri,delta_t2_peri
        double precision, parameter    :: small = 1.d-30
        
        ! Compute orbital period
        T = twopi / xn
        
        call TOF_since_pericenter(e,xn,theta_1,delta_t1_peri)
        call TOF_since_pericenter(e,xn,theta_2,delta_t2_peri)
        
        delta_t = delta_t2_peri - delta_t1_peri
        
        if (delta_t .gt. small) then     
           delta_t = delta_t
        elseif (delta_t .lt. -small) then 
           delta_t = delta_t + T
        endif
           
        return
        
    end subroutine TOF_between_2TA
 
!*************************************************************

Further comments and improvment suggestions are well accepted.