如何比较 Fortran 中的两个双精度实数？

Question

我编写了一个 Fortran 子程序来计算椭圆轨道上两点之间的飞行时间 (TOF)。显然这个 TOF 必须是正数。我用不同的数据测试了我的子例程，但在某些情况下，尽管我编写了一种可能的方法来解决这个问题，但我得到了负面结果。

这是我的子程序：

!*************************************************************
    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

您能否建议我一种方法来增加我的 soubrutine 的稳健性并保护我免受不良结果（即负数）的影响？我几乎可以肯定问题是当我尝试比较变量 theta1 和 theta2 时，比较两个实数的最聪明的方法是什么？

编辑 - 经过多次测试，我在这里报告工作例程：

!*************************************************************
    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
 
!*************************************************************

接受进一步的评论和改进建议。

Answer 1

您不需要进行明确的比较。您可以将 theta_vec 构造为

theta_vec = [min(theta1, theta2), max(theta1, theta2)]

这具有无分支的优势，因此运行比比较快得多。（并不是说这可能与代码的整体性能相关）。

我还应该指出行

xe = mod(xe,twopi)
if (xe .lt. 0.d0) xe = xe + twopi

可以换成

xe = modulo(xe, twopi)

利用 modulo(a, p) 的输出总是 non-negative 如果 p 为正，这与 mod 的输出不同。

如何比较 Fortran 中的两个双精度实数？

How to compare two double precision real numbers in Fortran?

floating-point

fortran

compare