Merge branch 'realistic' into visualize3d

python/ctem/ctem/NPF.py

@@ -0,0 +1,104 @@
+import numpy as np
+# The Neukum production function
+# Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86
+def N1lunar(T):
+	return 5.44e-14 * (np.exp(6.93*T) - 1) + 8.38e-4 * T
+
+
+def Nlunar(D):
+    # Lunar crater SFD
+    aL00 = -3.0876
+    aL01 = -3.557528
+    aL02 =  0.781027
+    aL03 =  1.021521
+    aL04 = -0.156012
+    aL05 = -0.444058
+    aL06 =  0.019977
+    aL07 =  0.086850
+    aL08 = -0.005874
+    aL09 = -0.006809
+    aL10 =  8.25e-4
+    aL11 =  5.54e-5
+
+    return np.where((D > 0.01) & (D < 1000.0),
+        10**(aL00
+        + aL01 * np.log10(D) ** 1
+        + aL02 * np.log10(D) ** 2
+        + aL03 * np.log10(D) ** 3
+        + aL04 * np.log10(D) ** 4
+        + aL05 * np.log10(D) ** 5
+        + aL06 * np.log10(D) ** 6
+        + aL07 * np.log10(D) ** 7
+        + aL08 * np.log10(D) ** 8
+        + aL09 * np.log10(D) ** 9
+        + aL10 * np.log10(D) ** 10
+        + aL11 * np.log10(D) ** 11),float('nan'))
+
+def Rproj(D):
+    #Projectile SFD
+    aP00 =  0.0
+    aP01 = -1.375458
+    aP02 =  1.272521e-1
+    aP03 = -1.282166
+    aP04 = -3.074558e-1
+    aP05 =  4.149280e-1
+    aP06 =  1.910668e-1
+    aP07 = -4.260980e-2
+    aP08 = -3.976305e-2
+    aP09 = -3.180179e-3
+    aP10 =  2.799369e-3
+    aP11 =  6.892223e-4
+    aP12 =  2.614385e-6
+    aP13 = -1.416178e-5
+    aP14 = -1.191124e-6
+
+    return np.where((D > 1e-4) & (D < 300),
+        10**(aP00
+        + aP01 * np.log10(D)**1
+        + aP02 * np.log10(D)**2
+        + aP03 * np.log10(D)**3
+        + aP04 * np.log10(D)**4
+        + aP05 * np.log10(D)**5
+        + aP06 * np.log10(D)**6
+        + aP07 * np.log10(D)**7
+        + aP08 * np.log10(D)**8
+        + aP09 * np.log10(D)**9
+        + aP10 * np.log10(D)**10
+        + aP11 * np.log10(D)**11
+        + aP12 * np.log10(D)**12
+        + aP13 * np.log10(D)**13
+        + aP14 * np.log10(D)**14),float('nan'))
+
+def N1mars(T):
+    # Mars crater SFD
+    # Ivanov (2001) SSR v. 96 pp. 87-104
+    return 2.68e-14 * (np.exp(6.93 * T) - 1) + 4.13e-4 * T
+
+def Nmars(D):
+    aM00 = -3.384
+    aM01 = -3.197
+    aM02 = 1.257
+    aM03 = 0.7915
+    aM04 = -0.4861
+    aM05 = -0.3630
+    aM06 = 0.1016
+    aM07 = 6.756e-2
+    aM08 = -1.181e-2
+    aM09 = -4.753e-3
+    aM10 = 6.233e-4
+    aM11 = 5.805e-5
+
+    return np.where((D > 0.01) & (D < 1000.),
+        10**(aM00
+        + aM01 * np.log10(D)**1
+        + aM02 * np.log10(D)**2
+        + aM03 * np.log10(D)**3
+        + aM04 * np.log10(D)**4
+        + aM05 * np.log10(D)**5
+        + aM06 * np.log10(D)**6
+        + aM07 * np.log10(D)**7
+        + aM08 * np.log10(D)**8
+        + aM09 * np.log10(D)**9
+        + aM10 * np.log10(D)**10
+        + aM11 * np.log10(D)**11),
+        np.full_like(D, np.nan))

python/readme.txt → python/ctem/readme.txt

