From b199af274797eb82f6b1d3bfadf129f0528470d5 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 13 Nov 2023 13:23:53 -0500 Subject: [PATCH] Check to ensure that the value passed to exp doesn't cause an underflow exception --- src/ejecta/ejecta_ray_func.f90 | 13 +++++++++---- 1 file changed, 9 insertions(+), 4 deletions(-) diff --git a/src/ejecta/ejecta_ray_func.f90 b/src/ejecta/ejecta_ray_func.f90 index 459fefde..25320bec 100644 --- a/src/ejecta/ejecta_ray_func.f90 +++ b/src/ejecta/ejecta_ray_func.f90 @@ -33,13 +33,18 @@ function ejecta_ray_func(theta,thetar,r,n,w) result(ans) real(DP) :: ans real(DP),intent(in) :: theta,thetar,r,w integer(I4B),intent(in) :: n - real(DP) :: thetap,thetapp,a,b,c,dtheta + real(DP) :: thetap,thetapp,a,b,c,dtheta,logval c = w / r b = thetar dtheta = min(2*pi - abs(theta - b),abs(theta - b)) - a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c))) !this is the intensity function - ans = a * exp(-dtheta**2 / (2 * c**2)) + logval = -dtheta**2 / (2 * c**2) + if (logval < LOGVSMALL) then + ans = 0.0_DP + else + a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c))) !this is the intensity function + ans = a * exp(logval) + end if -!return + return end function ejecta_ray_func \ No newline at end of file