Commit d310e9ea authored by William D. Fullmer's avatar William D. Fullmer
Browse files

simplify the R_d expression from Koch-Sangani to look more like F0 from BVK

parent ec4d4bb3
Loading
Loading
Loading
Loading
+18 −11
Original line number Diff line number Diff line
@@ -1215,19 +1215,26 @@ CONTAINS

! calculating the term phi*dRd/dphi
             dRdphi = zero
             IF((EP_SM > SMALL_NUMBER) .AND. (EP_SM <= 0.4d0)) THEN
                denom = one+0.681d0*EP_SM-8.48d0*EP_SM**2+8.16d0*EP_SM**3
                dRdphi = (1.5d0*dsqrt(EP_SM/2d0)+135d0/64d0*EP_SM*&
                   (dlog(EP_SM)+one)+ 17.14d0*EP_SM)/denom - EP_SM*&
                   (one+3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*&
                   dlog(EP_SM)+17.14*EP_SM)/denom**2 * &
                   (0.681d0-16.96d0*EP_SM+24.48d0*EP_SM**2)
             ELSEIF(EP_SM > 0.4d0) THEN
!wdf 
! missing leading phi from taking EP_SM inside derivative... 
                dRdphi = 10d0*EP_SM*(one+2d0*EP_SM)/(one-EP_SM)**4
!             IF((EP_SM > SMALL_NUMBER) .AND. (EP_SM <= 0.4d0)) THEN
!                denom = one+0.681d0*EP_SM-8.48d0*EP_SM**2+8.16d0*EP_SM**3
!                dRdphi = (1.5d0*dsqrt(EP_SM/2d0)+135d0/64d0*EP_SM*&
!                   (dlog(EP_SM)+one)+ 17.14d0*EP_SM)/denom - EP_SM*&
!                   (one+3d0*dsqrt(EP_SM/2d0) + 135d0/64d0*EP_SM*&
!                   dlog(EP_SM)+17.14*EP_SM)/denom**2 * &
!                   (0.681d0-16.96d0*EP_SM+24.48d0*EP_SM**2)
!             ELSEIF(EP_SM > 0.4d0) THEN
!!wdf 
!! missing leading phi from taking EP_SM inside derivative... 
!                dRdphi = 10d0*EP_SM*(one+2d0*EP_SM)/(one-EP_SM)**4
!!wdf
!             ENDIF
! use a simpler fit: 
              IF (EP_SM .GT. SMALL_NUMBER) dRdphi & !really phi*D(R_d)/D(phi) 
                = 10d0*EP_SM*(one+2d0*EP_SM)/(one-EP_SM)**4 &
                + 1.5d0*(1.0d0-EP_SM)*DSQRT(EP_SM) & 
                - EP_SM*(1.0d0 + 3.0d0*DSQRT(EP_SM))
!wdf
             ENDIF

! calculating the term phi*dS_star/dphi
             dSdphi = zero
+14 −7
Original line number Diff line number Diff line
@@ -807,13 +807,20 @@
      DOUBLE PRECISION, INTENT(IN) :: phi
!---------------------------------------------------------------------//

      R_d = 1.0d0  ! this avoids singularity at phi = 0.0
      if((phi > 1d-15) .and. (phi <= 0.4d0)) then
        R_d = (1d0+3d0*dsqrt(phi/2d0)+135d0/64d0*phi*dlog(phi)+17.14d0*phi) / &
              (1d0+0.681d0*phi-8.48*phi**2+8.16d0*phi**3)
      elseif(phi > 0.4d0) then
        R_d = 10d0*phi/(1d0-phi)**3 + 0.7d0
      endif
!wdf
!      R_d = 1.0d0  ! this avoids singularity at phi = 0.0
!      if((phi > 1d-15) .and. (phi <= 0.4d0)) then
!        R_d = (1d0+3d0*dsqrt(phi/2d0)+135d0/64d0*phi*dlog(phi)+17.14d0*phi) / &
!              (1d0+0.681d0*phi-8.48*phi**2+8.16d0*phi**3)
!      elseif(phi > 0.4d0) then
!        R_d = 10d0*phi/(1d0-phi)**3 + 0.7d0
!      endif
! use a simpler fit:
         R_d = 1.0d0
         IF (phi > 1.0d-15) &  
         R_d = 10d0*phi/(1d0-phi)**3 & 
             + (1.0d0 - phi)*(1.0d0 + 3.0d0*DSQRT(phi)) 
!wdf
      RETURN
      END FUNCTION R_d