From 058d20de58e691f1708972f2b74b99a7f5029ea3 Mon Sep 17 00:00:00 2001 From: Austin Michael Blevins Date: Tue, 15 Nov 2022 14:46:03 -0500 Subject: [PATCH] moved function ray() to its own file for debugging --- src/Makefile.am | 3 +- src/ejecta/ejecta_ray_pattern.f90 | 4 +- src/ejecta/ejecta_ray_pattern_function.f90 | 99 ---------------------- src/ejecta/module_ejecta.f90 | 14 ++- 4 files changed, 16 insertions(+), 104 deletions(-) delete mode 100644 src/ejecta/ejecta_ray_pattern_function.f90 diff --git a/src/Makefile.am b/src/Makefile.am index 8d53b108..56df03dd 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -65,7 +65,8 @@ io/io_write_pindex_map.f90\ io/io_write_age_depth.f90\ ejecta/ejecta_emplace.f90\ ejecta/ejecta_ray_pattern.f90\ -ejecta/ejecta_ray_pattern_function.f90\ +ejecta/ejecta_ray_pattern_func.f90\ +ejecta/ejecta_ray_ray_func.f90\ ejecta/ejecta_blanket.f90\ ejecta/ejecta_blanket_func.f90\ ejecta/ejecta_thickness.f90\ diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index bc280982..352cea4a 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -106,8 +106,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,diffdistribution, areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) r = sqrt(xbar**2 + ybar**2) / crater%frad theta = mod(atan2(ybar,xbar) + pi + rn * 2 * pi,2 * pi) - diffdistribution(i,j) = areafrac * ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,.false.) - ejdistribution(i,j) = areafrac * ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,.true.) + diffdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.false.) + ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.true.) end do end do !!$OMP END PARALLEL DO diff --git a/src/ejecta/ejecta_ray_pattern_function.f90 b/src/ejecta/ejecta_ray_pattern_function.f90 deleted file mode 100644 index 26d00f81..00000000 --- a/src/ejecta/ejecta_ray_pattern_function.f90 +++ /dev/null @@ -1,99 +0,0 @@ -!****f* regolith/regolith_streamtube_volume_func -! Name -! ejecta_ray_pattern_func -- Calculate ejecta ray pattern -! SYNOPSIS -! This uses -! * module_globals -! * module_ejecta -! -! ans = regolith_streamtube_volume_func() -! -! DESCRIPTION -! -! This function was separated into an independent function for debugging purposes. -! -! ARGUMENTS -! Input -! * theta -- -! * r -- -! * rmin -- -! * rmax -- -! * thetari -- -! * ej -- -! -! Output -! * ans -- THe result of this function, used in ejecta_ray_pattern.f90 -! -!*** - -!********************************************************************************************************** -function ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,ej) result(ans) - use module_globals - use module_ejecta, EXCEPT_THIS_ONE => ejecta_ray_pattern_function - implicit none - real(DP) :: ans - real(DP),intent(in) :: r,rmin,rmax,theta - real(DP),dimension(:),intent(in) :: thetari - logical,intent(in) :: ej - real(DP) :: a,c - real(DP) :: thetar,rw,rw0,rw1 - real(DP) :: f,rtrans,length,rpeak,minray,FF - integer(I4B) :: n,i - - - minray = rmin * 3 - - if (r > rmax) then - ans = 0._DP - else if (r < 1.0_DP) then - if (ej) then - ans = 1.0_DP - else - ans = 0.0_DP - end if - else - rw0 = rmin * pi / Nraymax / 2 - rw1 = 2 * pi / Nraymax - rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2)) - n = max(min(floor((Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray))**(1._DP/rayp)),Nraymax),1) ! Exponential decay of ray number with distance - ans = 0._DP - rtrans = r - 1.0_DP - c = rw / r - a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c))) - do i = 1,Nraymax - length = minray * exp(log(rray/minray) * ((Nraymax - i + 1)**rayp - 1_DP) / ((Nraymax**rayp - 1))) - rpeak = (length - 1_DP) * 0.5_DP - if (ej) then - FF = 1.0_DP - if (r > length) then - f = 0.0_DP - else - f = a - end if - else - FF = rayfmult * (20 / rmax)**(0.5_DP) * 0.25_DP - f = FF * fpeak * (rtrans / rpeak)**rayq * exp(1._DP / rayq * (1.0_DP - (rtrans / rpeak)**rayq)) - end if - ans = ans + ray(theta,thetari(i),r,n,rw) * f / a - end do - end if - !return - contains - - pure function ray(theta,thetar,r,n,w) result(ans) - implicit none - real(DP) :: ans - real(DP),intent(in) :: theta,thetar,r,w - integer(I4B),intent(in) :: n - real(DP) :: thetap,thetapp,a,b,c,dtheta - - 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))) - ans = a * exp(-dtheta**2 / (2 * c**2)) - - return - end function ray - -end function ejecta_ray_pattern_function \ No newline at end of file diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 578af4fb..0a56819f 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -55,14 +55,24 @@ end subroutine ejecta_ray_pattern end interface interface - function ejecta_ray_pattern_function(theta,r,rmin,rmax,thetari,ej) result(ans) + function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans) use module_globals implicit none real(DP) :: ans real(DP),intent(in) :: r,rmin,rmax,theta real(DP),dimension(:),intent(in) :: thetari logical,intent(in) :: ej - end function ejecta_ray_pattern_function + end function ejecta_ray_pattern_func + end interface + + interface + function ejecta_ray_ray_func(theta,thetar,r,n,w) result(ans) + use module_globals + implicit none + real(DP) :: ans + real(DP),intent(in) :: theta,thetar,r,w + integer(I4B),intent(in) :: n + end function ejecta_ray_ray_func end interface interface