Added lunar mare and highlands runs to tests directory

tests/Minton_et_al_2019_ray2/RAYPARAMS.txt

@@ -0,0 +1,8 @@
+rw0 = rmin * pi / Nraymax / 2
+ireal(DP),parameter :: rray = 16_DP
+integer(I4B),parameter :: Nraymax = 16
+real(DP),parameter :: fpeak = 8000_DP
+real(DP),parameter :: rayp = 2.0_DP 
+integer(I4B),parameter :: rayq = 4
+real(DP),parameter :: rayfmult  = (5)**(-4.0_DP / (1.2_DP))  
+

tests/Minton_et_al_2019_ray2/Robbins2014-A15-FassettCounts-binned.dat

@@ -0,0 +1,13 @@
+   5.1241684784040373       6.23948354994535758E-004   21.000000000000000     
+   6.9159685663443886       1.94218692856060401E-002   376.00000000000000     
+   9.1971717131274460       5.97003450726182050E-002   695.00000000000000     
+   12.816081239799573       6.62346057220059986E-002   403.00000000000000     
+   17.996125645322635       5.46997230809128931E-002   170.00000000000000     
+   25.800613087626900       6.43643621282638384E-002   96.000000000000000     
+   36.348639964664329       5.70035776365636210E-002   43.000000000000000     
+   50.608244555146413       5.81896598544816915E-002   23.000000000000000     
+   72.473302918894660       8.93145572505870955E-002   17.000000000000000     
+   101.84498213506704       8.24770394157215703E-002   8.0000000000000000     
+   147.19614402527847       0.13205344016811349        6.0000000000000000     
+   210.22410381342851        0.0000000000000000        0.0000000000000000     
+   262.53658899999999       6.24379576535417230E-002   1.0000000000000000     

tests/Minton_et_al_2019_ray2/cleanup

