Attempt to add length relations based on fits to Jake's data; current…

…ly bugged
MintonGroup · Nov 10, 2023 · 145881a · 145881a
1 parent 0d45f8c
commit 145881a
Show file tree

Hide file tree

Showing 6 changed files with 34 additions and 14 deletions.
diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90
@@ -51,13 +51,16 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
    integer(I8B),dimension(user%gridsize,user%gridsize) :: Ngrid
 
    ! Crater ray parameters
-   real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019)
-   integer(I4B) :: Nraymax = 14
+   real(DP) :: rray != 48_DP ! "L16" in Minton et al. (2019)
+   integer(I4B) :: Nraymax = 12
    real(DP) :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
    real(DP) :: rayp = 2.0_DP 
    integer(I4B) :: rayq = 4
    real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP))
+   real(DP) :: l1
 
+   rray = 11.95*crater%frad**1.32
+   l1 = 5.32*crater%frad**1.27
 
    ! Create box for soften calculation (will be no bigger than the grid itself)
    do j = 0,user%gridsize + 1
@@ -184,7 +187,7 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
 
             allocate(diffdistribution(imin:imax,jmin:jmax))
             allocate(ejdistribution(imin:imax,jmin:jmax))
-            call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
+            call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1)
             ! Loop over affected matrix area
             !!$OMP PARALLEL DO DEFAULT(SHARED) IF(inc > INCPAR) &
             !!$OMP FIRSTPRIVATE(jmin,jmax,imin,imax) &

diff --git a/src/crater/crater_superdomain.f90 b/src/crater/crater_superdomain.f90
@@ -52,12 +52,16 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval)
    real(SP) :: gglass, glrad
 
       ! Crater ray parameters
-   real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019)
+   real(DP) :: rray != 48_DP ! "L16" in Minton et al. (2019)
    integer(I4B) :: Nraymax = 14
    real(DP) :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
    real(DP) :: rayp = 2.0_DP 
    integer(I4B) :: rayq = 4
    real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP))
+   real(DP) :: l1
+
+   rray = 11.95*crater%frad**1.32
+   l1 = 5.32*crater%frad**1.27
 
    ! Create box for soften calculation (will be no bigger than the grid itself)
    do j = 0,user%gridsize + 1
@@ -156,7 +160,7 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval)
          allocate(diffdistribution(xi:xf,yi:yf))
          allocate(ejisray(xi:xf,yi:yf))
          ! Now generate ray pattern
-         call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
+         call ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1)
          ! Now if doregotrack is on, do melt zone calculation
          if (user%doregotrack) call regolith_melt_zone_superdomain(user,crater,domain,rm,depthb)
 

diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90
@@ -129,12 +129,18 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
    real(SP) :: age_mean
 
    ! Crater ray parameters
-   real(DP) :: rray = 48_DP ! "L16" in Minton et al. (2019)
+   real(DP) :: rray != 48_DP ! "L16" in Minton et al. (2019)
    integer(I4B) :: Nraymax = 12
    real(DP) :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4
    real(DP) :: rayp = 2.0_DP 
    integer(I4B) :: rayq = 4
    real(DP) :: rayfmult = (5)**(-4.0_DP / (1.2_DP))
+   real(DP) :: l1
+
+
+   rray = 11.95*crater%frad**1.32
+   l1 = 5.32*crater%frad**1.27
+   write(*,*) "L16 = ", rray, "; L1 = ", l1
 
 
    ! Executable code
@@ -190,7 +196,7 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ
 
    ! *************************** Layered Ejecta Rays *****************************!
    do i=1,Npatt
-      call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,rray,Nraymax+i,fpeak,rayp,rayq,rayfmult,tempdiff,tempej)
+      call ejecta_ray_pattern(user,surf,crater,inc,-inc,inc,-inc,inc,rray,Nraymax+i,fpeak,rayp,rayq,rayfmult,tempdiff,tempej,l1)
       diffdistribution(:,:) = diffdistribution(:,:) + frayreduction**(i-1) * tempdiff(:,:)
       ejdistribution(:,:) = ejdistribution(:,:) + frayreduction**(i-1) * tempej(:,:)
    end do

diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90
@@ -31,7 +31,7 @@
 !***
 
 !**********************************************************************************************************************************
-subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
+subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1)
    use module_globals
    use module_util
    use module_io
@@ -49,6 +49,7 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea
    integer(I4B),intent(in) :: Nraymax, rayq
    real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution
    real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution
+   real(DP),intent(in) :: l1
 
    ! Internal variables
    integer(I4B) :: nrays,i,j,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq
@@ -87,6 +88,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea
       call random_number(rn) ! randomize the ray orientation
       rmax = user%ejecta_truncation 
       rmin = crater%continuous / crater%frad
+      !rmax = rray / crater%frad
+      !rmin = l1 / crater%frad
       crater%fe = 10.0_DP ! Estimate the equivalent degradation radius
       !ejdistribution = 0.0_DP
       !diffdistribution = 0.0_DP
@@ -108,8 +111,8 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea
             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_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,.false.) 
-            ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,.true.) 
+            diffdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,.false.) 
+            ejdistribution(i,j) = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,.true.) 
          end do
       end do
       !!$OMP END PARALLEL DO

diff --git a/src/ejecta/ejecta_ray_pattern_func.f90 b/src/ejecta/ejecta_ray_pattern_func.f90
@@ -27,7 +27,7 @@
 !***
 
 !**********************************************************************************************************
-function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,ej) result(ans)
+function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,ej) result(ans)
     use module_globals
     use module_ejecta, EXCEPT_THIS_ONE => ejecta_ray_pattern_func
     implicit none
@@ -36,6 +36,7 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,ra
     real(DP),dimension(:),intent(in) :: thetari
     real(DP),intent(in) :: rray, fpeak, rayp, rayfmult
     integer(I4B),intent(in) :: Nraymax, rayq
+    real(DP),intent(in) :: l1
     logical,intent(in) :: ej
     real(DP) :: a,c
     real(DP) :: thetar,rw,rw0,rw1
@@ -45,7 +46,8 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,ra
 
 
     !minray = rmin * 3 !"L1" in Minton et al. (2019)
-    minray = rray / 2
+    !minray = rray / 2
+    minray = l1
 
     if (r > rmax) then
        ans = 0._DP

diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90
@@ -43,7 +43,7 @@ end subroutine ejecta_emplace
    end interface
 
    interface
-      subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution)
+      subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpeak,rayp,rayq,rayfmult,diffdistribution,ejdistribution,l1)
       use module_globals
       implicit none
       type(usertype),intent(in) :: user
@@ -54,19 +54,21 @@ subroutine ejecta_ray_pattern(user,surf,crater,inc,xi,xf,yi,yf,rray,Nraymax,fpea
       integer(I4B),intent(in) :: Nraymax, rayq
       real(DP),dimension(xi:xf,yi:yf),intent(out) :: diffdistribution
       real(DP),dimension(xi:xf,yi:yf),intent(out) :: ejdistribution
+      real(DP),intent(in) :: l1
 
       end subroutine ejecta_ray_pattern
    end interface
 
    interface
-      function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,ej) result(ans)
+      function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,rray,Nraymax,fpeak,rayp,rayq,rayfmult,l1,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
       real(DP),intent(in) :: rray, fpeak, rayp, rayfmult
       integer(I4B),intent(in) :: Nraymax, rayq
+      real(DP),intent(in) :: l1
       logical,intent(in) :: ej
       end function ejecta_ray_pattern_func
    end interface