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