Skip to content

Commit

Permalink
Check to ensure that the value passed to exp doesn't cause an underfl…
Browse files Browse the repository at this point in the history
…ow exception
  • Loading branch information
daminton committed Nov 13, 2023
1 parent e6ec31e commit b199af2
Showing 1 changed file with 9 additions and 4 deletions.
13 changes: 9 additions & 4 deletions src/ejecta/ejecta_ray_func.f90
Original file line number Diff line number Diff line change
Expand Up @@ -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

0 comments on commit b199af2

Please sign in to comment.