diff --git a/.gitignore b/.gitignore index e16dce6f..70c42081 100644 --- a/.gitignore +++ b/.gitignore @@ -7,7 +7,10 @@ !distclean.cmake !README.md !version.txt +!ctem/ !ctem/**.py !cmake/Modules/*.cmake +!*.bash + !pyproject.toml !LICENSE \ No newline at end of file diff --git a/CMakeLists.txt b/CMakeLists.txt index eb4fdb32..f6cd01ba 100644 --- a/CMakeLists.txt +++ b/CMakeLists.txt @@ -11,11 +11,12 @@ ################################################## # Define the project and the depencies that it has ################################################## -CMAKE_MINIMUM_REQUIRED(VERSION 3.6.0...3.27.1) +CMAKE_MINIMUM_REQUIRED(VERSION 3.24.0...3.27.1) +SET(SKBUILD_PROJECT_NAME "ctem" CACHE STRING "Name of project set by scikit-build") # Get version stored in text file FILE(READ "version.txt" VERSION) -PROJECT(CTEM LANGUAGES C Fortran VERSION ${VERSION}) +PROJECT(${SKBUILD_PROJECT_NAME} LANGUAGES C Fortran VERSION ${VERSION}) IF (CMAKE_Fortran_COMPILER_ID MATCHES "^Intel") SET(COMPILER_OPTIONS "Intel" CACHE STRING "Compiler identified as Intel") @@ -26,7 +27,6 @@ ELSE () ENDIF () # Set some options the user may choose -OPTION(USE_COARRAY "Use Coarray Fortran for parallelization of test particles" OFF) OPTION(USE_OPENMP "Use OpenMP for parallelization" ON) OPTION(USE_SIMD "Use SIMD vectorization" ON) OPTION(BUILD_SHARED_LIBS "Build using shared libraries" ON) @@ -39,30 +39,40 @@ ENDIF() # Ensure scikit-build modules FIND_PACKAGE(Python COMPONENTS Interpreter Development.Module REQUIRED) -# Add our local modules to the module path -FILE(TO_CMAKE_PATH "${CMAKE_SOURCE_DIR}/cmake/Modules" LOCAL_MODULE_PATH) -LIST(APPEND CMAKE_MODULE_PATH ${LOCAL_MODULE_PATH}) - # Define some directories that are important to the build SET(SRC "${CMAKE_SOURCE_DIR}/src") -SET(LIB "${CMAKE_SOURCE_DIR}/lib") -SET(BIN "${CMAKE_SOURCE_DIR}/bin") -SET(MOD "${CMAKE_SOURCE_DIR}/include") SET(PY "${CMAKE_SOURCE_DIR}/ctem") # Make sure paths are correct for Unix or Windows style FILE(TO_CMAKE_PATH ${SRC} SRC) -FILE(TO_CMAKE_PATH ${LIB} LIB) -FILE(TO_CMAKE_PATH ${BIN} BIN) -FILE(TO_CMAKE_PATH ${MOD} MOD) FILE(TO_CMAKE_PATH ${PY} PY) -set(CMAKE_ARCHIVE_OUTPUT_DIRECTORY ${LIB}) -set(CMAKE_LIBRARY_OUTPUT_DIRECTORY ${LIB}) -set(CMAKE_RUNTIME_OUTPUT_DIRECTORY ${BIN}) +INCLUDE(GNUInstallDirs) +IF (SKBUILD) + SET(INSTALL_BINDIR ${SKBUILD_SCRIPTS_DIR}) + SET(INSTALL_LIBDIR ${SKBUILD_DATA_DIR}/lib) + SET(INSTALL_INCLUDEDIR ${SKBUILD_HEADERS_DIR}) + SET(INSTALL_PYPROJ ${SKBUILD_PLATLIB_DIR}/${SKBUILD_PROJECT_NAME}) + IF (APPLE) + SET(CMAKE_INSTALL_RPATH "@loader_path;${CMAKE_BINARY_DIR}/bin") + ELSEIF (LINUX) + SET(CMAKE_INSTALL_RPATH "@ORIGIN;${CMAKE_BINARY_DIR}/bin") + ENDIF () +ELSE () + SET(INSTALL_PYPROJ ${PY}) + SET(INSTALL_BINDIR ${CMAKE_INSTALL_BINDIR}) + SET(INSTALL_LIBDIR ${CMAKE_INSTALL_LIBDIR}) + SET(INSTALL_INCLUDEDIR ${CMAKE_INSTALL_INCLUDEDIR}) +ENDIF () + +SET(CMAKE_INSTALL_RPATH_USE_LINK_PATH TRUE) # Have the .mod files placed in the lib folder -SET(CMAKE_Fortran_MODULE_DIRECTORY ${MOD}) +SET(CMAKE_Fortran_MODULE_DIRECTORY ${CMAKE_CURRENT_BINARY_DIR}/mod) + +# Add our local modules to the module path +FILE(TO_CMAKE_PATH "${CMAKE_SOURCE_DIR}/cmake/Modules" LOCAL_MODULE_PATH) +LIST(APPEND CMAKE_MODULE_PATH ${LOCAL_MODULE_PATH}) # Set the name of the ctem library SET(CTEM_LIBRARY ctem) diff --git a/cmake/Modules/SetFortranFlags.cmake b/cmake/Modules/SetFortranFlags.cmake index edf73442..585e7e3c 100644 --- a/cmake/Modules/SetFortranFlags.cmake +++ b/cmake/Modules/SetFortranFlags.cmake @@ -533,10 +533,10 @@ IF (CMAKE_BUILD_TYPE STREQUAL "RELEASE" OR CMAKE_BUILD_TYPE STREQUAL "PROFILE") ) # Tells the compiler to link to certain libraries in the Intel oneAPI Math Kernel Library (oneMKL). - SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "/Qmkl:cluster" # Intel Windows - "/Qmkl" # Intel Windows - ) + # SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + # Fortran "/Qmkl:cluster" # Intel Windows + # "/Qmkl" # Intel Windows + # ) # Enables additional interprocedural optimizations for a single file compilation SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" @@ -574,10 +574,10 @@ IF (CMAKE_BUILD_TYPE STREQUAL "RELEASE" OR CMAKE_BUILD_TYPE STREQUAL "PROFILE") ) # Tells the compiler to link to certain libraries in the Intel oneAPI Math Kernel Library (oneMKL). - SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-qmkl=cluster" # Intel - "-qmkl" # Intel - ) + # SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" + # Fortran "-qmkl=cluster" # Intel + # "-qmkl" # Intel + # ) # Enables additional interprocedural optimizations for a single file compilation SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_RELEASE "${CMAKE_Fortran_FLAGS_RELEASE}" @@ -671,11 +671,11 @@ IF (CMAKE_BUILD_TYPE STREQUAL "PROFILE") # Enables the optimization reports to be generated IF (WINOPT) SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "/O2 /Qopt-report:5 /traceback /Z7" # Intel Windows + Fortran "/O2 /Qopt-report:3 /traceback /Z7" # Intel Windows ) ELSE () SET_COMPILE_FLAG(CMAKE_Fortran_FLAGS_PROFILE "${CMAKE_Fortran_FLAGS_RELEASE}" - Fortran "-O2 -pg -qopt-report=5 -traceback -p -g3" # Intel + Fortran "-O2 -pg -qopt-report=3 -traceback -p -g3" # Intel ) ENDIF () ELSEIF (COMPILER_OPTIONS STREQUAL "GNU") diff --git a/cmake/Modules/SetParallelizationLibrary.cmake b/cmake/Modules/SetParallelizationLibrary.cmake index 505a77d6..71a382a9 100644 --- a/cmake/Modules/SetParallelizationLibrary.cmake +++ b/cmake/Modules/SetParallelizationLibrary.cmake @@ -17,20 +17,11 @@ IF (USE_OPENMP) ENDIF (NOT OpenMP_Fortran_FLAGS) ENDIF (USE_OPENMP) -IF (USE_COARRAY) - IF (NOT Coarray_Fortran_FLAGS) - FIND_PACKAGE (Coarray_Fortran) - IF (NOT Coarray_Fortran_FLAGS) - MESSAGE (FATAL_ERROR "Fortran compiler does not support Coarrays") - ENDIF (NOT Coarray_Fortran_FLAGS) - ENDIF (NOT Coarray_Fortran_FLAGS) -ENDIF (USE_COARRAY) -IF (NOT USE_OPENMP AND NOT USE_COARRAY) - # Turn off both OpenMP and CAF +IF (NOT USE_OPENMP) + # Turn off OpenMP SET (OMP_NUM_PROCS 0 CACHE STRING "Number of processors OpenMP may use" FORCE) UNSET (OpenMP_Fortran_FLAGS CACHE) - UNSET (Coarray_Fortran_FLAGS CACHE) UNSET (GOMP_Fortran_LINK_FLAGS CACHE) -ENDIF (NOT USE_OPENMP AND NOT USE_COARRAY) +ENDIF () diff --git a/ctem/ctem/NPF.py b/ctem/NPF.py similarity index 100% rename from ctem/ctem/NPF.py rename to ctem/NPF.py diff --git a/ctem/__init__.py b/ctem/__init__.py new file mode 100644 index 00000000..81f6944a --- /dev/null +++ b/ctem/__init__.py @@ -0,0 +1,2 @@ +from ctem.driver import * +from ctem import util \ No newline at end of file diff --git a/ctem/ctem/craterproduction.py b/ctem/craterproduction.py similarity index 100% rename from ctem/ctem/craterproduction.py rename to ctem/craterproduction.py diff --git a/ctem/ctem/__init__.py b/ctem/ctem/__init__.py deleted file mode 100644 index 2e8ce2a0..00000000 --- a/ctem/ctem/__init__.py +++ /dev/null @@ -1 +0,0 @@ -from ctem.driver import * \ No newline at end of file diff --git a/ctem/ctem/driver.py b/ctem/driver.py similarity index 96% rename from ctem/ctem/driver.py rename to ctem/driver.py index 0249c19d..661d18c8 100644 --- a/ctem/ctem/driver.py +++ b/ctem/driver.py @@ -52,14 +52,10 @@ def __init__(self, param_file="ctem.in", isnew=True): } # Get the location of the CTEM executable - self.ctem_executable = Path(_pyfile).parent.parent.parent / "bin" / "CTEM" - if not self.ctem_executable.exists(): - print(f"CTEM driver not found at {self.ctem_executable}. Trying current directory.") - self.ctem_executable = Path(currentdir) / "CTEM" - if not self.ctem_executable.exists(): - warnings.warn(f"Cannot find the CTEM driver {str(self.ctem_executable)}", stacklevel=2) - self.ctem_executable = None - + self.ctem_executable = shutil.which("CTEM") + if not self.ctem_executable: + warnings.warn(f"Cannot find the CTEM driver {str(self.ctem_executable)}", stacklevel=2) + self.ctem_executable = None self.user = util.read_user_input(self.user) @@ -169,7 +165,10 @@ def __init__(self, param_file="ctem.in", isnew=True): rclist[:,0] = np.exp(interp(np.log(rclist[:,0]))) #Convert age in Ga to "interval time" - rclist[:,5] = (self.user['interval'] * self.user['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') + if (self.user['runtype'].upper() == 'STATISTICAL'): + rclist[:,5] = (self.user['interval']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') + else: + rclist[:,5] = (self.user['interval'] * self.user['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') rclist = rclist[rclist[:,5].argsort()] #Export to dat file diff --git a/ctem/setup.py b/ctem/setup.py deleted file mode 100644 index f7003601..00000000 --- a/ctem/setup.py +++ /dev/null @@ -1,6 +0,0 @@ -from setuptools import setup, find_packages - -setup(name='ctem', - version='0.1', - author='David A. Minton', - packages=find_packages()) diff --git a/ctem/ctem/util.py b/ctem/util.py similarity index 99% rename from ctem/ctem/util.py rename to ctem/util.py index 00922421..0597c9a4 100644 --- a/ctem/ctem/util.py +++ b/ctem/util.py @@ -136,7 +136,7 @@ def image_regolith(user, regolith): height = user['gridsize'] / dpi width = height fig = plt.figure(figsize=(width, height), dpi=dpi) - fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower') + fig.figimage(regolith_scaled, cmap=cm.magma, origin='lower') plt.savefig(filename) plt.close() diff --git a/ctem/ctem/viewer3d.py b/ctem/viewer3d.py similarity index 100% rename from ctem/ctem/viewer3d.py rename to ctem/viewer3d.py diff --git a/debug_compile.bash b/debug_compile.bash new file mode 100644 index 00000000..05eb2280 --- /dev/null +++ b/debug_compile.bash @@ -0,0 +1,3 @@ +cmake -P distclean.cmake +cmake -B build -S . -DCMAKE_BUILD_TYPE=Debug +cmake --build build diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in index 7b6750a7..97925a0e 100644 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -2,21 +2,21 @@ ! 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 100000 ! Diameter of test impactor (m) +testflag T ! Set to T to create a single crater with user-defined impactor properties +testimp 15000.0 ! Diameter of test impactor (m) testvel 18.3e3 ! Velocity of test crater (m/s) -testang 45.0 ! Impact angle of test crater (deg) - Default 90.0 -testxoffset 0 ! x-axis offset of crater center from grid center (m) - Default 0.0 -testyoffset 0 ! y-axis offset of crater center from grid center (m) - Default 0.0 +testang 90.0 ! Impact angle of test crater (deg) - Default 90.0 +testxoffset 0.0 ! x-axis offset of crater center from grid center (m) - Default 0.0 +testyoffset 0.0 ! y-axis offset of crater center from grid center (m) - Default 0.0 tallyonly F ! Tally the craters without generating any craters testtally F quasimc F ! MC run constrained by non-MC 'real' craters given in a list -realcraterlist qmctest2.in ! list of 'real' craters for Quasi-MC runs +realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs ! IDL driver in uts -interval 10.0 +interval 1.0 numintervals 1 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single' restart F ! Restart a previous run impfile NPFextrap.dat ! Impactor SFD rate file (col 1: Dimp (m), col 2: ! impactors > D (m**(-2) y**(-1)) @@ -33,10 +33,10 @@ runtype single ! Run type: options are normal / ! statistical: surface is reset between intervals ! CTEM required inputs -seed 765 ! Random number generator seed -gridsize 500 ! Size of grid in pixels +seed 1111 ! Random number generator seed +gridsize 1000 ! Size of grid in pixels numlayers 10 ! Number of perched layers -pix 12.32e3 ! Pixel size (m) +pix 6.16e3 ! Pixel size (m) mat rock ! Material (rock or ice) ! Bedrock scaling parameters mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law @@ -60,11 +60,11 @@ doseismic F ! Perform seismic shaking calcul ! Optional inputF These have internally set default values that work reasonable well. Comment them out with deplimit 9e99 ! Depth limit for craters (m) - Default is to ignore. -maxcrat 0.15 ! Fraction of gridsize that maximum crater can be - Default 1.0 (0.15 for full MC; 3.15e-2 for Quasi-MC) +maxcrat 3.15e-2 ! Fraction of gridsize that maximum crater can be - Default 1.0 <--0.15 for full MC; 3.15e-2 for Quasi-MC killatmaxcrater F ! Stop the run if a crater larger than the maximum is produced - Default F -basinimp 35.0e3 ! Size of impactor to switch to lunar basin scaling law - Default is to ignore +basinimp 35.0e6 ! 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 +dosoftening F ! Do ejecta softening - Default T doangle T ! Vary the impact angle. Set to F to have only vertical impacts - Default T Kd1 0.0001 psi 2.000 @@ -74,5 +74,5 @@ doregotrack T dorays F superdomain F dorealistic F -domixing T -dotopodiffusion F \ No newline at end of file +dotopodiffusion F +domixing F \ No newline at end of file diff --git a/examples/global-lunar-bombardment/ctemdriver.py b/examples/global-lunar-bombardment/ctemdriver.py deleted file mode 100644 index 50392285..00000000 --- a/examples/global-lunar-bombardment/ctemdriver.py +++ /dev/null @@ -1,3 +0,0 @@ -import ctem -sim = ctem.Simulation() -sim.run() \ No newline at end of file diff --git a/examples/global-lunar-bombardment/temp.in b/examples/global-lunar-bombardment/temp.in deleted file mode 100644 index 328c21ed..00000000 --- a/examples/global-lunar-bombardment/temp.in +++ /dev/null @@ -1,76 +0,0 @@ -! CTEM Input file - - -! Testing input. These are used to perform non-Monte Carlo tests. -testflag T! F ! Set to T to create a single crater with user-defined impactor properties -testimp 10 ! 15000.0 ! Diameter of test impactor (m) -testvel 18.3e3 ! Velocity of test crater (m/s) -testang 45.0 ! Impact angle of test crater (deg) - Default 90.0 -testxoffset 0.0 ! x-axis offset of crater center from grid center (m) - Default 0.0 -testyoffset 0.0 ! y-axis offset of crater center from grid center (m) - Default 0.0 -tallyonly F ! Tally the craters without generating any craters -testtally F -quasimc F! T ! MC run constrained by non-MC 'real' craters given in a list -realcraterlist qmctest.in ! list of 'real' craters for Quasi-MC runs - - - -! IDL driver in uts -interval 1 ! 0.1 -numintervals 1 ! 1 ! Total number of interval 1 !s (total time = interval 1 ! * numintervals 1 !) <--when runtype is 'single' -restart F ! Restart a previous run -impfile NPFextrap.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 1 ! -saveshaded F ! Output shaded relief images -saverego T ! Output regolith map images -savepres F ! Output simplified console display images (presentation-compatible images) -savetruelist T ! Save the true cumulative crater distribution for each interval 1 ! (large file size) -sfdcompare LOLASethCraterCatalogv8gt20-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 interval 1 !s - ! statistical: surface is reset between interval 1 !s - -! CTEM required inputs -seed 76535 ! Random number generator seed -gridsize 2000 ! Size of grid in pixels -numlayers 10 ! Number of perched layers -pix 3.08e3 ! Pixel size (m) -mat rock ! Material (rock or ice) -! Bedrock scaling parameters -mu_b 0.55e0 ! Experimentally derived parameter for bedrock crater scaling law -kv_b 0.20e0 ! Experimentally derived parameter for bedrock crater scaling law -trho_b 2250.0e0 ! Target density (bedrock) (kg/m**3) -ybar_b 0.0e6 ! Bedrock strength (Pa) -! Regolith scaling parameters -mu_r 0.55e0 ! Experimentally derived parameter for regolith crater scaling law -kv_r 0.20e0 ! Experimentally derived parameter for regolith crater scaling law -trho_r 2250.0e0 ! Target density (regolith) (kg/m**3) -ybar_r 0.00e6 ! Regolith strength (Pa) -! Body parameters -gaccel 1.62e0 ! Gravitational acceleration at target (m/s**2) -trad 1737.35e3 ! Target radius (m) -prho 2500.0e0 ! 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 - -! Optional inputF These have internally set default values that work reasonable well. Comment them out with -deplimit 9e99 ! Depth limit for craters (m) - Default is to ignore. -maxcrat 3.15e-2 ! 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.0e3 ! 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.0001 -psi 2.000 -fe 4.00 -ejecta_truncation 4.0 -doregotrack T -dorays F -superdomain F -dorealistic F diff --git a/examples/quasimc-global/CTEM b/examples/quasimc-global/CTEM deleted file mode 120000 index 604970bd..00000000 --- a/examples/quasimc-global/CTEM +++ /dev/null @@ -1 +0,0 @@ -../../build/src/CTEM \ No newline at end of file diff --git a/examples/ray-generation-test-environment/ctem.in b/examples/ray-generation-test-environment/ctem.in index e9f20086..c297424a 100644 --- a/examples/ray-generation-test-environment/ctem.in +++ b/examples/ray-generation-test-environment/ctem.in @@ -20,9 +20,9 @@ runtype single ! Run type: options are normal / ! CTEM required inputs seed 1108267 ! Random number generator seed -gridsize 1000 ! Size of grid in pixels +gridsize 2000 ! Size of grid in pixels numlayers 10 ! Number of perched layers -pix 1.0e0 ! Pixel size (m) +pix 10.0 ! Pixel size (m) mat rock ! Material (rock or ice) ! Bedrock scaling parameters mu_b 0.41d0 ! Experimentally derived parameter for bedrock crater scaling law @@ -62,18 +62,23 @@ neff 5.0d-10 ! Impact energy seismic efficien ! Optional input. These have internally set default values that work reasonable well. Comment them out with countingmodel Fassett ! Crater counter model. Default is by Caleb Fassett. deplimit 10d3 ! Depth limit for craters (m) - Default is to ignore. -maxcrat 1.00d0 ! Fraction of gridsize that maximum crater can be - Default 1.0 +maxcrat 1.00d4 ! 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 -doregotrack F +doregotrack T +Kd1 0.00312669649143281 +psi 2.000 +fe 10.0 +dorays T +ejecta_truncation 60.0 ! Testing input. These are used to perform non-Monte Carlo tests. testflag T ! Set to T to create a single crater with user-defined impactor properties -testimp 1.0 ! Diameter of test impactor (m) -testvel 15.d3 ! Velocity of test crater (m/s) +testimp 15.0 ! Diameter of test impactor (m) +testvel 18.3e3 ! 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 diff --git a/examples/ray-generation-test-environment/testprofile.dat b/examples/ray-generation-test-environment/testprofile.dat deleted file mode 100644 index 09f1b660..00000000 --- a/examples/ray-generation-test-environment/testprofile.dat +++ /dev/null @@ -1,1000 +0,0 @@ - 1.0000000000000000 1.2525955388259645E-003 - 2.0000000000000000 1.2527493144024800E-003 - 3.0000000000000000 1.2530056375141877E-003 - 4.0000000000000000 1.2533645538970216E-003 - 5.0000000000000000 1.2538261276006199E-003 - 6.0000000000000000 1.2543904410076655E-003 - 7.0000000000000000 1.2550575948587633E-003 - 8.0000000000000000 1.2558277082828800E-003 - 9.0000000000000000 1.2567009188333385E-003 - 10.000000000000000 1.2576773825294040E-003 - 11.000000000000000 1.2587572739034647E-003 - 12.000000000000000 1.2599407860538437E-003 - 13.000000000000000 1.2612281307032500E-003 - 14.000000000000000 1.2626195382629108E-003 - 15.000000000000000 1.2641152579024033E-003 - 16.000000000000000 1.2657155576252257E-003 - 17.000000000000000 1.2674207243501341E-003 - 18.000000000000000 1.2692310639982928E-003 - 19.000000000000000 1.2711469015862706E-003 - 20.000000000000000 1.2731685813249238E-003 - 21.000000000000000 1.2752964667242241E-003 - 22.000000000000000 1.2775309407040661E-003 - 23.000000000000000 1.2798724057111101E-003 - 24.000000000000000 1.2823212838417191E-003 - 25.000000000000000 1.2848780169710419E-003 - 26.000000000000000 1.2875430668883000E-003 - 27.000000000000000 1.2903169154383515E-003 - 28.000000000000000 1.2932000646695841E-003 - 29.000000000000000 1.2961930369882163E-003 - 30.000000000000000 1.2992963753190782E-003 - 31.000000000000000 1.3025106432729382E-003 - 32.000000000000000 1.3058364253204667E-003 - 33.000000000000000 1.3092743269729082E-003 - 34.000000000000000 1.3128249749695513E-003 - 35.000000000000000 1.3164890174720873E-003 - 36.000000000000000 1.3202671242659397E-003 - 37.000000000000000 1.3241599869686746E-003 - 38.000000000000000 1.3281683192455729E-003 - 39.000000000000000 1.3322928570324861E-003 - 40.000000000000000 1.3365343587660638E-003 - 41.000000000000000 1.3408936056214739E-003 - 42.000000000000000 1.3453714017577245E-003 - 43.000000000000000 1.3499685745707120E-003 - 44.000000000000000 1.3546859749540999E-003 - 45.000000000000000 1.3595244775681770E-003 - 46.000000000000000 1.3644849811168085E-003 - 47.000000000000000 1.3695684086326135E-003 - 48.000000000000000 1.3747757077705212E-003 - 49.000000000000000 1.3801078511098319E-003 - 50.000000000000000 1.3855658364649425E-003 - 51.000000000000000 1.3911506872048797E-003 - 52.000000000000000 1.3968634525818040E-003 - 53.000000000000000 1.4027052080686439E-003 - 54.000000000000000 1.4086770557060206E-003 - 55.000000000000000 1.4147801244586476E-003 - 56.000000000000000 1.4210155705813678E-003 - 57.000000000000000 1.4273845779950217E-003 - 58.000000000000000 1.4338883586723259E-003 - 59.000000000000000 1.4405281530339581E-003 - 60.000000000000000 1.4473052303550491E-003 - 61.000000000000000 1.4542208891822866E-003 - 62.000000000000000 1.4612764577618328E-003 - 63.000000000000000 1.4684732944782823E-003 - 64.000000000000000 1.4758127883048875E-003 - 65.000000000000000 1.4832963592652608E-003 - 66.000000000000000 1.4909254589068135E-003 - 67.000000000000000 1.4987015707861599E-003 - 68.000000000000000 1.5066262109667365E-003 - 69.000000000000000 1.5147009285289086E-003 - 70.000000000000000 1.5229273060928051E-003 - 71.000000000000000 1.5313069603541743E-003 - 72.000000000000000 1.5398415426335313E-003 - 73.000000000000000 1.5485327394388858E-003 - 74.000000000000000 1.5573822730423456E-003 - 75.000000000000000 1.5663919020709037E-003 - 76.000000000000000 1.5755634221117198E-003 - 77.000000000000000 1.5848986663322198E-003 - 78.000000000000000 1.5943995061153374E-003 - 79.000000000000000 1.6040678517102553E-003 - 80.000000000000000 1.6139056528989722E-003 - 81.000000000000000 1.6239148996790818E-003 - 82.000000000000000 1.6340976229631109E-003 - 83.000000000000000 1.6444558952948157E-003 - 84.000000000000000 1.6549918315828139E-003 - 85.000000000000000 1.6657075898519639E-003 - 86.000000000000000 1.6766053720129090E-003 - 87.000000000000000 1.6876874246501943E-003 - 88.000000000000000 1.6989560398294184E-003 - 89.000000000000000 1.7104135559238519E-003 - 90.000000000000000 1.7220623584609926E-003 - 91.000000000000000 1.7339048809895405E-003 - 92.000000000000000 1.7459436059672681E-003 - 93.000000000000000 1.7581810656702973E-003 - 94.000000000000000 1.7706198431243151E-003 - 95.000000000000000 1.7832625730582283E-003 - 96.000000000000000 1.7961119428808376E-003 - 97.000000000000000 1.8091706936810869E-003 - 98.000000000000000 1.8224416212524442E-003 - 99.000000000000000 1.8359275771420555E-003 - 100.00000000000000 1.8496314697252390E-003 - 101.00000000000000 1.8635562653059900E-003 - 102.00000000000000 1.8777049892441181E-003 - 103.00000000000000 1.8920807271096969E-003 - 104.00000000000000 1.9066866258655249E-003 - 105.00000000000000 1.9215258950782826E-003 - 106.00000000000000 1.9366018081591314E-003 - 107.00000000000000 1.9519177036344986E-003 - 108.00000000000000 1.9674769864478194E-003 - 109.00000000000000 1.9832831292930296E-003 - 110.00000000000000 1.9993396739806271E-003 - 111.00000000000000 2.0156502328371496E-003 - 112.00000000000000 2.0322184901389144E-003 - 113.00000000000000 2.0490482035809407E-003 - 114.00000000000000 2.0661432057819409E-003 - 115.00000000000000 2.0835074058263544E-003 - 116.00000000000000 2.1011447908443766E-003 - 117.00000000000000 2.1190594276309835E-003 - 118.00000000000000 2.1372554643050061E-003 - 119.00000000000000 2.1557371320092872E-003 - 120.00000000000000 2.1745087466530297E-003 - 121.00000000000000 2.1935747106974638E-003 - 122.00000000000000 2.2129395149859928E-003 - 123.00000000000000 2.2326077406199878E-003 - 124.00000000000000 2.2525840608815201E-003 - 125.00000000000000 2.2728732432042223E-003 - 126.00000000000000 2.2934801511936545E-003 - 127.00000000000000 2.3144097466984742E-003 - 128.00000000000000 2.3356670919338234E-003 - 129.00000000000000 2.3572573516583463E-003 - 130.00000000000000 2.3791857954063218E-003 - 131.00000000000000 2.4014577997764183E-003 - 132.00000000000000 2.4240788507786464E-003 - 133.00000000000000 2.4470545462411133E-003 - 134.00000000000000 2.4703905982782432E-003 - 135.00000000000000 2.4940928358221821E-003 - 136.00000000000000 2.5181672072191520E-003 - 137.00000000000000 2.5426197828925666E-003 - 138.00000000000000 2.5674567580748066E-003 - 139.00000000000000 2.5926844556095809E-003 - 140.00000000000000 2.6183093288268604E-003 - 141.00000000000000 2.6443379644924738E-003 - 142.00000000000000 2.6707770858344676E-003 - 143.00000000000000 2.6976335556484383E-003 - 144.00000000000000 2.7249143794840724E-003 - 145.00000000000000 2.7526267089152795E-003 - 146.00000000000000 2.7807778448962765E-003 - 147.00000000000000 2.8093752412061343E-003 - 148.00000000000000 2.8384265079843445E-003 - 149.00000000000000 2.8679394153600731E-003 - 150.00000000000000 2.8979218971777942E-003 - 151.00000000000000 2.9283820548221678E-003 - 152.00000000000000 2.9593281611450524E-003 - 153.00000000000000 2.9907686644976583E-003 - 154.00000000000000 3.0227121928709582E-003 - 155.00000000000000 3.0551675581475608E-003 - 156.00000000000000 3.0881437604683331E-003 - 157.00000000000000 3.1216499927172472E-003 - 158.00000000000000 3.1556956451279272E-003 - 159.00000000000000 3.1902903100155676E-003 - 160.00000000000000 3.2254437866380244E-003 - 161.00000000000000 3.2611660861899218E-003 - 162.00000000000000 3.2974674369338304E-003 - 163.00000000000000 3.3343582894726923E-003 - 164.00000000000000 3.3718493221677566E-003 - 165.00000000000000 3.4099514467064965E-003 - 166.00000000000000 3.4486758138250920E-003 - 167.00000000000000 3.4880338191902345E-003 - 168.00000000000000 3.5280371094451681E-003 - 169.00000000000000 3.5686975884250227E-003 - 170.00000000000000 3.6100274235467147E-003 - 171.00000000000000 3.6520390523788592E-003 - 172.00000000000000 3.6947451893972437E-003 - 173.00000000000000 3.7381588329317848E-003 - 174.00000000000000 3.7822932723108639E-003 - 175.00000000000000 3.8271620952093398E-003 - 176.00000000000000 3.8727791952066653E-003 - 177.00000000000000 3.9191587795617197E-003 - 178.00000000000000 3.9663153772113192E-003 - 179.00000000000000 4.0142638469994902E-003 - 180.00000000000000 4.0630193861449105E-003 - 181.00000000000000 4.1125975389541799E-003 - 182.00000000000000 4.1730182948308721E-003 - 183.00000000000000 4.2243216838178509E-003 - 184.00000000000000 4.2764966192740906E-003 - 185.00000000000000 4.3295601267864465E-003 - 186.00000000000000 4.3835296265737102E-003 - 187.00000000000000 4.4384229438513287E-003 - 188.00000000000000 4.4942583195036415E-003 - 189.00000000000000 4.5510544210738031E-003 - 190.00000000000000 4.6088303540819170E-003 - 191.00000000000000 4.6676056736822073E-003 - 192.00000000000000 4.7274003966705615E-003 - 193.00000000000000 4.7882350138541037E-003 - 194.00000000000000 4.8501305027949596E-003 - 195.00000000000000 4.9131083409407516E-003 - 196.00000000000000 4.9771905191548999E-003 - 197.00000000000000 5.0423995556602109E-003 - 198.00000000000000 5.1087585104098209E-003 - 199.00000000000000 5.1762909999000378E-003 - 200.00000000000000 5.2450212124401652E-003 - 201.00000000000000 5.3149739238949898E-003 - 202.00000000000000 5.3861745139161931E-003 - 203.00000000000000 5.4586489826795757E-003 - 204.00000000000000 5.5324239681455627E-003 - 205.00000000000000 5.6075267638612796E-003 - 206.00000000000000 5.6839853373229464E-003 - 207.00000000000000 5.7618283489183318E-003 - 208.00000000000000 5.8410851714695454E-003 - 209.00000000000000 5.9217859103973891E-003 - 210.00000000000000 6.0039614245291544E-003 - 211.00000000000000 6.0876433475728145E-003 - 212.00000000000000 6.1728641102812435E-003 - 213.00000000000000 6.2596569633311689E-003 - 214.00000000000000 6.3480560009424600E-003 - 215.00000000000000 6.4380961852644313E-003 - 216.00000000000000 6.5298133715567929E-003 - 217.00000000000000 6.6232443341941475E-003 - 218.00000000000000 6.7184267935238486E-003 - 219.00000000000000 6.8153994436085368E-003 - 220.00000000000000 6.9142019808855985E-003 - 221.00000000000000 7.0148751337773766E-003 - 222.00000000000000 7.1174606932871524E-003 - 223.00000000000000 7.2220015446174245E-003 - 224.00000000000000 7.3285416998484692E-003 - 225.00000000000000 7.4371263317167639E-003 - 226.00000000000000 7.5478018085344111E-003 - 227.00000000000000 7.6606157302924541E-003 - 228.00000000000000 7.7756169659927770E-003 - 229.00000000000000 7.8928556922550196E-003 - 230.00000000000000 8.0123834332470646E-003 - 231.00000000000000 8.1342531019894544E-003 - 232.00000000000000 8.2585190430864604E-003 - 233.00000000000000 8.3852370769385079E-003 - 234.00000000000000 8.5144645454932541E-003 - 235.00000000000000 8.6462603595946833E-003 - 236.00000000000000 8.7806850479925959E-003 - 237.00000000000000 8.9178008080769821E-003 - 238.00000000000000 9.0576715584051136E-003 - 239.00000000000000 9.2003629930916733E-003 - 240.00000000000000 9.3459426381354450E-003 - 241.00000000000000 9.4944799097594658E-003 - 242.00000000000000 9.6460461748445157E-003 - 243.00000000000000 9.8007148135396769E-003 - 244.00000000000000 9.9585612841372653E-003 - 245.00000000000000 1.0119663190303147E-002 - 246.00000000000000 1.0284100350757738E-002 - 247.00000000000000 1.0451954871507015E-002 - 248.00000000000000 1.0623311220727558E-002 - 249.00000000000000 1.0798256306413921E-002 - 250.00000000000000 1.0976879556902018E-002 - 251.00000000000000 1.1159273004386849E-002 - 252.00000000000000 1.1345531371558734E-002 - 253.00000000000000 1.1535752161487469E-002 - 254.00000000000000 1.1730035750890079E-002 - 255.00000000000000 1.1928485486923908E-002 - 256.00000000000000 1.2131207787653352E-002 - 257.00000000000000 1.2338312246345502E-002 - 258.00000000000000 1.2549911739757138E-002 - 259.00000000000000 1.2766122540583029E-002 - 260.00000000000000 1.2987064434243741E-002 - 261.00000000000000 1.3212860840199205E-002 - 262.00000000000000 1.3443638937983457E-002 - 263.00000000000000 1.3679529798165002E-002 - 264.00000000000000 1.3920668518447202E-002 - 265.00000000000000 1.4167194365133235E-002 - 266.00000000000000 1.4419250920191131E-002 - 267.00000000000000 1.4676986234165679E-002 - 268.00000000000000 1.4940552985196020E-002 - 269.00000000000000 1.5210108644410412E-002 - 270.00000000000000 1.5485815647982959E-002 - 271.00000000000000 1.5767841576151079E-002 - 272.00000000000000 1.6056359339507314E-002 - 273.00000000000000 1.6351547372894584E-002 - 274.00000000000000 1.6653589837250633E-002 - 275.00000000000000 1.6962676829764385E-002 - 276.00000000000000 1.7279004602725481E-002 - 277.00000000000000 1.7602775791467590E-002 - 278.00000000000000 1.7934199651826136E-002 - 279.00000000000000 1.8273492307552812E-002 - 280.00000000000000 1.8620877008151869E-002 - 281.00000000000000 1.8976584397627255E-002 - 282.00000000000000 1.9340852794654442E-002 - 283.00000000000000 1.9713928484718531E-002 - 284.00000000000000 2.0096066024787265E-002 - 285.00000000000000 2.0487528561118414E-002 - 286.00000000000000 2.0888588160832215E-002 - 287.00000000000000 2.1299526157912838E-002 - 288.00000000000000 2.1720633514338408E-002 - 289.00000000000000 2.2152211197076520E-002 - 290.00000000000000 2.2594570571721605E-002 - 291.00000000000000 2.3048033813592564E-002 - 292.00000000000000 2.3512934337153655E-002 - 293.00000000000000 2.3989617244668075E-002 - 294.00000000000000 2.4478439795044493E-002 - 295.00000000000000 2.4979771893888478E-002 - 296.00000000000000 2.5493996605827855E-002 - 297.00000000000000 2.6021510690239679E-002 - 298.00000000000000 2.6562725161570114E-002 - 299.00000000000000 2.7118065875504745E-002 - 300.00000000000000 2.7687974142318595E-002 - 301.00000000000000 2.8272907368810076E-002 - 302.00000000000000 2.8873339730303021E-002 - 303.00000000000000 2.9489762874286656E-002 - 304.00000000000000 3.0122686657352959E-002 - 305.00000000000000 3.0772639917187820E-002 - 306.00000000000000 3.1440171281473826E-002 - 307.00000000000000 3.2125850015671850E-002 - 308.00000000000000 3.2830266911764008E-002 - 309.00000000000000 3.3554035220163626E-002 - 310.00000000000000 3.4297791627129350E-002 - 311.00000000000000 3.5062197280159559E-002 - 312.00000000000000 3.5847938863992790E-002 - 313.00000000000000 3.6655729729997970E-002 - 314.00000000000000 3.7486311081907003E-002 - 315.00000000000000 3.8340453221023932E-002 - 316.00000000000000 3.9218956854235126E-002 - 317.00000000000000 4.0122654468352767E-002 - 318.00000000000000 4.1052411774540720E-002 - 319.00000000000000 4.2009129226807745E-002 - 320.00000000000000 4.2993743618802205E-002 - 321.00000000000000 4.4007229763409865E-002 - 322.00000000000000 4.5050602259942674E-002 - 323.00000000000000 4.6124917354011044E-002 - 324.00000000000000 4.7231274895500608E-002 - 325.00000000000000 4.8370820400423085E-002 - 326.00000000000000 4.9544747222786539E-002 - 327.00000000000000 5.0754298843032057E-002 - 328.00000000000000 5.2000771280012958E-002 - 329.00000000000000 5.3285515633955075E-002 - 330.00000000000000 5.4609940768330839E-002 - 331.00000000000000 5.5975516139110490E-002 - 332.00000000000000 5.7383774780423016E-002 - 333.00000000000000 5.8836316456272353E-002 - 334.00000000000000 6.0334810988609065E-002 - 335.00000000000000 6.1881001772766543E-002 - 336.00000000000000 6.3476709492027106E-002 - 337.00000000000000 6.5123836043901817E-002 - 338.00000000000000 6.6824368691585215E-002 - 339.00000000000000 6.8580384454992482E-002 - 340.00000000000000 7.0394054756804908E-002 - 341.00000000000000 7.2267650340046452E-002 - 342.00000000000000 7.4203546474899101E-002 - 343.00000000000000 7.6204228473738117E-002 - 344.00000000000000 7.8272297534747376E-002 - 345.00000000000000 8.0410476935960026E-002 - 346.00000000000000 8.2621618603174990E-002 - 347.00000000000000 8.4908710076935068E-002 - 348.00000000000000 8.7274881905626517E-002 - 349.00000000000000 8.9723415493789913E-002 - 350.00000000000000 9.2257751436925217E-002 - 351.00000000000000 9.4881498376452864E-002 - 352.00000000000000 9.7598442411065217E-002 - 353.00000000000000 0.10041255710349541 - 354.00000000000000 0.10332801412475574 - 355.00000000000000 0.10634919458118347 - 356.00000000000000 0.10948070107319365 - 357.00000000000000 0.11272737053851653 - 358.00000000000000 0.11609428793690132 - 359.00000000000000 0.11956681095755470 - 360.00000000000000 0.12319058655672448 - 361.00000000000000 0.20215606231105379 - 362.00000000000000 0.20767836655079511 - 363.00000000000000 0.21339757104679685 - 364.00000000000000 0.21932309981141065 - 365.00000000000000 0.22547112302603478 - 366.00000000000000 0.23184378314292695 - 367.00000000000000 0.23845117169658930 - 368.00000000000000 0.24531076919203776 - 369.00000000000000 0.25243074132053034 - 370.00000000000000 0.25981965000064650 - 371.00000000000000 0.26749299089800815 - 372.00000000000000 0.27547045567711703 - 373.00000000000000 0.28375710030617890 - 374.00000000000000 0.29236763916203834 - 375.00000000000000 0.30132993283334525 - 376.00000000000000 0.31065146905631030 - 377.00000000000000 0.32034714134579778 - 378.00000000000000 0.33044261174769046 - 379.00000000000000 0.34096159348153232 - 380.00000000000000 0.35191439794373885 - 381.00000000000000 0.36332296755228577 - 382.00000000000000 0.37523253789844113 - 383.00000000000000 0.38764758539697330 - 384.00000000000000 0.40059394198574017 - 385.00000000000000 0.41411473244920344 - 386.00000000000000 0.42823772087934930 - 387.00000000000000 0.44298275905051820 - 388.00000000000000 0.45838796882698890 - 389.00000000000000 0.47451481465346507 - 390.00000000000000 0.49137327862230068 - 391.00000000000000 0.50900398863576979 - 392.00000000000000 0.52748101490987698 - 393.00000000000000 0.54683467109905082 - 394.00000000000000 0.56710244968828383 - 395.00000000000000 0.58835288386395801 - 396.00000000000000 0.61066753969252008 - 397.00000000000000 0.63406993565488778 - 398.00000000000000 0.65862489001281666 - 399.00000000000000 0.68446344414786442 - 400.00000000000000 0.71161260859447073 - 401.00000000000000 0.74014341279599116 - 402.00000000000000 0.77018232901078987 - 403.00000000000000 0.80183662937043043 - 404.00000000000000 0.83515728202995065 - 405.00000000000000 0.87025409194107461 - 406.00000000000000 0.90736005332093461 - 407.00000000000000 0.94648891155794024 - 408.00000000000000 0.98777520868187607 - 409.00000000000000 1.0314615474044391 - 410.00000000000000 1.0776831936316438 - 411.00000000000000 1.1265483490871420 - 412.00000000000000 1.1781379027335419 - 413.00000000000000 1.2322390642582450 - 414.00000000000000 1.2887605860153917 - 415.00000000000000 1.3478388596435678 - 416.00000000000000 1.4098492302337373 - 417.00000000000000 1.4748123692744954 - 418.00000000000000 1.5428240767018095 - 419.00000000000000 1.6141989538573758 - 420.00000000000000 1.6892819247025452 - 421.00000000000000 1.7680162980503829 - 422.00000000000000 1.8506228702063829 - 423.00000000000000 1.9378440080421810 - 424.00000000000000 2.0295018786282628 - 425.00000000000000 2.1258181935060092 - 426.00000000000000 2.2274413067773802 - 427.00000000000000 2.3347138818077733 - 428.00000000000000 2.4475928828467048 - 429.00000000000000 2.5664920279678136 - 430.00000000000000 2.6926181354813652 - 431.00000000000000 2.8254607617069070 - 432.00000000000000 2.9654303946603466 - 433.00000000000000 3.1137265704532910 - 434.00000000000000 3.2703892004237818 - 435.00000000000000 3.4354237516213146 - 436.00000000000000 3.6095859988127006 - 437.00000000000000 3.7939003021594333 - 438.00000000000000 3.9876027391801281 - 439.00000000000000 4.1910534762001381 - 440.00000000000000 4.4045901165285013 - 441.00000000000000 4.3237452455533312 - 442.00000000000000 3.7472280632877575 - 443.00000000000000 3.0641344129130803 - 444.00000000000000 2.3373994510541216 - 445.00000000000000 1.6235270526813124 - 446.00000000000000 0.92090939737629762 - 447.00000000000000 0.22150738307397194 - 448.00000000000000 -0.45377732487761319 - 449.00000000000000 -1.1161994693693704 - 450.00000000000000 -1.7657590504077436 - 451.00000000000000 -2.4024560679991263 - 452.00000000000000 -3.0262905221498708 - 453.00000000000000 -3.6372624128662738 - 454.00000000000000 -4.2353717401545845 - 455.00000000000000 -4.8206185040210014 - 456.00000000000000 -5.3930027044716811 - 457.00000000000000 -5.9525243415127242 - 458.00000000000000 -6.4991834151501902 - 459.00000000000000 -7.0329799253900882 - 460.00000000000000 -7.5539138722383834 - 461.00000000000000 -8.0619852557009928 - 462.00000000000000 -8.5571940757837890 - 463.00000000000000 -9.0395403324925923 - 464.00000000000000 -9.5090240258331917 - 465.00000000000000 -9.9656451558113179 - 466.00000000000000 -10.409403722432666 - 467.00000000000000 -10.840299725702879 - 468.00000000000000 -11.258333165627567 - 469.00000000000000 -11.663504042212283 - 470.00000000000000 -12.055812355462551 - 471.00000000000000 -12.435258105383843 - 472.00000000000000 -12.801841291981589 - 473.00000000000000 -13.155561915261183 - 474.00000000000000 -13.496419975227973 - 475.00000000000000 -13.824415471887267 - 476.00000000000000 -14.139548405244325 - 477.00000000000000 -14.441818775304380 - 478.00000000000000 -14.731226582072612 - 479.00000000000000 -15.007771825554169 - 480.00000000000000 -15.271454505754155 - 481.00000000000000 -15.522274622677635 - 482.00000000000000 -15.760232176329636 - 483.00000000000000 -15.985327166715141 - 484.00000000000000 -16.197559593839109 - 485.00000000000000 -16.396929457706445 - 486.00000000000000 -16.583436758322016 - 487.00000000000000 -16.757081495690663 - 488.00000000000000 -16.917863669817187 - 489.00000000000000 -17.065783280706341 - 490.00000000000000 -17.200840328362855 - 491.00000000000000 -17.323034812791413 - 492.00000000000000 -17.432366733996663 - 493.00000000000000 -17.528836091983223 - 494.00000000000000 -17.612442886755673 - 495.00000000000000 -17.683187118318560 - 496.00000000000000 -17.741068786676387 - 497.00000000000000 -17.786087891833624 - 498.00000000000000 -17.818244433794721 - 499.00000000000000 -17.837538412564072 - 500.00000000000000 -17.843969828146058 - 501.00000000000000 -17.837538680545002 - 502.00000000000000 -17.818244969765221 - 503.00000000000000 -17.786088695810978 - 504.00000000000000 -17.741069858686508 - 505.00000000000000 -17.683188458396021 - 506.00000000000000 -17.612444494943681 - 507.00000000000000 -17.528837968333633 - 508.00000000000000 -17.432368878569978 - 509.00000000000000 -17.323037225656801 - 510.00000000000000 -17.200843009598131 - 511.00000000000000 -17.065786230397993 - 512.00000000000000 -16.917866888060356 - 513.00000000000000 -16.757084982589177 - 514.00000000000000 -16.583440513988378 - 515.00000000000000 -16.396933482261844 - 516.00000000000000 -16.197563887413430 - 517.00000000000000 -15.985331729446967 - 518.00000000000000 -15.760237008366257 - 519.00000000000000 -15.522279724175064 - 520.00000000000000 -15.271459876877133 - 521.00000000000000 -15.007777466476169 - 522.00000000000000 -14.731232492975856 - 523.00000000000000 -14.441824956379850 - 524.00000000000000 -14.139554856691770 - 525.00000000000000 -13.824422193915218 - 526.00000000000000 -13.496426968053756 - 527.00000000000000 -13.155569179110930 - 528.00000000000000 -12.801848827090250 - 529.00000000000000 -12.435265911995202 - 530.00000000000000 -12.055820433829242 - 531.00000000000000 -11.663512392595802 - 532.00000000000000 -11.258341788298289 - 533.00000000000000 -10.840308620940073 - 534.00000000000000 -10.409412890524514 - 535.00000000000000 -9.9656545970549288 - 536.00000000000000 -9.5090337405346190 - 537.00000000000000 -9.0395503209668586 - 538.00000000000000 -8.5572043383548966 - 539.00000000000000 -8.0619957927019481 - 540.00000000000000 -7.5539246840112124 - 541.00000000000000 -7.0329910122858621 - 542.00000000000000 -6.4991947775290413 - 543.00000000000000 -5.9525359797438711 - 544.00000000000000 -5.3930146189334502 - 545.00000000000000 -4.8206306951008475 - 546.00000000000000 -4.2353842082491147 - 547.00000000000000 -3.6372751583812724 - 548.00000000000000 -3.0263035455003222 - 549.00000000000000 -2.4024693696092405 - 550.00000000000000 -1.7657726307109816 - 551.00000000000000 -1.1162133288084697 - 552.00000000000000 -0.45379146390461472 - 553.00000000000000 0.22149296399769688 - 554.00000000000000 0.92089469760988818 - 555.00000000000000 1.6235120718872194 - 556.00000000000000 2.3373841887665656 - 557.00000000000000 3.0641188686330736 - 558.00000000000000 3.7472122989084702 - 559.00000000000000 4.3237293854665744 - 560.00000000000000 4.4045747218253233 - 561.00000000000000 4.1910388896988220 - 562.00000000000000 3.9875889581923225 - 563.00000000000000 3.7938873221348377 - 564.00000000000000 3.6095738133337378 - 565.00000000000000 3.4354123523976066 - 566.00000000000000 3.2703785772863987 - 567.00000000000000 3.1137167113487938 - 568.00000000000000 2.9654212856446143 - 569.00000000000000 2.8254523869387467 - 570.00000000000000 2.6926104772159052 - 571.00000000000000 2.5664850665496783 - 572.00000000000000 2.4475865967029722 - 573.00000000000000 2.3347082474407097 - 574.00000000000000 2.2274362987569556 - 575.00000000000000 2.1258137844621876 - 576.00000000000000 2.0294980392430695 - 577.00000000000000 1.9378407070415469 - 578.00000000000000 1.8506200743517576 - 579.00000000000000 1.7680139721301240 - 580.00000000000000 1.6892800315230763 - 581.00000000000000 1.6141974542341133 - 582.00000000000000 1.5428229294498581 - 583.00000000000000 1.4748115311991106 - 584.00000000000000 1.4098486561205510 - 585.00000000000000 1.3478385022485782 - 586.00000000000000 1.2887603960547158 - 587.00000000000000 1.2322389903976294 - 588.00000000000000 1.1781378915776335 - 589.00000000000000 1.1265483490871420 - 590.00000000000000 1.0776831936316438 - 591.00000000000000 1.0314615474044393 - 592.00000000000000 0.98777520868187607 - 593.00000000000000 0.94648891155794013 - 594.00000000000000 0.90736005332093472 - 595.00000000000000 0.87025409194107450 - 596.00000000000000 0.83515728202995076 - 597.00000000000000 0.80183662937043032 - 598.00000000000000 0.77018232901078976 - 599.00000000000000 0.74014341279599116 - 600.00000000000000 0.71161260859447073 - 601.00000000000000 0.68446344414786442 - 602.00000000000000 0.65862489001281654 - 603.00000000000000 0.63406993565488778 - 604.00000000000000 0.61066753969252008 - 605.00000000000000 0.58835288386395801 - 606.00000000000000 0.56710244968828372 - 607.00000000000000 0.54683467109905082 - 608.00000000000000 0.52748101490987698 - 609.00000000000000 0.50900398863576968 - 610.00000000000000 0.49137327862230062 - 611.00000000000000 0.47451481465346507 - 612.00000000000000 0.45838796882698890 - 613.00000000000000 0.44298275905051826 - 614.00000000000000 0.42823772087934930 - 615.00000000000000 0.41411473244920344 - 616.00000000000000 0.40059394198574022 - 617.00000000000000 0.38764758539697330 - 618.00000000000000 0.37523253789844113 - 619.00000000000000 0.36332296755228582 - 620.00000000000000 0.35191439794373891 - 621.00000000000000 0.34096159348153232 - 622.00000000000000 0.33044261174769046 - 623.00000000000000 0.32034714134579784 - 624.00000000000000 0.31065146905631036 - 625.00000000000000 0.30132993283334530 - 626.00000000000000 0.29236763916203834 - 627.00000000000000 0.28375710030617884 - 628.00000000000000 0.27547045567711698 - 629.00000000000000 0.26749299089800815 - 630.00000000000000 0.25981965000064655 - 631.00000000000000 0.25243074132053034 - 632.00000000000000 0.24531076919203776 - 633.00000000000000 0.23845117169658933 - 634.00000000000000 0.23184378314292695 - 635.00000000000000 0.22547112302603478 - 636.00000000000000 0.21932309981141065 - 637.00000000000000 0.21339757104679688 - 638.00000000000000 0.20767836655079513 - 639.00000000000000 0.20215606231105382 - 640.00000000000000 0.19682582761625819 - 641.00000000000000 0.19167982330775035 - 642.00000000000000 0.18672696654771398 - 643.00000000000000 0.18192054139684583 - 644.00000000000000 0.17727890418311293 - 645.00000000000000 0.17278991703524316 - 646.00000000000000 0.16844748741350876 - 647.00000000000000 0.16424900989967869 - 648.00000000000000 0.16018736546821363 - 649.00000000000000 0.15625542794898403 - 650.00000000000000 0.15244910084499205 - 651.00000000000000 0.14876591644195133 - 652.00000000000000 0.14519782951023819 - 653.00000000000000 0.14174046701721671 - 654.00000000000000 0.13839255824457214 - 655.00000000000000 0.13514781451803515 - 656.00000000000000 0.13200165216490886 - 657.00000000000000 0.12895160852273427 - 658.00000000000000 0.12599492383504082 - 659.00000000000000 0.12312622039848933 - 660.00000000000000 0.12034230949484220 - 661.00000000000000 0.11764271949271111 - 662.00000000000000 0.11502197461611340 - 663.00000000000000 0.11247713539610175 - 664.00000000000000 0.11000677535126910 - 665.00000000000000 0.10760812700883783 - 666.00000000000000 0.10527757823651729 - 667.00000000000000 0.10301298705578783 - 668.00000000000000 0.10081370407939616 - 669.00000000000000 9.8675646695142380E-002 - 670.00000000000000 9.6596733867987283E-002 - 671.00000000000000 9.4576199022144605E-002 - 672.00000000000000 9.2611410021550877E-002 - 673.00000000000000 9.0699938520709733E-002 - 674.00000000000000 8.8840429484792899E-002 - 675.00000000000000 8.7031917266284592E-002 - 676.00000000000000 8.5271563056124128E-002 - 677.00000000000000 8.3557798664696439E-002 - 678.00000000000000 8.1890277424310459E-002 - 679.00000000000000 8.0266567842262979E-002 - 680.00000000000000 7.8685043792767945E-002 - 681.00000000000000 7.7144874253372214E-002 - 682.00000000000000 7.5644950036665298E-002 - 683.00000000000000 7.4183291895265444E-002 - 684.00000000000000 7.2758704195004742E-002 - 685.00000000000000 7.1371105514012864E-002 - 686.00000000000000 7.0018287782061567E-002 - 687.00000000000000 6.8699173155846380E-002 - 688.00000000000000 6.7413280674158596E-002 - 689.00000000000000 6.6159459959854769E-002 - 690.00000000000000 6.4936331692535265E-002 - 691.00000000000000 6.3743100248271650E-002 - 692.00000000000000 6.2579471415829382E-002 - 693.00000000000000 6.1443838156918514E-002 - 694.00000000000000 6.0335373868345567E-002 - 695.00000000000000 5.9253831366600347E-002 - 696.00000000000000 5.8198078528077887E-002 - 697.00000000000000 5.7167158217017247E-002 - 698.00000000000000 5.6160558893561835E-002 - 699.00000000000000 5.5177845798565002E-002 - 700.00000000000000 5.4217863533364788E-002 - 701.00000000000000 5.3279965979651928E-002 - 702.00000000000000 5.2364061336419972E-002 - 703.00000000000000 5.1469068883887095E-002 - 704.00000000000000 5.0594330570570840E-002 - 705.00000000000000 4.9739532188371835E-002 - 706.00000000000000 4.8904164548527411E-002 - 707.00000000000000 4.8087397979924298E-002 - 708.00000000000000 4.7288739231721132E-002 - 709.00000000000000 4.6508149579767247E-002 - 710.00000000000000 4.5744665849155201E-002 - 711.00000000000000 4.4997826565984569E-002 - 712.00000000000000 4.4267454594403767E-002 - 713.00000000000000 4.3553008838545411E-002 - 714.00000000000000 4.2853896091344032E-002 - 715.00000000000000 4.2169786340945391E-002 - 716.00000000000000 4.1500529857782770E-002 - 717.00000000000000 4.0845416044326281E-002 - 718.00000000000000 4.0204076575993876E-002 - 719.00000000000000 3.9576429233948190E-002 - 720.00000000000000 3.8961928056838478E-002 - 721.00000000000000 3.8360152418914592E-002 - 722.00000000000000 3.7770889438067598E-002 - 723.00000000000000 3.7193920290172600E-002 - 724.00000000000000 3.6628721584536783E-002 - 725.00000000000000 3.6074997087473720E-002 - 726.00000000000000 3.5532729838403804E-002 - 727.00000000000000 3.5001385848684771E-002 - 728.00000000000000 3.4480670966863129E-002 - 729.00000000000000 3.3970456374860895E-002 - 730.00000000000000 3.3470482814120253E-002 - 731.00000000000000 3.2980366390782077E-002 - 732.00000000000000 3.2499885615723638E-002 - 733.00000000000000 3.2029012398702816E-002 - 734.00000000000000 3.1567290256334243E-002 - 735.00000000000000 3.1114499504827434E-002 - 736.00000000000000 3.0670571534666773E-002 - 737.00000000000000 3.0235226193148244E-002 - 738.00000000000000 2.9808183988412906E-002 - 739.00000000000000 2.9389296068666325E-002 - 740.00000000000000 2.8978478688023938E-002 - 741.00000000000000 2.8575388204418833E-002 - 742.00000000000000 2.8179844963048331E-002 - 743.00000000000000 2.7791823266279037E-002 - 744.00000000000000 2.7411035296325978E-002 - 745.00000000000000 2.7037279964544809E-002 - 746.00000000000000 2.6670462006151918E-002 - 747.00000000000000 2.6310461102025968E-002 - 748.00000000000000 2.5957019581810542E-002 - 749.00000000000000 2.5609989616343247E-002 - 750.00000000000000 2.5269375829575676E-002 - 751.00000000000000 2.4934891582875485E-002 - 752.00000000000000 2.4606394807498121E-002 - 753.00000000000000 2.4283829718863833E-002 - 754.00000000000000 2.3967053046550134E-002 - 755.00000000000000 2.3655872860561922E-002 - 756.00000000000000 2.3350182231881594E-002 - 757.00000000000000 2.3049960648894439E-002 - 758.00000000000000 2.2754974710119519E-002 - 759.00000000000000 2.2465111410781242E-002 - 760.00000000000000 2.2180343954422831E-002 - 761.00000000000000 2.1900515822827798E-002 - 762.00000000000000 2.1625485501298048E-002 - 763.00000000000000 2.1355181376686905E-002 - 764.00000000000000 2.1089552989145281E-002 - 765.00000000000000 2.0828421932410007E-002 - 766.00000000000000 2.0571693948829572E-002 - 767.00000000000000 2.0319363398577487E-002 - 768.00000000000000 2.0071267394891165E-002 - 769.00000000000000 1.9827303166308288E-002 - 770.00000000000000 1.9587425750335268E-002 - 771.00000000000000 1.9351564164693955E-002 - 772.00000000000000 1.9119582845056663E-002 - 773.00000000000000 1.8891403180723781E-002 - 774.00000000000000 1.8667033597476981E-002 - 775.00000000000000 1.8446310906346790E-002 - 776.00000000000000 1.8229161735692140E-002 - 777.00000000000000 1.8015561186808252E-002 - 778.00000000000000 1.7805424714368349E-002 - 779.00000000000000 1.7598650211370175E-002 - 780.00000000000000 1.7395183022348493E-002 - 781.00000000000000 1.7195007972322103E-002 - 782.00000000000000 1.6997998197970241E-002 - 783.00000000000000 1.6804091769382998E-002 - 784.00000000000000 1.6613278828341193E-002 - 785.00000000000000 1.6425466396661478E-002 - 786.00000000000000 1.6240578543692329E-002 - 787.00000000000000 1.6058579067756845E-002 - 788.00000000000000 1.5879435742429140E-002 - 789.00000000000000 1.5703050465128657E-002 - 790.00000000000000 1.5529370730227140E-002 - 791.00000000000000 1.5358397879165835E-002 - 792.00000000000000 1.5190034269608002E-002 - 793.00000000000000 1.5024224474121998E-002 - 794.00000000000000 1.4860946483311042E-002 - 795.00000000000000 1.4700155988965073E-002 - 796.00000000000000 1.4541777709063007E-002 - 797.00000000000000 1.4385769352419501E-002 - 798.00000000000000 1.4232133279758852E-002 - 799.00000000000000 1.4080777175136588E-002 - 800.00000000000000 1.3931659264211801E-002 - 801.00000000000000 1.3784768419054131E-002 - 802.00000000000000 1.3640052020327663E-002 - 803.00000000000000 1.3497452922874676E-002 - 804.00000000000000 1.3356941903011775E-002 - 805.00000000000000 1.3218507588744160E-002 - 806.00000000000000 1.3082077478306328E-002 - 807.00000000000000 1.2947615802451686E-002 - 808.00000000000000 1.2815119736236631E-002 - 809.00000000000000 1.2684531206552191E-002 - 810.00000000000000 1.2555807498173440E-002 - 811.00000000000000 1.2428929623850334E-002 - 812.00000000000000 1.2303876041117205E-002 - 813.00000000000000 1.2180590192056199E-002 - 814.00000000000000 1.2059041335108758E-002 - 815.00000000000000 1.1939232936385482E-002 - 816.00000000000000 1.1821103628580416E-002 - 817.00000000000000 1.1704622182492984E-002 - 818.00000000000000 1.1589777611929437E-002 - 819.00000000000000 1.1466568696423448E-002 - 820.00000000000000 1.1354927962697966E-002 - 821.00000000000000 1.1244827629165965E-002 - 822.00000000000000 1.1136267419006149E-002 - 823.00000000000000 1.1029192759102888E-002 - 824.00000000000000 1.0923578740838847E-002 - 825.00000000000000 1.0819420610610727E-002 - 826.00000000000000 1.0716684160802958E-002 - 827.00000000000000 1.0615336047055288E-002 - 828.00000000000000 1.0515360090319071E-002 - 829.00000000000000 1.0416747671761928E-002 - 830.00000000000000 1.0319455596454462E-002 - 831.00000000000000 1.0223462279344854E-002 - 832.00000000000000 1.0128767796266234E-002 - 833.00000000000000 1.0035334342908380E-002 - 834.00000000000000 9.9431369401379172E-003 - 835.00000000000000 9.8521653803766446E-003 - 836.00000000000000 9.7624047371550511E-003 - 837.00000000000000 9.6738210986040789E-003 - 838.00000000000000 9.5863957017491268E-003 - 839.00000000000000 9.5001323383326463E-003 - 840.00000000000000 9.4149908942391496E-003 - 841.00000000000000 9.3309531418945338E-003 - 842.00000000000000 9.2480136075845411E-003 - 843.00000000000000 9.1661526665213454E-003 - 844.00000000000000 9.0853439879362823E-003 - 845.00000000000000 9.0055738683715945E-003 - 846.00000000000000 8.9268409680820643E-003 - 847.00000000000000 8.8491118180304817E-003 - 848.00000000000000 8.7723709823741593E-003 - 849.00000000000000 8.6966167262425478E-003 - 850.00000000000000 8.6218259994250666E-003 - 851.00000000000000 8.5479786659088206E-003 - 852.00000000000000 8.4750655406575694E-003 - 853.00000000000000 8.4030800293105300E-003 - 854.00000000000000 8.3319954891620095E-003 - 855.00000000000000 8.2617984107957974E-003 - 856.00000000000000 8.1924899971071666E-003 - 857.00000000000000 8.1240447784226479E-003 - 858.00000000000000 8.0564476814267566E-003 - 859.00000000000000 2.6443379644924738E-003 - 860.00000000000000 2.6183093288268604E-003 - 861.00000000000000 2.5926844556095809E-003 - 862.00000000000000 2.5674567580748066E-003 - 863.00000000000000 2.5426197828925666E-003 - 864.00000000000000 2.5181672072191520E-003 - 865.00000000000000 2.4940928358221817E-003 - 866.00000000000000 2.4703905982782432E-003 - 867.00000000000000 2.4470545462411133E-003 - 868.00000000000000 2.4240788507786464E-003 - 869.00000000000000 2.4014577997764183E-003 - 870.00000000000000 2.3791857954063222E-003 - 871.00000000000000 2.3572573516583463E-003 - 872.00000000000000 2.3356670919338234E-003 - 873.00000000000000 2.3144097466984742E-003 - 874.00000000000000 2.2934801511936545E-003 - 875.00000000000000 2.2728732432042219E-003 - 876.00000000000000 2.2525840608815196E-003 - 877.00000000000000 2.2326077406199878E-003 - 878.00000000000000 2.2129395149859928E-003 - 879.00000000000000 2.1935747106974642E-003 - 880.00000000000000 2.1745087466530297E-003 - 881.00000000000000 2.1557371320092872E-003 - 882.00000000000000 2.1372554643050061E-003 - 883.00000000000000 2.1190594276309835E-003 - 884.00000000000000 2.1011447908443766E-003 - 885.00000000000000 2.0835074058263544E-003 - 886.00000000000000 2.0661432057819409E-003 - 887.00000000000000 2.0490482035809407E-003 - 888.00000000000000 2.0322184901389148E-003 - 889.00000000000000 2.0156502328371492E-003 - 890.00000000000000 1.9993396739806271E-003 - 891.00000000000000 1.9832831292930296E-003 - 892.00000000000000 1.9674769864478198E-003 - 893.00000000000000 1.9519177036344986E-003 - 894.00000000000000 1.9366018081591314E-003 - 895.00000000000000 1.9215258950782824E-003 - 896.00000000000000 1.9066866258655247E-003 - 897.00000000000000 1.8920807271096969E-003 - 898.00000000000000 1.8777049892441181E-003 - 899.00000000000000 1.8635562653059900E-003 - 900.00000000000000 1.8496314697252390E-003 - 901.00000000000000 1.8359275771420555E-003 - 902.00000000000000 1.8224416212524442E-003 - 903.00000000000000 1.8091706936810869E-003 - 904.00000000000000 1.7961119428808378E-003 - 905.00000000000000 1.7832625730582283E-003 - 906.00000000000000 1.7706198431243151E-003 - 907.00000000000000 1.7581810656702973E-003 - 908.00000000000000 1.7459436059672679E-003 - 909.00000000000000 1.7339048809895405E-003 - 910.00000000000000 1.7220623584609926E-003 - 911.00000000000000 1.7104135559238519E-003 - 912.00000000000000 1.6989560398294184E-003 - 913.00000000000000 1.6876874246501943E-003 - 914.00000000000000 1.6766053720129090E-003 - 915.00000000000000 1.6657075898519639E-003 - 916.00000000000000 1.6549918315828139E-003 - 917.00000000000000 1.6444558952948157E-003 - 918.00000000000000 1.6340976229631112E-003 - 919.00000000000000 1.6239148996790818E-003 - 920.00000000000000 1.6139056528989722E-003 - 921.00000000000000 1.6040678517102553E-003 - 922.00000000000000 1.5943995061153376E-003 - 923.00000000000000 1.5848986663322198E-003 - 924.00000000000000 1.5755634221117198E-003 - 925.00000000000000 1.5663919020709040E-003 - 926.00000000000000 1.5573822730423456E-003 - 927.00000000000000 1.5485327394388858E-003 - 928.00000000000000 1.5398415426335313E-003 - 929.00000000000000 1.5313069603541743E-003 - 930.00000000000000 1.5229273060928053E-003 - 931.00000000000000 1.5147009285289086E-003 - 932.00000000000000 1.5066262109667365E-003 - 933.00000000000000 1.4987015707861599E-003 - 934.00000000000000 1.4909254589068135E-003 - 935.00000000000000 1.4832963592652608E-003 - 936.00000000000000 1.4758127883048875E-003 - 937.00000000000000 1.4684732944782826E-003 - 938.00000000000000 1.4612764577618328E-003 - 939.00000000000000 1.4542208891822866E-003 - 940.00000000000000 1.4473052303550491E-003 - 941.00000000000000 1.4405281530339581E-003 - 942.00000000000000 1.4338883586723259E-003 - 943.00000000000000 1.4273845779950217E-003 - 944.00000000000000 1.4210155705813678E-003 - 945.00000000000000 1.4147801244586476E-003 - 946.00000000000000 1.4086770557060206E-003 - 947.00000000000000 1.4027052080686437E-003 - 948.00000000000000 1.3968634525818040E-003 - 949.00000000000000 1.3911506872048797E-003 - 950.00000000000000 1.3855658364649425E-003 - 951.00000000000000 1.3801078511098319E-003 - 952.00000000000000 1.3747757077705212E-003 - 953.00000000000000 1.3695684086326135E-003 - 954.00000000000000 1.3644849811168085E-003 - 955.00000000000000 1.3595244775681770E-003 - 956.00000000000000 1.3546859749540999E-003 - 957.00000000000000 1.3499685745707120E-003 - 958.00000000000000 1.3453714017577245E-003 - 959.00000000000000 1.3408936056214739E-003 - 960.00000000000000 1.3365343587660638E-003 - 961.00000000000000 1.3322928570324861E-003 - 962.00000000000000 1.3281683192455727E-003 - 963.00000000000000 1.3241599869686746E-003 - 964.00000000000000 1.3202671242659397E-003 - 965.00000000000000 1.3164890174720873E-003 - 966.00000000000000 1.3128249749695513E-003 - 967.00000000000000 1.3092743269729082E-003 - 968.00000000000000 1.3058364253204665E-003 - 969.00000000000000 1.3025106432729382E-003 - 970.00000000000000 1.2992963753190782E-003 - 971.00000000000000 1.2961930369882163E-003 - 972.00000000000000 1.2932000646695839E-003 - 973.00000000000000 1.2903169154383515E-003 - 974.00000000000000 1.2875430668882997E-003 - 975.00000000000000 1.2848780169710419E-003 - 976.00000000000000 1.2823212838417191E-003 - 977.00000000000000 1.2798724057111101E-003 - 978.00000000000000 1.2775309407040661E-003 - 979.00000000000000 1.2752964667242241E-003 - 980.00000000000000 1.2731685813249238E-003 - 981.00000000000000 1.2711469015862706E-003 - 982.00000000000000 1.2692310639982930E-003 - 983.00000000000000 1.2674207243501341E-003 - 984.00000000000000 1.2657155576252257E-003 - 985.00000000000000 1.2641152579024033E-003 - 986.00000000000000 1.2626195382629108E-003 - 987.00000000000000 1.2612281307032500E-003 - 988.00000000000000 1.2599407860538437E-003 - 989.00000000000000 1.2587572739034647E-003 - 990.00000000000000 1.2576773825294040E-003 - 991.00000000000000 1.2567009188333385E-003 - 992.00000000000000 1.2558277082828800E-003 - 993.00000000000000 1.2550575948587633E-003 - 994.00000000000000 1.2543904410076655E-003 - 995.00000000000000 1.2538261276006199E-003 - 996.00000000000000 1.2533645538970216E-003 - 997.00000000000000 1.2530056375141877E-003 - 998.00000000000000 1.2527493144024800E-003 - 999.00000000000000 1.2525955388259645E-003 - 1000.0000000000000 1.2525442833486017E-003 diff --git a/production_compile.bash b/production_compile.bash new file mode 100644 index 00000000..8f07ea27 --- /dev/null +++ b/production_compile.bash @@ -0,0 +1,3 @@ +cmake -P distclean.cmake +cmake -B build -S . +cmake --build build diff --git a/pyproject.toml b/pyproject.toml index ba12d7c1..d1d84406 100644 --- a/pyproject.toml +++ b/pyproject.toml @@ -1,6 +1,6 @@ [project] name = "ctem" -version = "2023.10.0" +version = "2024.2.0" authors=[ {name = 'David A. Minton', email='daminton@purdue.edu'}, {name = 'James E. Richardson'}, diff --git a/src/CMakeLists.txt b/src/CMakeLists.txt index 99ff6ff0..bc1aebc7 100644 --- a/src/CMakeLists.txt +++ b/src/CMakeLists.txt @@ -239,9 +239,9 @@ ELSE() SET(CMAKE_INSTALL_PREFIX /usr/local) ENDIF(WIN32) INSTALL(TARGETS ${CTEM_DRIVER} ${CTEM_LIBRARY} - RUNTIME DESTINATION bin - LIBRARY DESTINATION lib - ARCHIVE DESTINATION lib - INCLUDES DESTINATION include + RUNTIME DESTINATION ${INSTALL_BINDIR} + LIBRARY DESTINATION ${INSTALL_LIBDIR} + ARCHIVE DESTINATION ${INSTALL_LIBDIR} + INCLUDES DESTINATION ${INSTALL_INCLUDEDIR} ) diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 2c42e765..eb858938 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -145,20 +145,26 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! Reset age clock = 0.0_DP finterval = 1.0_DP / real(ntotcrat,kind=DP) - maxage = user%interval * user%numintervals - if (maxage < 0._DP ) then - write(*,*) "MAJOR ERROR: Negative age!" - stop - else if (maxage < 2330._DP) then - maxageGa = util_t_from_scale(maxage,1e-11_DP,4.5_DP) - else - maxageGa = 4.5_DP !util_t_from_scale only supports ages <4.5 Ga + if (user%doregotrack) then + if (user%runtype .eq. 'STATISTICAL') then + maxage = user%interval + else + maxage = user%interval * user%numintervals + end if + if (maxage < 0._DP ) then + write(*,*) "MAJOR ERROR: Negative age!" + stop + else if (maxage < 2330._DP) then + maxageGa = util_t_from_scale(maxage,1e-11_DP,4.5_DP) + else + maxageGa = 4.5_DP !util_t_from_scale only supports ages <4.5 Ga + end if + age_resolution = maxageGa / real(MAXAGEBINS) + write(*,*) "Age resolution: ", age_resolution, " Ga." + do i = 1,MAXAGEBINS + domain%age_bin_times(i) = maxageGa-(i*age_resolution) + end do end if - age_resolution = maxageGa / real(MAXAGEBINS) - write(*,*) "Age resolution: ", age_resolution, " Ga." - do i = 1,MAXAGEBINS - domain%age_bin_times(i) = maxageGa-(i*age_resolution) - end do domain%age_counter = 1 oldGa = 0._DP @@ -175,25 +181,27 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt timestamp_old = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=DP) icrater = icrater + 1 crater%timestamp = real(curyear + real(icrater,kind=DP) / real(ntotcrat,kind=DP) * user%interval,kind=DP) - if (icrater .eq. 1) then - agemin = crater%timestamp * 0.9_DP - end if - if (crater%timestamp < 2330._DP) then - if (oldGa > 0._DP) then - if (user%numintervals .eq. 1) then - crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,agemin,oldGa) + if (user%doregotrack) then + if (icrater .eq. 1) then + agemin = crater%timestamp * 0.9_DP + end if + if (crater%timestamp < 2330._DP) then + if (oldGa > 0._DP) then + if ((user%numintervals .eq. 1) .or. (user%runtype .eq. 'STATISTICAL')) then + crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,agemin,oldGa) + else + crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,1e-10_DP,oldGa) + end if else - crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,1e-10_DP,oldGa) + if ((user%numintervals .eq. 1) .or. (user%runtype .eq. 'STATISTICAL')) then + crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,agemin,maxageGa) + else + crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,1e-10_DP,maxageGa) + end if end if else - if (user%numintervals .eq. 1) then - crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,agemin,maxageGa) - else - crater%timestampGa = util_t_from_scale(maxage-crater%timestamp,1e-10_DP,maxageGa) - end if + crater%timestampGa = 4.5_DP end if - else - crater%timestampGa = 4.5_DP end if pbarpos = nint(real(icrater) / real(ntotcrat) * PBARRES) if (crater%timestampGa < domain%age_bin_times(domain%age_counter)) then @@ -307,10 +315,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt if (crater%ejdis > domain%smallest_ejecta) then ! Estimated size is big enough, so proceed with precise calculation if (user%doregotrack) then call ejecta_table_define(user,crater,domain,ejb,ejtble,melt) - !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) else call ejecta_table_define(user,crater,domain,ejb,ejtble) - !call ejecta_interpolate(crater,domain,crater%frad,ejb(1:ejtble),ejtble,crater%ejrim) end if call ejecta_emplace(user,surf,crater,domain,ejb(1:ejtble),ejtble,ejbmass,& ejecta_dem,nmeltsheet,vmeltsheet) diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index 3da7b696..63a69bfe 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -45,12 +45,20 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) real(DP),dimension(3) :: rn real(DP) :: superlen,rayfrac,cutout,fe,fd type(cratertype) :: crater - real(DP) :: eji, diffi + real(DP),dimension(:,:),allocatable :: diffdistribution,ejdistribution integer(I4B),dimension(:,:),allocatable :: ejisray real(DP) :: xbar,ybar,dD,xp,yp,areafrac,krad,supersize,lrad 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 = 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 + ! Create box for soften calculation (will be no bigger than the grid itself) do j = 0,user%gridsize + 1 do i = 0,user%gridsize + 1 @@ -174,6 +182,12 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) jmax = ypi - crater%ylpx + allocate(diffdistribution(imin:imax,jmin:jmax)) + allocate(ejdistribution(imin:imax,jmin:jmax)) + rray = 11.95_DP*crater%frad**1.32 + l1 = 5.32_DP*crater%frad**1.27 + + call ejecta_ray_pattern(user,surf,crater,inc,imin,imax,jmin,jmax,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) & @@ -182,17 +196,17 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) do i = imin,imax xpi = crater%xlpx + i ypi = crater%ylpx + j - call ejecta_ray_pattern(user,crater,i,j,diffi,eji) - kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffi + kdiff(xpi,ypi) = kdiff(xpi,ypi) + dKdN * diffdistribution(i,j) !TEMP xp = xpi * user%pix yp = ypi * user%pix lrad = sqrt((xp - crater%xl)**2 + (yp - crater%yl)**2) - surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + eji & + surf(xpi,ypi)%ejcov = surf(xpi,ypi)%ejcov + ejdistribution(i,j) & * 0.14_DP * crater%frad**(0.74_DP) * (lrad / crater%frad)**(-3) end do end do !!$OMP END PARALLEL DO + deallocate(diffdistribution,ejdistribution) end do end if diff --git a/src/crater/crater_superdomain.f90 b/src/crater/crater_superdomain.f90 index 997bd473..65c94e09 100644 --- a/src/crater/crater_superdomain.f90 +++ b/src/crater/crater_superdomain.f90 @@ -42,8 +42,8 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) real(DP),dimension(2) :: rn real(DP) :: superlen,rayfrac type(cratertype) :: crater - real(DP) :: diffi, eji - logical :: ejisray + real(DP),dimension(:,:),allocatable :: ejdistribution, diffdistribution + integer(I4B),dimension(:,:),allocatable :: ejisray ! Melt or glassy ray test real(DP) :: xp, yp, lradsq, lrad, erad @@ -51,6 +51,18 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) real(DP) :: rm, depthb, maxgcrat, maxgcratkm real(SP) :: gglass, glrad + ! Crater ray parameters + 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 do i = 0,user%gridsize + 1 @@ -144,26 +156,40 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) lrad = minval(lradif) if (lrad > crater%ejdis) cycle + allocate(ejdistribution(xi:xf,yi:yf)) + 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,l1) ! Now if doregotrack is on, do melt zone calculation if (user%doregotrack) call regolith_melt_zone_superdomain(user,crater,domain,rm,depthb) + do j = yi,yf + do i = xi,xf + ! Ejecta ray distribution + if (ejdistribution(i,j) > 0.0001_DP) then + ejisray(i,j) = 1 + else + ejisray(i,j) = 0 + end if + end do + end do + if (user%doregotrack) then do j = 1, user%gridsize do i = 1, user%gridsize xpi = i - crater%xlpx ypi = j - crater%ylpx if ((abs(xpi) > inc) .or. (abs(ypi) > inc)) cycle - call ejecta_ray_pattern(user,crater,i,j,diffi,eji) - ejisray = (eji > 0.0001_DP) - if (.not.ejisray) cycle - call regolith_superdomain(user,crater,domain,surf(i,j)%regolayer,eji,& + if (ejisray(xpi,ypi) == 0) cycle + call regolith_superdomain(user,crater,domain,surf(i,j)%regolayer,ejdistribution(xpi,ypi),& i,j,rm,depthb) end do end do end if + deallocate(ejdistribution,diffdistribution,ejisray) end do diff --git a/src/crater/module_crater.f90 b/src/crater/module_crater.f90 index 1c8e9253..d5108399 100644 --- a/src/crater/module_crater.f90 +++ b/src/crater/module_crater.f90 @@ -338,62 +338,4 @@ subroutine crater_superdomain(user,surf,prod,nflux,domain,finterval) end subroutine crater_superdomain end interface - - - -! new subroutines added by jundu on 10/15/2022 - - -! Rim_crest: - - ! interface - - ! subroutine Calculate_am_wl_phase_from_diameter(psd_1D,amplitude,wavelength,phase) - ! use module_globals - ! implicit none - ! ! in and out - ! type(psdtype),intent(inout) :: psd_1D - ! real(DP),dimension(:), allocatable,intent(out) :: amplitude,wavelength,phase - ! end subroutine Calculate_am_wl_phase_from_diameter - - ! subroutine Calculate_breakpoint_slope_from_diameter(psd_1D) - ! use module_globals - ! implicit none - ! ! in and out - ! type(psdtype),intent(inout) :: psd_1D - ! end subroutine Calculate_breakpoint_slope_from_diameter - - ! subroutine Calculate_targetPSD_from_breakpoint_slope(psd_1D,wavelength,psd) - ! use module_globals - ! implicit none - ! !in and out - ! type(psdtype), intent(in) :: psd_1D - ! real(DP) ,dimension(:),allocatable,intent(out) :: wavelength,psd - ! end subroutine Calculate_targetPSD_from_breakpoint_slope - - ! subroutine Calculate_am_wl_phase_from_targetPSD(psd_1D,wavelength,psd,amplitude,phase) - ! use module_globals - ! implicit none - ! !in and out - ! type(psdtype), intent(in) :: psd_1D - ! real(DP) ,dimension(:),intent(in) :: wavelength,psd - ! real(DP) ,dimension(:),allocatable,intent(out) :: amplitude,phase - ! end subroutine Calculate_am_wl_phase_from_targetPSD - - - ! subroutine Create_rim(arc_length,psd_1D,amplitude,wavelength,phase,rim_parameter) - ! use module_globals - ! implicit none - ! ! in and out - ! real(DP),intent(in) :: arc_length - ! type(psdtype), intent(in) :: psd_1D - ! real(DP),dimension(:),intent(in) :: amplitude - ! real(DP),dimension(:),intent(in) :: wavelength - ! real(DP),dimension(:),intent(in) :: phase - ! real(DP),intent(out) :: rim_parameter - ! end subroutine Create_rim - - ! end interface - - end module diff --git a/src/ejecta/ejecta_emplace.f90 b/src/ejecta/ejecta_emplace.f90 index 737be2b3..5186f6cd 100644 --- a/src/ejecta/ejecta_emplace.f90 +++ b/src/ejecta/ejecta_emplace.f90 @@ -1,4 +1,4 @@ -!***** ejecta/ejecta_emplace +!****f* ejecta/ejecta_emplace ! Name ! ejecta_emplace -- Calculate ejecta mass during excavation stage. ! SYNOPSIS @@ -100,18 +100,25 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ ! Internal variables real(DP) :: lrad,lradsq integer(I4B),parameter :: MAXLOOP = 100 ! Maximum number of times to loop the ejecta angle correction calculation - integer(I4B) :: xpi,ypi,i,j,n,inc,incsq,iradsq,idistorted,jdistorted + integer(I4B) :: xpi,ypi,i,j,k,n,inc,incsq,iradsq,idistorted,jdistorted real(DP) :: xp,yp,fradsq,fradpxsq,radsq,ebh,ejdissq,ejbmass,fmasscons,areafrac,xbar,ybar,krad,kdiffmax - real(DP),dimension(:,:),allocatable :: kdiff,cel - integer(I4B),dimension(:,:,:),allocatable :: indarray - integer(I4B) :: maxhits,nin,nnot,dradsq + real(DP),dimension(:,:),allocatable :: big_cumulative_elchange,kdiff,big_kdiff,cel,big_cel + integer(I4B),dimension(:,:,:),allocatable :: indarray,big_indarray + real(DP),dimension(:,:),allocatable :: ejdistribution,diffdistribution,maxdiff, maxej + real(DP),dimension(:,:),allocatable :: tempdiff,tempej + integer(I4B) :: bigi,bigj,maxhits,nin,nnot,dradsq character(len=MESSAGESIZE) :: message ! message for the progress bar real(DP) :: vmelt, totmelt, volm - real(DP) :: diffi,eji - logical :: bigej + real(DP) :: frayreduction = 0.5_DP ! Factor to apply to reduce the relative thickness of the ray for each subsequent pattern + integer(I4B), parameter :: Npatt = 8 ! Number of times to call ray pattern + + + + ! Ray mixing model variables + real(DP) :: dsc ! Melt zone's radius - real(DP) :: rm, dm, melt + real(DP) :: rm, dm, melt, eradc ! Ejecta pattern distortion parameters real(DP) :: distance,erad,craterslope,landslope,baseline,lrange,frac,ejheight,ebh0,maxdistance @@ -119,6 +126,22 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ real(DP) :: vsq, ejtheta integer(I4B) :: ind,klo + ! Age + real(SP) :: age_mean + + ! Crater ray parameters + real(DP) :: rray != 48.0_DP ! "L16" in Minton et al. (2019) + integer(I4B) :: Nraymax = 5 + 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 + + + l1 = (5.32_DP*(crater%frad/1000)**1.27)/(crater%frad/1000) + rray = user%ejecta_truncation + ! Executable code if (user%doregotrack) call regolith_melt_zone(user,crater,crater%imp,crater%impvel,rm,dm,totmelt) @@ -152,39 +175,56 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ incsq = inc**2 - bigej = (inc >= user%gridsize / 2) - if (bigej) then - allocate(cumulative_elchange(0:user%gridsize+1,0:user%gridsize+1)) - allocate(cel(0:user%gridsize+1,0:user%gridsize+1)) - allocate(kdiff(0:user%gridsize+1,0:user%gridsize+1)) - allocate(indarray(2,0:user%gridsize+1,0:user%gridsize+1)) - cumulative_elchange(:,:) = 0.0_DP - cel(:,:) = 0.0_DP - kdiff(:,:) = 0.0_DP - + if (inc >= user%gridsize / 2) then if (user%testflag) then - write(*,*) 'Big ejecta: fcrat =',crater%fcrat, ' Ej/S =',(crater%ejdispx*user%pix)/domain%side, ' Ejrim =', crater%ejrim - else + write(*,*) 'Big ejecta: fcrat =',crater%fcrat, ' Ej/S =',(crater%ejdispx*user%pix)/domain%side, ' Ejrim =', crater%ejrim + write(*,*) 'Rray = ', rray, 'L1 = ', l1 + else write(message,'("Ejb: Dc=",ES9.2," Ej/S=",F0.3)') crater%fcrat,(crater%ejdispx*user%pix)/domain%side call io_updatePbar(message) - end if - else - allocate(cumulative_elchange(-inc:inc,-inc:inc)) - allocate(cel(-inc:inc,-inc:inc)) - allocate(kdiff(-inc:inc,-inc:inc)) - allocate(indarray(2,-inc:inc,-inc:inc)) - end if + end if + endif + + allocate(ejdistribution(-inc:inc,-inc:inc)) + allocate(diffdistribution(-inc:inc,-inc:inc)) + allocate(tempdiff(-inc:inc,-inc:inc)) + allocate(tempej(-inc:inc,-inc:inc)) + + ejdistribution(:,:) = 0.0_DP + diffdistribution(:,:) = 0.0_DP + + ! *************************** 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,l1) + diffdistribution(:,:) = diffdistribution(:,:) + frayreduction**(i-1) * tempdiff(:,:) + ejdistribution(:,:) = ejdistribution(:,:) + frayreduction**(i-1) * tempej(:,:) + end do + + diffdistribution(:,:) = diffdistribution(:,:) / maxval(ejdistribution) + ejdistribution(:,:) = ejdistribution(:,:) / maxval(diffdistribution) + + deallocate(tempdiff,tempej) + ! *****************************************************************************! + + allocate(cumulative_elchange(-inc:inc,-inc:inc)) + allocate(cel(-inc:inc,-inc:inc)) + allocate(kdiff(-inc:inc,-inc:inc)) + allocate(indarray(2,-inc:inc,-inc:inc)) cumulative_elchange = 0.0_DP kdiff = 0.0_DP indarray = inc - 1 ! initialize this array to point to a corner (this should have 0 elevation change since we're only doing work ! within a circle of radius irad - ejbmass = 0.0_DP nin = 0 nnot = 0 + !!$OMP PARALLEL DO DEFAULT(PRIVATE) IF(inc > INCPAR) & + !!$OMP SHARED(user,domain,crater,surf,ejb,ejtble) & + !!$OMP SHARED(inc,incsq) & + !!$OMP SHARED(cumulative_elchange,kdiff,kdiffmax,indarray,ejdistribution,diffdistribution) + !open(74, file='meltvserad.csv', status='replace') do j = -inc,inc do i = -inc,inc ! find distance from crater center @@ -199,10 +239,8 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) - if (.not.bigej) then - indarray(1,i,j) = xpi - indarray(2,i,j) = ypi - end if + indarray(1,i,j) = xpi + indarray(2,i,j) = ypi lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 lrad = sqrt(lradsq) @@ -214,11 +252,17 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ maxslp = -huge(maxslp) klo = int((log(lrad) - log(crater%ejrad)) / domain%ejbres) do n = 1,MAXLOOP - call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) - if ((n > 1).and.((abs(ebh0 - ebh) / ebh0) < domain%small)) exit + if (user%doregotrack) then + call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad,melt=melt) + else + call ejecta_interpolate(crater,domain,distance,ejb,ejtble,ebh,vsq=vsq,theta=ejtheta,erad=erad) + end if + if (ebh < VSMALL) exit + if (n > 1) then + if ((abs(ebh0 - ebh) / ebh0) < domain%small) exit + endif ebh0 = ebh - erad = exp(erad) lrange = lrad - erad baseline = ((i * crater%xslp) + (j * crater%yslp)) * user%pix @@ -237,77 +281,61 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ ! Find out where in the table this new velocity corresponds to ind = 1 - call util_search(ejb%vesq,ind,ejtble,vsq,klo) + call util_search(ejb%vesq,ind,ejtble,log(vsq),klo) klo = min(max(klo,1),ejtble-1) ! Interpolate on the table to find the flat plane equivalent landing distance for this velocity - frac = (vsq - ejb(klo)%vesq) / (ejb(klo+1)%vesq - ejb(klo)%vesq) + frac = (vsq - exp(ejb(klo)%vesq)) / (exp(ejb(klo+1)%vesq) - exp(ejb(klo)%vesq)) distance = exp(ejb(klo)%lrad) + frac * (exp(ejb(klo+1)%lrad) - exp(ejb(klo)%lrad)) end do if (vsq < 0.0_DP) cycle if (distance /= distance) cycle - idistorted = int(i * distance / lrad) if (abs(idistorted) > inc) cycle jdistorted = int(j * distance / lrad) if (abs(jdistorted) > inc) cycle - + iradsq = idistorted**2 + jdistorted**2 if ((iradsq > incsq).or.(distance <= crater%ejrad)) cycle - call ejecta_ray_pattern(user,crater,idistorted,jdistorted,diffi,eji) ! we need to cut a hole out from the inside of the crater xbar = xpi * user%pix - crater%xl ybar = ypi * user%pix - crater%yl areafrac = (1.0_DP - util_area_intersection(crater%ejrad,xbar,ybar,user%pix)) - ebh = areafrac * eji * ebh - if (bigej) then - cumulative_elchange(xpi,ypi) = cumulative_elchange(xpi,ypi) + ebh + crater_profile(user, crater, lrad) - else - cumulative_elchange(i,j) = ebh + crater_profile(user, crater, lrad) - end if + ebh = areafrac * ejdistribution(idistorted,jdistorted) * ebh + cumulative_elchange(i,j) = ebh + crater_profile(user, crater, lrad) if (user%dosoftening) then ! Do extra diffusive degradation over ejecta region areafrac = (1.0_DP - util_area_intersection(crater%frad,xbar,ybar,user%pix)) areafrac = areafrac * util_area_intersection(crater%fe * crater%frad,xbar,ybar,user%pix) - if (bigej) then - kdiff(xpi,ypi) = kdiff(xpi,ypi) + areafrac * diffi * kdiffmax - else - kdiff(i,j) = areafrac * diffi * kdiffmax - end if - + kdiff(i,j) = areafrac * diffdistribution(idistorted,jdistorted) * kdiffmax end if end do end do + !close(74) + !!$OMP END PARALLEL DO + ! if(user%doregotrack .and. user%testflag) then + ! write(*,*) 'Ejected Melt: ', vmelt + ! write(*,*) 'Total Melt: ', totmelt + ! write(*,*) 'ejected / total melt:', vmelt/totmelt + ! end if ejbmass = sum(cumulative_elchange) ! Create buffer to prevent infinite hole bug - if (bigej) then - kdiff(0,:) = 0.0_DP - kdiff(user%gridsize+1,:) = 0.0_DP - kdiff(:,0) = 0.0_DP - kdiff(:,user%gridsize+1) = 0.0_DP - - cumulative_elchange(0,:) = 0.0_DP - cumulative_elchange(user%gridsize+1,:) = 0.0_DP - cumulative_elchange(:,0) = 0.0_DP - cumulative_elchange(:,user%gridsize+1) = 0.0_DP - else - kdiff(-inc,:) = 0.0_DP - kdiff(inc,:) = 0.0_DP - kdiff(:,-inc) = 0.0_DP - kdiff(:,inc) = 0.0_DP - - cumulative_elchange(-inc,:) = 0.0_DP - cumulative_elchange(inc,:) = 0.0_DP - cumulative_elchange(:,-inc) = 0.0_DP - cumulative_elchange(:,inc) = 0.0_DP - end if + kdiff(-inc,:) = 0.0_DP + kdiff(inc,:) = 0.0_DP + kdiff(:,-inc) = 0.0_DP + kdiff(:,inc) = 0.0_DP + + cumulative_elchange(-inc,:) = 0.0_DP + cumulative_elchange(inc,:) = 0.0_DP + cumulative_elchange(:,-inc) = 0.0_DP + cumulative_elchange(:,inc) = 0.0_DP ! Do mass conservation by adjusting ejecta thickness fmasscons = (-deltaMtot)/ ejbmass @@ -331,67 +359,42 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ ! periodic boundary conditions call util_periodic(xpi,ypi,user%gridsize) - - if (.not.bigej) then - indarray(1,i,j) = xpi - indarray(2,i,j) = ypi - end if + + indarray(1,i,j) = xpi + indarray(2,i,j) = ypi lradsq = (crater%xl - xp)**2 + (crater%yl - yp)**2 lrad = sqrt(lradsq) if (lrad < crater%ejrad) cycle - - if (bigej) then - ebh = cumulative_elchange(xpi,ypi) - crater_profile(user, crater, lrad) - else - ebh = cumulative_elchange(i,j) - crater_profile(user, crater, lrad) - end if - - if (user%doregotrack .and. ebh>1.0e-8_DP) then - call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm) - vmelt = vmelt + volm - end if + + + ebh = cumulative_elchange(i,j) - crater_profile(user, crater, lrad) + + + if (user%doregotrack .and. ebh>1.0e-8_DP) then + call regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm) + vmelt = vmelt + volm + !write(74,*) erad, surf(xpi,ypi)%regolayer%regodata%meltfrac + end if end do end do end if - if (totmelt > vmelt) then - vmeltsheet = totmelt - vmelt - else !give the craters a melt sheet of 1m - vmeltsheet = 1.0_DP * user%pix * user%pix * nmeltsheet + if (user%doregotrack) then + if (totmelt > vmelt) then + vmeltsheet = totmelt - vmelt + else !give the craters a melt sheet of 1m + vmeltsheet = 1.0_DP * user%pix * user%pix * nmeltsheet + end if end if - ! extra soften calculation - if (user%dosoftening) then - cel = 0.0_DP - - if (bigej) then - call util_diffusion_solver(user,surf,user%gridsize + 2,indarray,kdiff,cel,maxhits) - do ypi = 1,user%gridsize - do xpi = 1,user%gridsize - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cel(xpi,ypi) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cel(xpi,ypi),0.0_DP) - indarray(1,xpi,ypi) = xpi - indarray(2,xpi,ypi) = ypi - end do - end do - indarray(1:2,0,:) = user%gridsize - indarray(1:2,user%gridsize+1,:) = 1 - indarray(1:2,:,0) = user%gridsize - indarray(1:2,:,user%gridsize+1) = 1 - - call ejecta_soften(user,surf,user%gridsize + 2,indarray,cumulative_elchange) - - ! Add the ejecta back to the DEM - do ypi = 1,user%gridsize - do xpi = 1,user%gridsize - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(xpi,ypi) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(xpi,ypi), 0.0_DP) - end do - end do + ! Create box for soften calculation (will be no bigger than the grid itself) + if (2 * inc + 1 < user%gridsize) then - else + ! extra soften calculation + if (user%dosoftening) then + cel = 0.0_DP call util_diffusion_solver(user,surf,2 * inc + 1,indarray,kdiff,cel,maxhits) do j = -inc,inc do i = -inc,inc @@ -401,22 +404,75 @@ subroutine ejecta_emplace(user,surf,crater,domain,ejb,ejtble,deltaMtot,cumulativ surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cel(i,j),0.0_DP) end do end do + end if - call ejecta_soften(user,surf,2 * inc + 1,indarray,cumulative_elchange) + call ejecta_soften(user,surf,2 * inc + 1,indarray,cumulative_elchange) - ! Add the ejecta back to the DEM - do j = -inc,inc - do i = -inc,inc - xpi = indarray(1,i,j) - ypi = indarray(2,i,j) - surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j) - surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j), 0.0_DP) + ! Add the ejecta back to the DEM + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + surf(xpi,ypi)%dem = surf(xpi,ypi)%dem + cumulative_elchange(i,j) + surf(xpi,ypi)%ejcov = max(surf(xpi,ypi)%ejcov + cumulative_elchange(i,j), 0.0_DP) + end do + end do + + + else ! Ejecta wraps around the grid. + ! We will therefore send in the whole grid with the total ejecta thickness added to each pixel + allocate(big_cumulative_elchange(0:user%gridsize+1,0:user%gridsize+1)) + allocate(big_cel(0:user%gridsize+1,0:user%gridsize+1)) + allocate(big_indarray(2,0:user%gridsize+1,0:user%gridsize+1)) + allocate(big_kdiff(0:user%gridsize+1,0:user%gridsize+1)) + + do bigj = 0,user%gridsize + 1 + do bigi = 0,user%gridsize + 1 + xpi = bigi + ypi = bigj + call util_periodic(xpi,ypi,user%gridsize) + big_indarray(1,bigi,bigj) = xpi + big_indarray(2,bigi,bigj) = ypi + big_cumulative_elchange(bigi,bigj) = 0.0_DP + end do + end do + + big_kdiff = 0.0_DP + + do j = -inc,inc + do i = -inc,inc + xpi = indarray(1,i,j) + ypi = indarray(2,i,j) + big_cumulative_elchange(xpi,ypi) = big_cumulative_elchange(xpi,ypi) + cumulative_elchange(i,j) + big_kdiff(xpi,ypi) = big_kdiff(xpi,ypi) + kdiff(i,j) + end do + end do + + if (user%dosoftening) then + big_cel = 0.0_DP + call util_diffusion_solver(user,surf,user%gridsize + 2,big_indarray,big_kdiff,big_cel,maxhits) + + do j = 1,user%gridsize + do i = 1,user%gridsize + surf(i,j)%dem = surf(i,j)%dem + big_cel(i,j) + surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cel(i,j),0.0_DP) end do end do end if + + call ejecta_soften(user,surf,user%gridsize + 2,big_indarray,big_cumulative_elchange) + + do i = 1,user%gridsize + do j = 1,user%gridsize + surf(i,j)%dem = surf(i,j)%dem + big_cumulative_elchange(i,j) + surf(i,j)%ejcov = max(surf(i,j)%ejcov + big_cumulative_elchange(i,j),0.0_DP) + end do + end do + + deallocate(big_cumulative_elchange,big_indarray,big_kdiff,big_cel) end if - deallocate(indarray,kdiff,cel) + deallocate(indarray,diffdistribution,ejdistribution,kdiff,cel) return end subroutine ejecta_emplace diff --git a/src/ejecta/ejecta_interpolate.f90 b/src/ejecta/ejecta_interpolate.f90 index c71d50a0..1b181ee2 100644 --- a/src/ejecta/ejecta_interpolate.f90 +++ b/src/ejecta/ejecta_interpolate.f90 @@ -69,9 +69,48 @@ subroutine ejecta_interpolate(crater,domain,lrad,ejb,ejtble,ebh,vsq,theta,erad,m if (present(melt)) melt = ejb(k)%meltfrac - ((ejb(k)%meltfrac - ejb(k+1)%meltfrac) * frac) if (present(erad)) erad = ejb(k)%erad - ((ejb(k)%erad - ejb(k+1)%erad) * frac) end if - ebh = exp(ebh) - if (lrad > crater%ejdis) ebh = 0._DP - + if (ebh < LOGVSMALL) then + ebh = 0.0_DP + else + ebh = exp(ebh) + end if + if (lrad > crater%ejdis) then + ebh = 0._DP + if (present(vsq)) vsq = 0._DP + if (present(theta)) theta = 0._DP + if (present(melt)) melt = 0._DP + if (present(erad)) erad = 0._DP + else + if (present(vsq)) then + if (vsq < LOGVSMALL) then + vsq = 0.0_DP + else + vsq = exp(vsq) + end if + end if + if (present(theta)) then + if (theta < LOGVSMALL) then + theta = 0.0_DP + else + theta = exp(theta) + end if + end if + if (present(melt)) then + if (melt < LOGVSMALL) then + melt = 0.0_DP + else + melt = exp(melt) + end if + end if + if (present(erad)) then + if (erad < LOGVSMALL) then + erad = 0.0_DP + else + erad = exp(erad) + end if + end if + end if + return end subroutine ejecta_interpolate 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 diff --git a/src/ejecta/ejecta_ray_pattern.f90 b/src/ejecta/ejecta_ray_pattern.f90 index 663a216d..41f26bbe 100644 --- a/src/ejecta/ejecta_ray_pattern.f90 +++ b/src/ejecta/ejecta_ray_pattern.f90 @@ -26,11 +26,12 @@ ! * domain -- Simulation domain variable container ! ! Output +! * ejdist -- logical array containing the pixels that contain ejecta ! !*** !********************************************************************************************************************************** -subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) +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 @@ -41,22 +42,27 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) ! Arguments type(usertype),intent(in) :: user + type(surftype),dimension(:,:),intent(in) :: surf type(cratertype),intent(inout) :: crater - integer(I4B),intent(in) :: i,j - real(DP),intent(out) :: diffi - real(DP),intent(out) :: eji + integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),intent(in) :: rray, fpeak, rayp, rayfmult + 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,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq + integer(I4B) :: nrays,i,j,k,n,nef,incsq,iradsq,xpi,ypi,ejpxsq real(DP) :: frac,mef,lrad,lradsq,xp,yp,binres,areafrac,xbar,ybar real(DP) :: rn real(DP) :: theta, lradp, maxdistance real(DP), parameter :: n1 = 4.0_DP real(DP) :: n2, mag + real(DP),dimension(xi:xf,yi:yf) :: isray real(DP),dimension(:),allocatable :: numinray,totnum real(DP),dimension(:),allocatable :: mefarray real(DP) :: ans - logical :: bigej + real(DP) :: C1,C2,p ! Fits to fe vs fd equation for the new ray pattern real(DP) :: rmin,rmax,r @@ -71,51 +77,79 @@ subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) logical :: ej real(DP),dimension(Nraymax) :: thetari - eji = 0.0_DP - diffi = 0.0_DP + !TEMPORARY - if (user%dorays) then - do k = 1,Nraymax - thetari(k) = 2 * pi * k / Nraymax + if (user%dorays .and. crater%fcrat<1500000._DP) then + do i = 1,Nraymax + thetari(i) = 2 * pi * i / Nraymax end do call shuffle(thetari) ! randomize the ray pattern call random_number(rn) ! randomize the ray orientation - rmax = user%ejecta_truncation + rmax = rray rmin = crater%continuous / crater%frad + !rmax = rray / crater%frad + !rmin = l1 / crater%frad crater%fe = 10.0_DP ! Estimate the equivalent degradation radius - xpi = crater%xlpx + i - ypi = crater%ylpx + j - - ! Find distance from crater center to current pixel center in real space - xp = xpi * user%pix - yp = ypi * user%pix - - xbar = xp - crater%xl - ybar = yp - crater%yl - - 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) - diffi = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.false.) - eji = areafrac * ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,.true.) + !ejdistribution = 0.0_DP + !diffdistribution = 0.0_DP + !!$OMP PARALLEL DO DEFAULT(PRIVATE) & + !!$OMP SHARED(user,crater) & + !!$OMP SHARED(xi,xf,yi,yf,rn,diffdistribution,ejdistribution,thetari,rmin) + do j = yi,yf + do i = xi,xf + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! Find distance from crater center to current pixel center in real space + xp = xpi * user%pix + yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl + + areafrac = util_area_intersection(rray * 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,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 + else + !Do simple circular region + incsq = inc**2 + ejdistribution = 0.0_DP + diffdistribution = 0.0_DP + !$OMP PARALLEL DO DEFAULT(PRIVATE) & + !$OMP SHARED(user,crater) & + !$OMP SHARED(inc,incsq,xi,xf,yi,yf,diffdistribution,ejdistribution) + do j = yi,yf + do i = xi,xf + iradsq = i*i + j*j + + if (iradsq < incsq) then + + xpi = crater%xlpx + i + ypi = crater%ylpx + j + + ! Find distance from crater center to current pixel center in real space + xp = xpi * user%pix + yp = ypi * user%pix + + xbar = xp - crater%xl + ybar = yp - crater%yl + areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) ! uniform circular + diffdistribution(i,j) = areafrac + ejdistribution(i,j) = areafrac + end if + end do + end do + !$OMP END PARALLEL DO - xpi = crater%xlpx + i - ypi = crater%ylpx + j - - ! Find distance from crater center to current pixel center in real space - xp = xpi * user%pix - yp = ypi * user%pix - - xbar = xp - crater%xl - ybar = yp - crater%yl - areafrac = util_area_intersection(user%ejecta_truncation * crater%frad,xbar,ybar,user%pix) ! uniform circular - diffi = areafrac - eji = areafrac end if - return contains diff --git a/src/ejecta/ejecta_ray_pattern_func.f90 b/src/ejecta/ejecta_ray_pattern_func.f90 index 8f83a704..7843a2e7 100644 --- a/src/ejecta/ejecta_ray_pattern_func.f90 +++ b/src/ejecta/ejecta_ray_pattern_func.f90 @@ -27,21 +27,27 @@ !*** !********************************************************************************************************** -function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,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 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 real(DP) :: a,c real(DP) :: thetar,rw,rw0,rw1 real(DP) :: f,rtrans,length,rpeak,minray,FF integer(I4B) :: n,i + real(DP) :: tmp - minray = rmin * 3 + !minray = rmin * 3 !"L1" in Minton et al. (2019) + !minray = rray / 2 + minray = l1 if (r > rmax) then ans = 0._DP @@ -55,18 +61,23 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans) rw0 = rmin * pi / Nraymax / 2 rw1 = 2 * pi / Nraymax rw = rw0 * (1._DP - (1.0_DP - rw1 / rw0) * exp(1._DP - (r / rmin)**2)) ! equation 40 Minton et al. 2019 - 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 + tmp = (Nraymax**rayp - (Nraymax**rayp - 1) * log(r/minray) / log(rray/minray)) + if (tmp < 0.0_DP) then + n = Nraymax ! "Nrays" in Minton et al. (2019) + else + 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 + end if ans = 0._DP rtrans = r - 1.0_DP c = rw / r a = sqrt(2 * pi) / (n * c * erf(pi / (2 *sqrt(2._DP) * c))) !equation 39 Minton et al., 2019 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 + length = minray * exp(log(rray/minray) * ((Nraymax - i + 1)**rayp - 1.0_DP) / ((Nraymax**rayp - 1))) + rpeak = (length - 1.0_DP) * 0.5_DP if (ej) then FF = 1.0_DP if (r > length) then - f = 0.0_DP + cycle ! Don't add any material beyond the length of the ray else f = a end if @@ -74,7 +85,8 @@ function ejecta_ray_pattern_func(theta,r,rmin,rmax,thetari,ej) result(ans) 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)) !equation 42 Minton et al. 2019 end if - ans = ans + ejecta_ray_func(theta,thetari(i),r,n,rw) * f / a + tmp = ejecta_ray_func(theta,thetari(i),r,n,rw) + if (tmp > epsilon(ans) .and. (f/a > epsilon(ans))) ans = ans + tmp * f / a ! Ensure that we don't get an underflow end do end if diff --git a/src/ejecta/ejecta_table_define.f90 b/src/ejecta/ejecta_table_define.f90 index cf259e2d..408d9c7e 100644 --- a/src/ejecta/ejecta_table_define.f90 +++ b/src/ejecta/ejecta_table_define.f90 @@ -80,12 +80,12 @@ subroutine ejecta_table_define(user,crater,domain,ejb,ejtble,melt) thick = max(crater_profile(user,crater,r),VSMALL) end if ejb(k)%thick = log(thick) - ejb(k)%vesq = vejsq - ejb(k)%angle = ejang + ejb(k)%vesq = log(vejsq) + ejb(k)%angle = log(ejang) ejb(k)%erad = log(erad) if (present(melt)) then call regolith_melt_fraction(dimp,depthb,erad,eradold,rmelt,melt) - ejb(k)%meltfrac = melt + ejb(k)%meltfrac = log(melt) end if if ((thick <= VSMALL) .or. (abs(eradold - erad) < VSMALL)) then ejtble = k diff --git a/src/ejecta/module_ejecta.f90 b/src/ejecta/module_ejecta.f90 index 9e75eb4d..e2d8747a 100644 --- a/src/ejecta/module_ejecta.f90 +++ b/src/ejecta/module_ejecta.f90 @@ -43,24 +43,32 @@ end subroutine ejecta_emplace end interface interface - subroutine ejecta_ray_pattern(user,crater,i,j,diffi,eji) + 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 + type(surftype),dimension(:,:),intent(in) :: surf type(cratertype),intent(inout) :: crater - integer(I4B),intent(in) :: i,j - real(DP),intent(out) :: diffi - real(DP),intent(out) :: eji + integer(I4B),intent(in) :: inc,xi,xf,yi,yf + real(DP),intent(in) :: rray, fpeak, rayp, rayfmult + 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,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 diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 3ce3f5cd..a1d38273 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -45,7 +45,7 @@ module module_globals integer(I4B), parameter :: UPPERCASE_OFFSET = iachar('A') - iachar('a') ! Miscellaneous constants: -real(DP),parameter :: VSMALL = tiny(1._DP) ! Very small number +real(DP),parameter :: VSMALL = 10*tiny(1._DP) ! Very small number real(DP),parameter :: LOGVSMALL = log(VSMALL) ! log of a very small number real(DP),parameter :: VBIG = huge(1._DP) ! Very big number real(DP),parameter :: SMALLFAC = 1e-5_DP ! Smallest unit of measurement proportional to pixel size @@ -221,7 +221,7 @@ module module_globals ! Ejecta softening variables logical :: dosoftening ! Set T to use the extra crater softening model - real(DP) :: ejecta_truncation ! Set the number of crater diameters to truncate the ejecta + real(DP) :: ejecta_truncation ! Set the number of crater radii to truncate the ejecta logical :: dorays ! Set T to use ray model logical :: superdomain ! Set T to include the superdomain @@ -383,12 +383,4 @@ module module_globals real(DP),parameter :: GFAC = 0.487_DP ! gravitational acceleration exponent real(DP),parameter :: DFAC = 0.556_DP ! impact distance exponent -! Crater ray parameters -real(DP),parameter :: rray = 16_DP -integer(I4B),parameter :: Nraymax = 16 -real(DP),parameter :: fpeak = 8000_DP ! narrow ray: rw0 propto 1/4 -real(DP),parameter :: rayp = 2.0_DP -integer(I4B),parameter :: rayq = 4 -real(DP),parameter :: rayfmult = (5)**(-4.0_DP / (1.2_DP)) - end module module_globals diff --git a/src/init/init_dist.f90 b/src/init/init_dist.f90 index fc17a1bd..72e96af5 100644 --- a/src/init/init_dist.f90 +++ b/src/init/init_dist.f90 @@ -33,7 +33,7 @@ subroutine init_dist(user,domain) ! Executable code ! Get size of crater production function array - open(unit=LUN,file=user%sfdfile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%sfdfile)),status="old",iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ", trim(user%sfdfile) stop @@ -47,12 +47,12 @@ subroutine init_dist(user,domain) end do if (domain%pnum == 0) then - write(*,*) "No valid entries in ",trim(user%sfdfile) + write(*,*) "No valid entries in ",trim(adjustl(user%sfdfile)) end if close(LUN) ! Get size of velocity distribution array - open(unit=LUN,file=user%velfile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%velfile)),status="old",iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ",trim(user%velfile) stop @@ -65,15 +65,15 @@ subroutine init_dist(user,domain) domain%vnum = domain%vnum + 1 end do if (domain%vnum == 0) then - write(*,*) "No valid entries in ",trim(user%velfile) + write(*,*) "No valid entries in ",trim(adjustl(user%velfile)) end if close(LUN) ! Get size of real crater list array if (user%doquasimc) then - open(unit=LUN,file=rcfile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(rcfile)),status="old",iostat=ierr) if (ierr /= 0) then - write(*,*) "Unable to open file ", trim(rcfile) + write(*,*) "Unable to open file ", trim(adjustl(rcfile)) stop end if diff --git a/src/init/init_domain.f90 b/src/init/init_domain.f90 index a1010f28..79317177 100644 --- a/src/init/init_domain.f90 +++ b/src/init/init_domain.f90 @@ -75,8 +75,8 @@ subroutine init_domain(user,crater,domain,prod,pdist,vdist,crtscl,nflux) end select ! Preliminary seismic property calculations - user%seisk = THIRD * user%tvel * user%tfrac - user%cohaccel = user%regcoh / user%trho_r + ! user%seisk = THIRD * user%tvel * user%tfrac + ! user%cohaccel = user%regcoh / user%trho_r ! Now we build an idealized production population domain%initialize = .true. diff --git a/src/init/init_regolith_stack.f90 b/src/init/init_regolith_stack.f90 index 7115a63c..3538909f 100644 --- a/src/init/init_regolith_stack.f90 +++ b/src/init/init_regolith_stack.f90 @@ -37,9 +37,10 @@ subroutine init_regolith_stack(user,surf,domain) !======================================= ! Initialize the grid space !======================================= - - do concurrent(xp=1:user%gridsize,yp=1:user%gridsize) - call util_init_array(user,surf(xp,yp)%regolayer,domain,initstat) + do yp=1,user%gridsize + do xp=1,user%gridsize + call util_init_array(user,surf(xp,yp)%regolayer,domain,initstat) + end do end do return diff --git a/src/io/io_ejecta_table.f90 b/src/io/io_ejecta_table.f90 index f5e1287a..f2f178a0 100644 --- a/src/io/io_ejecta_table.f90 +++ b/src/io/io_ejecta_table.f90 @@ -35,12 +35,12 @@ subroutine io_ejecta_table(crater,domain,ejb,ejtble,filename) integer(I4B) :: k ! Executable code - open(LUN, FILE=filename, status='replace') + open(LUN, FILE=trim(adjustl(filename)), status='replace') write(LUN,'("# trad = ",ES12.5, " frad = ",ES12.5)') crater%rad,crater%frad write(LUN,'("# ejrim = ",ES12.5, " ejdis = ",ES12.5," imp = ",ES12.5)') crater%ejrim,crater%ejdis,crater%imp write(LUN,'(A63)') '# "r (m)" "h (m)" "v (m/s)" "ang (deg)" "erad (m)"' do k=1,ejtble - write(LUN,'(5(ES13.5E3,1X))') exp(ejb(k)%lrad),exp(ejb(k)%thick),sqrt(ejb(k)%vesq),ejb(k)%angle/DEG2RAD, & + write(LUN,'(5(ES13.5E3,1X))') exp(ejb(k)%lrad),exp(ejb(k)%thick),sqrt(exp(ejb(k)%vesq)),exp(ejb(k)%angle)/DEG2RAD, & exp(ejb(k)%erad) end do close(LUN) diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index b89b2612..4a528246 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -97,7 +97,7 @@ subroutine io_input(infile,user) user%dotopodiffusion = .true. write(user%sfdfile,*) trim(adjustl(SFDFILE)) - open(unit=LUN,file=infile,status="old",iostat=ierr) + open(unit=LUN,file=trim(adjustl(infile)),status="old",iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ",trim(infile) stop diff --git a/src/io/io_read_craterlist.f90 b/src/io/io_read_craterlist.f90 index ec17ff58..35b4e38d 100644 --- a/src/io/io_read_craterlist.f90 +++ b/src/io/io_read_craterlist.f90 @@ -35,16 +35,16 @@ subroutine io_read_craterlist(rclist, user, domain) ! Read the next crater from the craterlist - open(unit=LUN,file=rcfile,status='old',iostat=ierr) + open(unit=LUN,file=trim(adjustl(rcfile)),status='old',iostat=ierr) if (ierr /= 0) then - write(*,*) "Unable to open file ",trim(rcfile) + write(*,*) "Unable to open file ",trim(adjustl(rcfile)) stop end if do i=1,domain%rcnum read(LUN,*,iostat=ierr) rclist(1:6,i) if (ierr/=0) then - write(*,*) "Unable to read file ",trim(rcfile) + write(*,*) "Unable to read file ",trim(adjustl(rcfile)) stop end if end do diff --git a/src/io/io_read_prod.f90 b/src/io/io_read_prod.f90 index 2eb47516..8a845980 100644 --- a/src/io/io_read_prod.f90 +++ b/src/io/io_read_prod.f90 @@ -35,7 +35,7 @@ subroutine io_read_prod(prod,user,domain) ! Executable code ! Read in the size-frequency distribution file - open(unit=LUN,file=user%sfdfile,status='old',iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%sfdfile)),status='old',iostat=ierr) if (ierr /= 0) then write(*,*) "Unable to open file ",trim(user%sfdfile) stop diff --git a/src/io/io_read_vdist.f90 b/src/io/io_read_vdist.f90 index a15dc8dd..84309feb 100644 --- a/src/io/io_read_vdist.f90 +++ b/src/io/io_read_vdist.f90 @@ -35,16 +35,16 @@ subroutine io_read_vdist(vdist,user,domain) ! Executable code ! Read in velocity distribution file - open(unit=LUN,file=user%velfile,status='old',iostat=ierr) + open(unit=LUN,file=trim(adjustl(user%velfile)),status='old',iostat=ierr) if (ierr /= 0) then - write(*,*) "Unable to open file ",trim(user%velfile) + write(*,*) "Unable to open file ",trim(adjustl(user%velfile)) stop end if do i=1,domain%vnum read(LUN,*,iostat=ierr) vdist(1:3,i) if (ierr/=0) then - write(*,*) "Unable to read file ",trim(user%velfile) + write(*,*) "Unable to read file ",trim(adjustl(user%velfile)) stop end if end do diff --git a/src/regolith/regolith_melt_fraction.f90 b/src/regolith/regolith_melt_fraction.f90 index 4886d9fc..01dac4ac 100755 --- a/src/regolith/regolith_melt_fraction.f90 +++ b/src/regolith/regolith_melt_fraction.f90 @@ -102,5 +102,11 @@ subroutine regolith_melt_fraction(dimp,depthb,erad1,erad2,rmelt,meltfrac) write(*,*) 'regolith_melt_fraction: this is a bug!' end if + if (meltfrac < epsilon(1.0_DP)) then + meltfrac = 0.0_DP + else if (meltfrac > 1.0_DP) then + meltfrac = 1.0_DP + end if + return end subroutine regolith_melt_fraction diff --git a/src/regolith/regolith_streamtube.f90 b/src/regolith/regolith_streamtube.f90 index 28c28d08..e6d04702 100644 --- a/src/regolith/regolith_streamtube.f90 +++ b/src/regolith/regolith_streamtube.f90 @@ -53,7 +53,7 @@ ! Arguments : surf :: surface ! ! -! Notes : +! Notes : 'eradc' is the center of the ejection radius. 'eradi' is inner, 'erado' is outer. 'cnt' is counting number. 'cmax', 'ri', and 'rip1' concern streamtube geometry. ! !********************************************************************************************************************************** subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi,lrad,ebh,rm,vsq,volm) @@ -154,7 +154,7 @@ subroutine regolith_streamtube(user,surf,crater,domain,ejb,ejtble,xp,yp,xpi,ypi, end if if (eradc<=0.0_DP) then - write(*,*) k,ebh,crater%ejdis,lrad,exp(ejb(k-1)%lrad),exp(ejb(k)%lrad),eradc,ejb(k-1)%erad,ejb(k)%erad + write(*,*) k,ebh,crater%ejdis,lrad,exp(ejb(k-1)%lrad),exp(ejb(k)%lrad),eradc,exp(ejb(k-1)%erad),exp(ejb(k)%erad) stop end if diff --git a/src/regolith/regolith_transport.f90 b/src/regolith/regolith_transport.f90 index 01ca20a1..ca832102 100644 --- a/src/regolith/regolith_transport.f90 +++ b/src/regolith/regolith_transport.f90 @@ -63,6 +63,11 @@ subroutine regolith_transport(user,surfi,crater,domain,ejb,ejtble,lrad,ebh,newla frac = (loglrad - logtablerad) / logdelta melt = ejb(k)%meltfrac - ((ejb(k)%meltfrac - ejb(k+1)%meltfrac) * frac) end if + if (melt < LOGVSMALL) then + melt = 0.0_DP + else + melt = exp(melt) + end if newlayer%meltvolume = melt * newlayer%totvolume diff --git a/src/util/util_sort_layer.f90 b/src/util/util_sort_layer.f90 index c3d055b1..65cd70d0 100644 --- a/src/util/util_sort_layer.f90 +++ b/src/util/util_sort_layer.f90 @@ -57,7 +57,7 @@ subroutine util_sort_layer(user,surf,crater) temptime = surf(mx,my)%timestamp(1:user%numlayers) ! Sort the layers by crater diameter - call util_mrgrnk(surf(mx,my)%diam(1:user%numlayers),ind) + call util_mrgrnk(tempdiam,ind) do k=1,user%numlayers surf(mx,my)%diam(k)=tempdiam(ind(k)) diff --git a/version.txt b/version.txt index 3bae6081..032892e0 100644 --- a/version.txt +++ b/version.txt @@ -1 +1 @@ -2023.10.0 \ No newline at end of file +2024.2.0 \ No newline at end of file