@@ -0,0 +1,32 @@
+#!/bin/bash
+rm console/console*.*
+rm surf/surf*.*
+rm rego/rego*.*
+rm shaded/shaded*.*
+rm dist/odist_*.dat
+rm dist/tdist_*.dat
+rm dist/pdist_*.dat
+rm dist/ocum_*.dat
+rm dist/tcum_*.dat
+rm misc/mass_*.dat
+rm misc/ctem_*.dat
+rm presentation/presentation*.jpg
+rm surface_*.dat
+rm odistribution.dat
+rm pdistribution.dat
+rm tdistribution.dat
+rm tcumulative.dat
+rm ocumulative.dat
+rm craterscale.dat
+rm ctem.dat
+rm fracdone.dat
+rm production.dat
+rm ejecta_table_max.dat
+rm ejecta_table_min.dat
+rm impactmass.dat
+rm frames/*
+rm regolithdepth.dat
+rm avgkdiff.dat
+rm *.o*
+rm term.out
+rm time.out

tests/Minton_et_al_2019_ray2/ctem.in

@@ -0,0 +1,84 @@
+! CTEM Input file: Compatible with v. 1.1
+
+! IDL driver in uts
+interval     0.899606 
+numintervals 300                                 ! Total number of intervals
+restart      F                                  ! Restart a previous run
+impfile      eta-3.2.dat          ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1))
+popupconsole F                                  ! Pop up console window every output interval 
+saveshaded   F                                  ! Output shaded relief images  
+saverego     F                                  ! Output regolith map images 
+savepres     F                                  ! Output simplified console display images (presentation-compatible images)
+savetruelist F                                  ! Save the true cumulative crater distribution for each interval (large file size)
+sfdcompare   Robbins2014-A15-FassettCounts-binned.dat         ! File name for the SFD comparison file used in the console display
+shadedminh   -85.0                              ! Minimum height for shaded relief map (m) (Default - automatic)
+shadedmaxh   85.0                               ! Maximum height for shaded relief map (m) (Default - automatic)
+runtype      single                             ! Run type: options are normal / statistical
+                                                !    single:      craters accumulate in successive intervals
+                                                !    statistical: surface is reset between intervals
+
+! CTEM required inputs
+seed         21286                    ! Random number generator seed
+gridsize     1000                          ! Size of grid in pixels
+numlayers    10                                 ! Number of perched layers
+pix          1.0                              ! Pixel size (m)
+mat          rock                               ! Material (rock or ice)
+! Bedrock scaling parameters
+mu_b         0.55d0                             ! Experimentally derived parameter for bedrock crater scaling law 
+kv_b         0.20d0                             ! Experimentally derived parameter for bedrock crater scaling law 
+trho_b       2250.0d0                           ! Target density (bedrock) (kg/m**3)
+ybar_b       0.0d6                               ! Bedrock strength (Pa)
+! Regolith scaling parameters
+mu_r         0.55d0                             ! Experimentally derived parameter for regolith crater scaling law 
+kv_r         0.20d0                             ! Experimentally derived parameter for regolith crater scaling law
+trho_r       2250.0d0                           ! Target density (regolith) (kg/m**3)
+ybar_r       0.00d6                              ! Regolith strength (Pa)
+! Body parameters
+gaccel       1.62d0                             ! Gravitational acceleration at target (m/s**2)
+trad         1737.35d3                          ! Target radius (m)
+prho         2500.0d0                           ! Projectile density (kg/m**3)
+sfdfile      production.dat                     ! Impactor SFD file (col 1: Dimp (m), col 2: ! impactors > D
+velfile      lunar-MBA-impactor-velocities.dat  ! Impactor velocity dist file
+
+! Seismic shaking input (required if seismic shaking is set to T, otherwise ignored)
+doseismic    F                                  ! Perform seismic shaking calculations with each impact - Default F
+seisq        2000.0d0                           ! Seimsic energy attenuation quality factor, Q 
+                                                ! Value taken from Richardson et al. (2005)
+
+tvel         100.d0                             ! Target p-wave speed (m/s)
+                                                ! Value taken from Watkins and Kovach (1973) LPSC 
+
+tfrac         2000.d0                             ! Mean free path for seismic wave scattering in target (m)
+                                                ! Value taken from Dainty et al. (1974) The Moon v. 9 p. 11
+
+regcoh       520.d0                             ! Regolith layer cohesion (Pa)
+                                                ! Value taken from Colwell et al. (2007) Rev. of Geophys. v. 45 p. RG2006
+                                                ! For upper 0-60 cm
+
+neff         5.0d-10                            ! Impact energy seismic efficiency factor 
+
+
+! Optional inputF These have internally set default values that work reasonable well. Comment them out with
+deplimit        10d3                            ! Depth limit for craters (m) - Default is to ignore.
+maxcrat         0.33d0                          ! Fraction of gridsize that maximum crater can be - Default 1.0
+killatmaxcrater F                               ! Stop the run if a crater larger than the maximum is produced - Default F      
+basinimp        35.0d3                          ! Size of impactor to switch to lunar basin scaling law - Default is to ignore
+docollapse      T                               ! Do slope collapse - Default T
+dosoftening     T                               ! Do ejecta softening - Default T
+doangle         T                               ! Vary the impact angle. Set to F to have only vertical impacts - Default T
+Kd1             0.00312669649143281
+psi             2.000
+fe              10.0
+ejecta_truncation 20.0
+dorays  T 
+
+! Testing input. These are used to perform non-Monte Carlo tests.
+testflag     F                                 ! Set to T to create a single crater with user-defined impactor properties
+testimp      20.0                              ! Diameter of test impactor (m)
+testvel      15.d3                             ! Velocity of test crater (m/s)
+testang      90.0                              ! Impact angle of test crater (deg) - Default 90.0
+testxoffset  0.0d0                             ! x-axis offset of crater center from grid center (m) - Default 0.0
+testyoffset  0.0d0                             ! y-axis offset of crater center from grid center (m) - Default 0.0
+tallyonly    F                                 ! Tally the craters without generating any craters
+testtally    F                                 ! Set to T to count all non-cookie cut craters, regardless of score
+

tests/Minton_et_al_2019_ray2/ctem_driver.pro

@@ -0,0 +1,233 @@
+pro ctem_driver
+
+;-------------------------------------------------------------------------
+;  Jim Richardson, Arecibo Observatory
+;  David Minton, Purdue University Dept. of Earth, Atmospheric, & Planetary Sciences
+;  July 2014
+;  Cratered Terrain Evolution Model display module
+;
+;  Inputs are read in from the ctem.in file
+;
+;-------------------------------------------------------------------------
+;------------- Initial Setup ----------------
+;-------------------------------------------------------------------------
+Compile_Opt DEFINT32
+!EXCEPT=2
+
+; ----------- input file -----------------------
+infilename = 'ctem.in'
+DATFILE='ctem.dat'
+
+; ---------- reading input files ----------
+seedarr = lon64arr(100)
+seedn = 1
+totalimpacts = long64(0)
+ncount = long64(0)
+curyear = 0.d0
+restart = "F" 
+fracdone = 1.0d0
+masstot = 0.d0
+
+ctem_io_read_input,infilename,interval,numintervals,gridsize,pix,seed,numlayers,sfdfile,impfile,maxcrat,ph1,shadedmaxhdefault,shadedminhdefault,shadedminh,shadedmaxh,restart,runtype,popupconsole,saveshaded,saverego,savepres,savetruelist
+
+seedarr(0) = seed
+area = (gridsize * pix)^2
+
+;read input data
+pnum = file_lines(impfile)
+production = dblarr(2,pnum)
+productionfunction = dblarr(2,pnum)
+openr,LUN,impfile,/GET_LUN
+readf,LUN,productionfunction
+free_lun,LUN
+
+;create impactor production population
+production(0,*) = productionfunction(0,*)
+production(1,*) = productionfunction(1,*)*area*interval
+
+;write out corrected production population
+openw,1,sfdfile
+printf,1,production
+close,1
+free_lun,1
+
+;set up cratering surface grid and display-only grid
+surface_dem = dblarr(gridsize,gridsize)
+regolith = dblarr(gridsize,gridsize)
+
+;set up temporary distribution bins
+distl = 1
+pdistl = 1
+odist = dblarr(6,distl)
+tdist = dblarr(6,distl)
+pdist = dblarr(6,pdistl)
+pdisttotal = dblarr(6,pdistl)
+
+datformat = "(I17,1X,I12,1X,E19.12,1X,A1,1X,F9.6,1X,E19.12)"
+
+
+if strmatch(restart,'F',/fold_case) then begin ; Start with a clean slate
+	print, 'Starting a new run'
+	curyear = 0.0d0
+	totalimpacts = 0
+	masstot = 0.d0
+	fracdone = 1.0d0
+
+	if strmatch(runtype,'statistical',/fold_case) then begin
+		ncount = 1
+		openw,LUN,DATFILE,/GET_LUN
+		printf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat
+		for n=0,seedn-1 do begin
+			printf,LUN,seedarr(n),format='(I12)'
+		endfor
+		free_lun,LUN
+	endif else begin
+		ncount = 0
+	endelse
+
+	surface_dem(*,*) = 0.0d0
+	regolith(*,*) = 0.0d0
+
+	file_delete, 'tdistribution.dat',/allow_nonexistent
+
+endif else begin ; continue an old run
+	print, 'Continuing a previous run'
+
+	ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass
+
+	;read in constants file
+	openr,LUN,DATFILE,/GET_LUN
+	readf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat
+	seedn = 0
+	while ~ eof(LUN) do begin 
+		readf,LUN,iseed
+		seedarr(seedn) = long64(iseed)
+	   seedn=seedn+1
+	endwhile	
+
+endelse
+
+openw,FRACDONEFILE,'fracdone.dat',/GET_LUN
+openw,REGODEPTHFILE,'regolithdepth.dat',/GET_LUN
+
+;-------------------------------------------------------------------------
+; ---------- begin loops ----------
+;-------------------------------------------------------------------------
+print, 'Beginning loops'
+
+;define number of loop iterations and begin
+;numintervals=ceil(endyear/interval)-ncount
+while (ncount le numintervals) do begin
+
+
+	; ---------- creating crater population ----------
+	if (ncount gt 0) then begin
+
+		fnum = string(ncount,format='(I6.6)')
+		if (file_test('misc',/DIRECTORY) eq 0) then begin
+			file_mkdir,'misc'
+		endif
+		; save a copy of the ctem.dat file 
+		fname = 'misc/ctem_' + fnum + '.dat'
+		file_copy, 'ctem.dat', fname, /OVERWRITE
+
+		print, ncount, '  Calling FORTRAN routine'
+		;call fortran program to create & count craters
+		spawn, './CTEM',/noshell
+
+		; ---------- reading FORTRAN output ----------
+		print, ncount, '  Reading FORTRAN output'
+		ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass
+
+		;read in constants file
+		openr,LUN,DATFILE,/GET_LUN
+		readf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat
+		seedn = 0
+		while ~ eof(LUN) do begin 
+			readf,LUN,iseed
+			seedarr(seedn) = long64(iseed)
+			seedn = seedn + 1
+		endwhile	
+		free_lun,LUN
+
+		curyear = curyear + fracdone * interval
+		masstot = masstot + mass
+		printf,FRACDONEFILE,fracdone,curyear
+		flush,FRACDONEFILE
+
+		printf,REGODEPTHFILE,curyear,mean(regolith),max(regolith),min(regolith)
+		flush,REGODEPTHFILE
+
+		;save a copy of the binned observed crater distribution
+		if (file_test('dist',/DIRECTORY) eq 0) then begin
+			file_mkdir,'dist'
+		endif
+		fname = 'dist/odist_' + fnum + '.dat'
+		file_copy, 'odistribution.dat', fname, /OVERWRITE
+
+		; save a copy of the cumulative observed crater distribution 
+		fname = 'dist/ocum_' + fnum + '.dat'
+		file_copy, 'ocumulative.dat', fname, /OVERWRITE
+
+		;save a copy of the binned true distribution
+		fname = 'dist/tdist_' + fnum + '.dat'
+		file_copy, 'tdistribution.dat', fname, /OVERWRITE
+
+		;save a copy of the binned idealized production function
+		fname = 'dist/pdist_' + fnum + '.dat'
+		file_copy, 'pdistribution.dat', fname, /OVERWRITE
+
+		; save a copy of the cumulative true crater distribution if the user requests it
+		if strmatch(savetruelist,'T',/fold_case) then begin
+			fname = 'dist/tcum_' + fnum + '.dat'
+			file_copy, 'tcumulative.dat', fname, /OVERWRITE
+		endif
+
+		; save a copy of the impacted mass 
+		fname = 'misc/mass_' + fnum + '.dat'
+		file_copy, 'impactmass.dat', fname, /OVERWRITE
+
+
+	endif
+
+	; Get the accumulated production function
+	pdisttotal = pdist 
+	pdisttotal(3:5,*) = pdist(3:5,*) * curyear / interval
+
+	; ---------- displaying results ----------
+	print, ncount, '  Displaying results'
+
+	ctem_image_dem,ncount,gridsize,pix,surface_dem,surface_dem_image
+	if strmatch(saverego,'T',/fold_case) then ctem_image_regolith,ncount,gridsize,pix,regolith,regolith_image
+	if strmatch(saveshaded,'T',/fold_case) then begin
+		if (shadedminhdefault eq 1) then shadedminh = min(surface_dem)
+		if (shadedmaxhdefault eq 1) then shadedmaxh = max(surface_dem)
+		ctem_image_shaded_relief,ncount,gridsize,pix,surface_dem,surface_dem,shadedminh,shadedmaxh,'shaded',shaded_image
+	endif
+	if strmatch(savepres,'T',/fold_case) then ctem_image_presentation,ncount,gridsize,pix,curyear,odist,pdisttotal,tdist,ph1,surface_dem_image
+	ctem_window_display,ncount,totalimpacts,gridsize,pix,curyear,masstot,odist,pdisttotal,tdist,ph1,surface_dem,regolith,surface_dem_image,popupconsole
+
+	ncount = ncount + 1
+
+	;write out the current data file
+	if (strmatch(runtype,'statistical',/fold_case)) || (ncount eq 1) then begin 
+		restart = 'F'
+		curyear = 0.0d0
+		totalimpacts = 0
+		masstot = 0.d0
+		file_delete, 'tdistribution.dat',/allow_nonexistent
+	endif else begin
+		restart = 'T'
+	endelse
+
+	openw,LUN,DATFILE,/GET_LUN
+	printf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat
+	for n=0,seedn-1 do begin
+		printf,LUN,seedarr(n),format='(I12)'
+	endfor
+	free_lun,LUN
+endwhile
+free_lun,FRACDONEFILE
+free_lun,REGODEPTHFILE
+
+end