@@ -1,22 +1,22 @@
-1)	Install a distribution of Python 2.x to your machine. This source code
-	has been tested on Anaconda Python 4.1.1.
-
-2)  Insure required modules are installed:
-	numpy
-	os
-	scipy
-	shutil
-	subprocess
-	
-3)	Place the .py files in this directory into the directory you wish
-	to execute from, with the required initialization files.  For this,
-	the IDL directory may be copied, and ctem.in replaced by the one found
-	in this directory.
-	
-4)  Inside Anaconda Python, the program can be run through the Spyder2 GUI.
-	Outside of Anaconda, the progam can be begun by configuring the ctem.in
-	file, then typing
-	
-	python < ctem_driver.py
-	
-	at the command line.
+1)	Install a distribution of Python 2.x to your machine. This source code
+	has been tested on Anaconda Python 4.1.1.
+
+2)  Insure required modules are installed:
+	numpy
+	os
+	scipy
+	shutil
+	subprocess
+
+3)	Place the .py files in this directory into the directory you wish
+	to execute from, with the required initialization files.  For this,
+	the IDL directory may be copied, and ctem.in replaced by the one found
+	in this directory.
+
+4)  Inside Anaconda Python, the program can be run through the Spyder2 GUI.
+	Outside of Anaconda, the progam can be begun by configuring the ctem.in
+	file, then typing
+
+	python < ctem_driver.py
+
+	at the command line.

src/crater/crater_dimensions.f90

@@ -68,8 +68,8 @@ subroutine crater_dimensions(user,crater,domain)
 
    ! Calculate the radius where the inner wall meets the original pre-existing surface
    ! This is used to demark the location where excavation transitions to deposition
-   crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad)
    crater%ejrim = 0.14_DP * (crater%fcrat * 0.5_DP)**(0.74_DP) ! McGetchin et al. (1973) Thickness of ejecta at rim
+   crater%ejrad = max(crater_profile_find_r_inner_wall(user,crater) * crater%frad, crater%rad)
 
    !find rim for counting purposes
    crater%frim = RIMFAC * crater%frad

src/crater/crater_realistic_topography.f90

@@ -314,10 +314,10 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot)
    integer(I4B), parameter :: num_octaves  = 10   ! Number of Perlin noise octaves
    integer(I4B), parameter :: offset = 4000 ! Scales the random xy-offset so that each crater's random noise is unique 
    real(DP), parameter :: xy_noise_fac = 3.0_DP  ! Spatial "size" of noise features at the first octave
-   real(DP), parameter :: noise_height = 7.0e-3_DP  ! Spatial "size" of noise features at the first octave
+   real(DP), parameter :: noise_height = 3.0e-3_DP  ! Spatial "size" of noise features at the first octave
    real(DP), parameter :: freq = 2.0_DP     ! Spatial size scale factor multiplier at each octave level
    real(DP), parameter :: pers = 1.20_DP  ! The relative size scaling at each octave level
-   real(DP), parameter :: outer_wall_size = 2.1_DP
+   real(DP), parameter :: outer_wall_size = 1.2_DP
 
    !Executable code
    call random_number(rn)
@@ -344,7 +344,7 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot)
          areafrac = 1.0 - util_area_intersection(0.3_DP * crater%floordiam,xbar,ybar,user%pix)
          areafrac = areafrac * util_area_intersection(outer_wall_size * crater%frad,xbar,ybar,user%pix)
          areafrac = areafrac * min((r / flr)**12,1.0_DP) ! Smooth out interface between wall and floor
-         areafrac = areafrac * max(min(2._DP - r,1.0_DP),0.0_DP) ! Smooth out region outside of the rim
+         areafrac = areafrac * max(min(outer_wall_size- r,1.0_DP),0.0_DP) ! Smooth out region outside of the rim
 
          ! Add some roughness to the walls
          noise = 0.0_DP
@@ -356,7 +356,6 @@ subroutine complex_wall_texture(user,surf,crater,domain,deltaMtot)
          end do
          newdem = newdem + noise * areafrac
          if (r < flr) newdem = max(newdem,crater%melev - crater%floordepth)
-         if (r > 1.1_DP) newdem = max(newdem,newdem + areafrac * crater%ejrim * r**(-3))
 
          elchange  = newdem - surf(xpi,ypi)%dem
          deltaMtot = deltaMtot + elchange
@@ -405,6 +404,7 @@ subroutine complex_terrace(user,surf,crater,deltaMtot)
    real(DP)                      :: router          ! The radius of the terrace outer edge
    real(DP)                      :: rinner          ! The radius of the terrace outer edge
    real(DP)                      :: upshift,dfloor
+   real(DP)                      :: h_scallop_profile ! Profile of the scallop wall
    integer(I4B)                  :: num_oct_tfloor 
    real(DP)                      :: noise_height_tfloor
    real(DP)                      :: freq_tfloor
@@ -424,7 +424,9 @@ subroutine complex_terrace(user,surf,crater,deltaMtot)
    nscallops = 16 
    scallop_p = 0.6_DP 
    scallop_width = 0.20_DP 
-   rimfloor = crater%melev
+   h_scallop_profile = 2.0_DP
+   rimfloor = crater%melev + 0.5_DP * crater%rimheight
+   upshift = 0.0_DP
 
    num_oct_tfloor  = 4
    noise_height_tfloor = crater%floordepth / nterraces
@@ -441,13 +443,12 @@ subroutine complex_terrace(user,surf,crater,deltaMtot)
    !^^^^^^^^^^^^^^^
 
    rad = 2.0_DP * crater%frad
-   upshift = 0._DP
 
    inc = max(min(nint(rad / user%pix),PBCLIM*user%gridsize),1) + 1
    crater%maxinc = max(crater%maxinc,inc)
 
    flr = crater%floordiam / crater%fcrat
-   do terrace = 1, nterraces 
+   do terrace = 1, nterraces
       router = flr + (1._DP - flr) * (terrace / real(nterraces,kind=DP))**(terracefac) ! The radius of the outer edge of the terrace
       rinner = flr + (1._DP - flr) * ((terrace - 1) / real(nterraces,kind=DP))**(terracefac) ! The radius of the inner edge of the terrace
 
@@ -484,7 +485,11 @@ subroutine complex_terrace(user,surf,crater,deltaMtot)
 
             isterrace = 1.0_DP - noise < hprof
 
-            tprof = crater_profile(user,crater,rinner) * (r/rinner)**2 + tnoise ! This is the floor profile that replaces the old one at each terrace
+            if (r < 1.0_DP) then
+               tprof = crater_profile(user,crater,rinner) * (r/rinner)**h_scallop_profile + tnoise ! This is the floor profile that replaces the old one at each terrace
+            else
+               tprof = (crater_profile(user,crater,rinner) + crater%ejrim * (rinner)**(-EJPROFILE)) + tnoise ! This is the floor profile that replaces the old one at each terrace
+            end if
 
             if (isterrace) then
                newdem = min(tprof,newdem) 
@@ -517,7 +522,7 @@ subroutine complex_terrace(user,surf,crater,deltaMtot)
          r = sqrt(xbar**2 + ybar**2) / crater%frad
 
          ! Make scalloped rim
-         znoise = scallop_width**(1._DP / (2 * scallop_p))
+         znoise = (scallop_width)**(1._DP / (2 * scallop_p))
          xynoise = (nscallops / PI) / crater%fcrat
          dnoise = util_perlin_noise(xynoise * xbar + offset * rn(1), &
                                     xynoise * ybar + offset * rn(2)) * znoise
@@ -563,7 +568,7 @@ subroutine complex_terrace(user,surf,crater,deltaMtot)
          r = sqrt(xbar**2 + ybar**2) / crater%frad
 
          if (r > flr) then
-            newdem = newdem + upshift * max(min(3.0_DP - 2 * r,1.0_DP),0.0_DP)
+            newdem = newdem + upshift * max(min(3.0_DP - 2 * r**2,1.0_DP),0.0_DP)
             elchange  = newdem - surf(xpi,ypi)%dem
             deltaMtot = deltaMtot + elchange
             surf(xpi,ypi)%dem = newdem

src/util/util_diffusion_solver.f90

@@ -22,6 +22,7 @@
 subroutine util_diffusion_solver(user,surf,N,indarray,kdiff,cumulative_elchange,maxhits)
    use module_globals
    use module_util, EXCEPT_THIS_ONE => util_diffusion_solver
+   use, intrinsic :: ieee_exceptions
    implicit none
 
    ! Arguments
@@ -40,7 +41,12 @@ subroutine util_diffusion_solver(user,surf,N,indarray,kdiff,cumulative_elchange,
    integer(I4B) :: mi,mip1,mim1,mj,mjp1,mjm1,nloops
    real(DP),dimension(N,N) :: elchange
    integer(I4B),dimension(6,N,N) :: indarray_extended
+   logical, dimension(size(IEEE_ALL))   :: fpe_halting_modes, fpe_quiet_modes
+   logical, dimension(size(IEEE_USUAL)) :: fpe_flag 
 
+   call ieee_get_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Save the current halting modes so we can turn them off temporarily
+   fpe_quiet_modes(:) = .false.
+   call ieee_set_halting_mode(IEEE_ALL,fpe_quiet_modes)
 
    cumulative_elchange = 0.0_DP
    kdiffmax = maxval(kdiff)
@@ -122,6 +128,8 @@ subroutine util_diffusion_solver(user,surf,N,indarray,kdiff,cumulative_elchange,
       cumulative_elchange = cumulative_elchange + elchange 
    end do
 
+   call ieee_set_halting_mode(IEEE_ALL,fpe_halting_modes)  ! Restore the original FPE halting modes
+
 return
 end subroutine util_diffusion_solver