From b9e63d8f3911884cd408cc16a0783e1dc2cd30e6 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Tue, 8 Mar 2022 15:45:49 -0500 Subject: [PATCH 01/25] quasi-MC now works with the updated python!! However, the time conversion may be messed up --- compile.bash | 4 ++ .../.idea/global-lunar-bombardment.iml | 2 +- .../global-lunar-bombardment/.idea/misc.xml | 4 -- examples/global-lunar-bombardment/ctem.in | 4 +- .../{driver.py => ctemdriver.py} | 0 python/ctem/ctem/driver.py | 29 ++++++++++-- python/ctem/ctem/util.py | 44 ++++++++++++++++++- 7 files changed, 74 insertions(+), 13 deletions(-) create mode 100644 compile.bash delete mode 100644 examples/global-lunar-bombardment/.idea/misc.xml mode change 100755 => 100644 examples/global-lunar-bombardment/ctem.in rename examples/global-lunar-bombardment/{driver.py => ctemdriver.py} (100%) diff --git a/compile.bash b/compile.bash new file mode 100644 index 00000000..8bf27cde --- /dev/null +++ b/compile.bash @@ -0,0 +1,4 @@ +./autogen.sh +cd build +../configure +make clean; make \ No newline at end of file diff --git a/examples/global-lunar-bombardment/.idea/global-lunar-bombardment.iml b/examples/global-lunar-bombardment/.idea/global-lunar-bombardment.iml index d0876a78..701b9669 100644 --- a/examples/global-lunar-bombardment/.idea/global-lunar-bombardment.iml +++ b/examples/global-lunar-bombardment/.idea/global-lunar-bombardment.iml @@ -2,7 +2,7 @@ - + \ No newline at end of file diff --git a/examples/global-lunar-bombardment/.idea/misc.xml b/examples/global-lunar-bombardment/.idea/misc.xml deleted file mode 100644 index a2e120dc..00000000 --- a/examples/global-lunar-bombardment/.idea/misc.xml +++ /dev/null @@ -1,4 +0,0 @@ - - - - \ No newline at end of file diff --git a/examples/global-lunar-bombardment/ctem.in b/examples/global-lunar-bombardment/ctem.in old mode 100755 new mode 100644 index ca50a7fa..470aaee4 --- a/examples/global-lunar-bombardment/ctem.in +++ b/examples/global-lunar-bombardment/ctem.in @@ -34,9 +34,9 @@ runtype single ! Run type: options are normal / ! CTEM required inputs seed 76535 ! Random number generator seed -gridsize 200 ! Size of grid in pixels +gridsize 2000 ! Size of grid in pixels numlayers 10 ! Number of perched layers -pix 3.08e4 ! Pixel size (m) +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 diff --git a/examples/global-lunar-bombardment/driver.py b/examples/global-lunar-bombardment/ctemdriver.py similarity index 100% rename from examples/global-lunar-bombardment/driver.py rename to examples/global-lunar-bombardment/ctemdriver.py diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 8f565f1e..d843de5a 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -46,8 +46,9 @@ def __init__(self, param_file="ctem.in", isnew=True): 'sfdcompare': None, 'sfdfile': None, 'quasimc': None, - 'realcraterlist': None + 'realcraterlist': None, } + self.user = util.read_user_input(self.user) # This is containsall files generated by the main Fortran CTEM program plus the ctem.dat file that @@ -113,8 +114,28 @@ def __init__(self, param_file="ctem.in", isnew=True): #Read list of real craters print("quasi-MC mode is ON") - # Use self.compute_one_interval() to generate craterlist.dat rclist = util.read_formatted_ascii(self.user['realcraterlist'], skip_lines = 0) + self.tempfile = os.path.join(currentdir, 'temp.in') + + # Generate craterlist.dat + with open(self.user['ctemfile']) as perm, open(self.tempfile, "w") as temp: + for line in perm: + temp.write(line) + + #Write a temporary input file to generate the necessary quasimc files + self.permfile = self.user['ctemfile'] + self.temp = util.write_temp_input(self.tempfile) + self.user['ctemfile'] = self.tempfile + self.user = util.read_user_input(self.user) + self.permin = os.path.join(currentdir, 'perm.in') + self.ctemin = os.path.join(currentdir, 'ctem.in') + os.rename(self.permfile, self.permin) + os.rename(self.tempfile, self.ctemin) + self.run() + + os.replace(self.permin, self.ctemin) + self.user['ctemfile'] = self.permfile + self.user = util.read_user_input(self.user) #Interpolate craterscale.dat to get impactor sizes from crater sizes given df = pandas.read_csv(self.output_filenames['craterscale'], sep='\s+') @@ -130,7 +151,7 @@ def __init__(self, param_file="ctem.in", isnew=True): rclist = rclist[rclist[:,5].argsort()] #Export to dat file - util.write_realcraters(user, rclist) + util.write_realcraters(self.output_filenames['craterlist'], rclist) util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) else: @@ -179,7 +200,7 @@ def run(self): # Read in output files self.read_output() - # Process the o utput files + # Process the output files self.process_output() # Update ncount diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index 8cab22e3..773235eb 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -4,6 +4,8 @@ from matplotlib.colors import LightSource import matplotlib.cm as cm import matplotlib.pyplot as plt +import re +from tempfile import mkstemp # Set pixel scaling common for image writing, at 1 pixel/ array element dpi = 72.0 @@ -301,6 +303,7 @@ def read_user_input(user): if ('impfile' == fields[0].lower()): user['impfile'] = os.path.join(user['workingdir'],fields[1]) if ('maxcrat' == fields[0].lower()): user['maxcrat'] = real2float(fields[1]) if ('sfdcompare' == fields[0].lower()): user['sfdcompare'] = os.path.join(user['workingdir'], fields[1]) + if ('realcraterlist' == fields[0].lower()): user['realcraterlist'] = os.path.join(user['workingdir'], fields[1]) if ('interval' == fields[0].lower()): user['interval'] = real2float(fields[1]) if ('numintervals' == fields[0].lower()): user['numintervals'] = int(fields[1]) if ('popupconsole' == fields[0].lower()): user['popupconsole'] = fields[1] @@ -378,6 +381,33 @@ def real2float(realstr): """ return float(realstr.replace('d', 'E').replace('D', 'E')) +def sed(pattern, replace, source, count=0): + """Python implementation of unix sed command; not fully functional sed.""" + + fin = open(source, 'r') + num_replaced = 0 + + fd, name = mkstemp() + fout = open(name, 'w') + + for line in fin: + out = re.sub(pattern, replace, line) + fout.write(out) + + if out != line: + num_replaced += 1 + if count and num_replaced > count: + break + + fout.writelines(fin.readlines()) + + + fin.close() + fout.close() + + shutil.move(name, source) + + return def write_datfile(user, filename, seedarr): # Write various user and random number seeds into ctem.dat file @@ -402,10 +432,20 @@ def write_production(user, production): return -def write_realcraters(user, realcraters): +def write_realcraters(filename, realcraters): """Writes file of real craters for use in quasi-MC runs""" - filename = user['craterlist'] np.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t') + return + +def write_temp_input(filename): + """Makes changes to a temporary input file for use when generating craterlist.dat for quasimc runs""" + + sed('testflag', 'testflag T!', filename) + sed('testimp', 'testimp 10 !', filename) + sed('quasimc', 'quasimc F!', filename) + sed('interval', 'interval 1 !', filename) + sed('numinterval 1 !s', 'numintervals 1 !', filename) + return \ No newline at end of file From 12a318e5a0a81155f1f12f89381a5844e38c2226 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 9 Mar 2022 13:24:00 -0500 Subject: [PATCH 02/25] Added ability to read in the name of the input file as a command line argument. Default is ctem.in --- src/main/CTEM.f90 | 10 ++++++++-- 1 file changed, 8 insertions(+), 2 deletions(-) diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index 986356c8..db31e955 100644 --- a/src/main/CTEM.f90 +++ b/src/main/CTEM.f90 @@ -63,11 +63,17 @@ program CTEM real(DP) :: lambda !$ real(DP) :: t1,t2 real(DP),dimension(:,:),allocatable :: nflux +integer(I4B) :: narg, ierr !$ t1 = omp_get_wtime() call io_splash() -!write(*,*) 'Reading input files' -infile="ctem.in" +narg = command_argument_count() +if (narg /= 0) then + call get_command_argument(1, infile, status = ierr) +else + infile="ctem.in" +end if +write(*,*) 'Reading input file ',trim(adjustl(infile)) call io_input(infile,user) ! Initialize distribution arrays (crater size, number) From e15622e7a2f2bb12bd810f68d686c948c133246b Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Mar 2022 13:27:29 -0500 Subject: [PATCH 03/25] Added ability of ctem driver to supply alternative input files --- python/ctem/ctem/driver.py | 7 ++++--- 1 file changed, 4 insertions(+), 3 deletions(-) diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 3c198a56..08999077 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -231,17 +231,18 @@ def run(self): return - def compute_one_interval(self): + def compute_one_interval(self, ctemin="ctem.in"): """ Executes the Fortran code to generate one interval of simulation output. """ # Create crater population and display CTEM progress on screen print(self.user['ncount'], ' Calling FORTRAN routine') try: - p = subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')], + p = subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM'), ctemin], stdout=subprocess.PIPE, stderr=subprocess.PIPE, - universal_newlines=True) + universal_newlines=True, + shell=False) for line in p.stdout: print(line, end='') res = p.communicate() From 36b0bf20f8b94a86ce2e12c85fe429259d12a1f3 Mon Sep 17 00:00:00 2001 From: Austin Michael Blevins Date: Wed, 9 Mar 2022 14:54:05 -0500 Subject: [PATCH 04/25] quasimc works with the correct time intervals --- INSTALL | 322 +++++++++++++++++++------------------ python/ctem/ctem/driver.py | 3 +- src/Makefile.am | 2 +- 3 files changed, 165 insertions(+), 162 deletions(-) diff --git a/INSTALL b/INSTALL index 8865734f..007e9396 100644 --- a/INSTALL +++ b/INSTALL @@ -1,8 +1,8 @@ Installation Instructions ************************* - Copyright (C) 1994-1996, 1999-2002, 2004-2016 Free Software -Foundation, Inc. +Copyright (C) 1994-1996, 1999-2002, 2004-2013 Free Software Foundation, +Inc. Copying and distribution of this file, with or without modification, are permitted in any medium without royalty provided the copyright @@ -12,96 +12,97 @@ without warranty of any kind. Basic Installation ================== - Briefly, the shell command './configure && make && make install' -should configure, build, and install this package. The following -more-detailed instructions are generic; see the 'README' file for + Briefly, the shell commands `./configure; make; make install' should +configure, build, and install this package. The following +more-detailed instructions are generic; see the `README' file for instructions specific to this package. Some packages provide this -'INSTALL' file but do not implement all of the features documented +`INSTALL' file but do not implement all of the features documented below. The lack of an optional feature in a given package is not necessarily a bug. More recommendations for GNU packages can be found in *note Makefile Conventions: (standards)Makefile Conventions. - The 'configure' shell script attempts to guess correct values for + The `configure' shell script attempts to guess correct values for various system-dependent variables used during compilation. It uses -those values to create a 'Makefile' in each directory of the package. -It may also create one or more '.h' files containing system-dependent -definitions. Finally, it creates a shell script 'config.status' that +those values to create a `Makefile' in each directory of the package. +It may also create one or more `.h' files containing system-dependent +definitions. Finally, it creates a shell script `config.status' that you can run in the future to recreate the current configuration, and a -file 'config.log' containing compiler output (useful mainly for -debugging 'configure'). +file `config.log' containing compiler output (useful mainly for +debugging `configure'). - It can also use an optional file (typically called 'config.cache' and -enabled with '--cache-file=config.cache' or simply '-C') that saves the -results of its tests to speed up reconfiguring. Caching is disabled by -default to prevent problems with accidental use of stale cache files. + It can also use an optional file (typically called `config.cache' +and enabled with `--cache-file=config.cache' or simply `-C') that saves +the results of its tests to speed up reconfiguring. Caching is +disabled by default to prevent problems with accidental use of stale +cache files. If you need to do unusual things to compile the package, please try -to figure out how 'configure' could check whether to do them, and mail -diffs or instructions to the address given in the 'README' so they can +to figure out how `configure' could check whether to do them, and mail +diffs or instructions to the address given in the `README' so they can be considered for the next release. If you are using the cache, and at -some point 'config.cache' contains results you don't want to keep, you +some point `config.cache' contains results you don't want to keep, you may remove or edit it. - The file 'configure.ac' (or 'configure.in') is used to create -'configure' by a program called 'autoconf'. You need 'configure.ac' if -you want to change it or regenerate 'configure' using a newer version of -'autoconf'. + The file `configure.ac' (or `configure.in') is used to create +`configure' by a program called `autoconf'. You need `configure.ac' if +you want to change it or regenerate `configure' using a newer version +of `autoconf'. The simplest way to compile this package is: - 1. 'cd' to the directory containing the package's source code and type - './configure' to configure the package for your system. + 1. `cd' to the directory containing the package's source code and type + `./configure' to configure the package for your system. - Running 'configure' might take a while. While running, it prints + Running `configure' might take a while. While running, it prints some messages telling which features it is checking for. - 2. Type 'make' to compile the package. + 2. Type `make' to compile the package. - 3. Optionally, type 'make check' to run any self-tests that come with + 3. Optionally, type `make check' to run any self-tests that come with the package, generally using the just-built uninstalled binaries. - 4. Type 'make install' to install the programs and any data files and + 4. Type `make install' to install the programs and any data files and documentation. When installing into a prefix owned by root, it is recommended that the package be configured and built as a regular - user, and only the 'make install' phase executed with root + user, and only the `make install' phase executed with root privileges. - 5. Optionally, type 'make installcheck' to repeat any self-tests, but + 5. Optionally, type `make installcheck' to repeat any self-tests, but this time using the binaries in their final installed location. This target does not install anything. Running this target as a - regular user, particularly if the prior 'make install' required + regular user, particularly if the prior `make install' required root privileges, verifies that the installation completed correctly. 6. You can remove the program binaries and object files from the - source code directory by typing 'make clean'. To also remove the - files that 'configure' created (so you can compile the package for - a different kind of computer), type 'make distclean'. There is - also a 'make maintainer-clean' target, but that is intended mainly + source code directory by typing `make clean'. To also remove the + files that `configure' created (so you can compile the package for + a different kind of computer), type `make distclean'. There is + also a `make maintainer-clean' target, but that is intended mainly for the package's developers. If you use it, you may have to get all sorts of other programs in order to regenerate files that came with the distribution. - 7. Often, you can also type 'make uninstall' to remove the installed + 7. Often, you can also type `make uninstall' to remove the installed files again. In practice, not all packages have tested that uninstallation works correctly, even though it is required by the GNU Coding Standards. - 8. Some packages, particularly those that use Automake, provide 'make + 8. Some packages, particularly those that use Automake, provide `make distcheck', which can by used by developers to test that all other - targets like 'make install' and 'make uninstall' work correctly. + targets like `make install' and `make uninstall' work correctly. This target is generally not run by end users. Compilers and Options ===================== Some systems require unusual options for compilation or linking that -the 'configure' script does not know about. Run './configure --help' +the `configure' script does not know about. Run `./configure --help' for details on some of the pertinent environment variables. - You can give 'configure' initial values for configuration parameters -by setting variables in the command line or in the environment. Here is -an example: + You can give `configure' initial values for configuration parameters +by setting variables in the command line or in the environment. Here +is an example: ./configure CC=c99 CFLAGS=-g LIBS=-lposix @@ -112,21 +113,21 @@ Compiling For Multiple Architectures You can compile the package for more than one kind of computer at the same time, by placing the object files for each architecture in their -own directory. To do this, you can use GNU 'make'. 'cd' to the +own directory. To do this, you can use GNU `make'. `cd' to the directory where you want the object files and executables to go and run -the 'configure' script. 'configure' automatically checks for the source -code in the directory that 'configure' is in and in '..'. This is known -as a "VPATH" build. +the `configure' script. `configure' automatically checks for the +source code in the directory that `configure' is in and in `..'. This +is known as a "VPATH" build. - With a non-GNU 'make', it is safer to compile the package for one + With a non-GNU `make', it is safer to compile the package for one architecture at a time in the source code directory. After you have -installed the package for one architecture, use 'make distclean' before +installed the package for one architecture, use `make distclean' before reconfiguring for another architecture. On MacOS X 10.5 and later systems, you can create libraries and executables that work on multiple system types--known as "fat" or -"universal" binaries--by specifying multiple '-arch' options to the -compiler but only a single '-arch' option to the preprocessor. Like +"universal" binaries--by specifying multiple `-arch' options to the +compiler but only a single `-arch' option to the preprocessor. Like this: ./configure CC="gcc -arch i386 -arch x86_64 -arch ppc -arch ppc64" \ @@ -135,104 +136,105 @@ this: This is not guaranteed to produce working output in all cases, you may have to build one architecture at a time and combine the results -using the 'lipo' tool if you have problems. +using the `lipo' tool if you have problems. Installation Names ================== - By default, 'make install' installs the package's commands under -'/usr/local/bin', include files under '/usr/local/include', etc. You -can specify an installation prefix other than '/usr/local' by giving -'configure' the option '--prefix=PREFIX', where PREFIX must be an + By default, `make install' installs the package's commands under +`/usr/local/bin', include files under `/usr/local/include', etc. You +can specify an installation prefix other than `/usr/local' by giving +`configure' the option `--prefix=PREFIX', where PREFIX must be an absolute file name. You can specify separate installation prefixes for architecture-specific files and architecture-independent files. If you -pass the option '--exec-prefix=PREFIX' to 'configure', the package uses +pass the option `--exec-prefix=PREFIX' to `configure', the package uses PREFIX as the prefix for installing programs and libraries. Documentation and other data files still use the regular prefix. In addition, if you use an unusual directory layout you can give -options like '--bindir=DIR' to specify different values for particular -kinds of files. Run 'configure --help' for a list of the directories -you can set and what kinds of files go in them. In general, the default -for these options is expressed in terms of '${prefix}', so that -specifying just '--prefix' will affect all of the other directory +options like `--bindir=DIR' to specify different values for particular +kinds of files. Run `configure --help' for a list of the directories +you can set and what kinds of files go in them. In general, the +default for these options is expressed in terms of `${prefix}', so that +specifying just `--prefix' will affect all of the other directory specifications that were not explicitly provided. The most portable way to affect installation locations is to pass the -correct locations to 'configure'; however, many packages provide one or +correct locations to `configure'; however, many packages provide one or both of the following shortcuts of passing variable assignments to the -'make install' command line to change installation locations without +`make install' command line to change installation locations without having to reconfigure or recompile. The first method involves providing an override variable for each -affected directory. For example, 'make install +affected directory. For example, `make install prefix=/alternate/directory' will choose an alternate location for all directory configuration variables that were expressed in terms of -'${prefix}'. Any directories that were specified during 'configure', -but not in terms of '${prefix}', must each be overridden at install time -for the entire installation to be relocated. The approach of makefile -variable overrides for each directory variable is required by the GNU -Coding Standards, and ideally causes no recompilation. However, some -platforms have known limitations with the semantics of shared libraries -that end up requiring recompilation when using this method, particularly -noticeable in packages that use GNU Libtool. - - The second method involves providing the 'DESTDIR' variable. For -example, 'make install DESTDIR=/alternate/directory' will prepend -'/alternate/directory' before all installation names. The approach of -'DESTDIR' overrides is not required by the GNU Coding Standards, and +`${prefix}'. Any directories that were specified during `configure', +but not in terms of `${prefix}', must each be overridden at install +time for the entire installation to be relocated. The approach of +makefile variable overrides for each directory variable is required by +the GNU Coding Standards, and ideally causes no recompilation. +However, some platforms have known limitations with the semantics of +shared libraries that end up requiring recompilation when using this +method, particularly noticeable in packages that use GNU Libtool. + + The second method involves providing the `DESTDIR' variable. For +example, `make install DESTDIR=/alternate/directory' will prepend +`/alternate/directory' before all installation names. The approach of +`DESTDIR' overrides is not required by the GNU Coding Standards, and does not work on platforms that have drive letters. On the other hand, it does better at avoiding recompilation issues, and works well even -when some directory options were not specified in terms of '${prefix}' -at 'configure' time. +when some directory options were not specified in terms of `${prefix}' +at `configure' time. Optional Features ================= If the package supports it, you can cause programs to be installed -with an extra prefix or suffix on their names by giving 'configure' the -option '--program-prefix=PREFIX' or '--program-suffix=SUFFIX'. - - Some packages pay attention to '--enable-FEATURE' options to -'configure', where FEATURE indicates an optional part of the package. -They may also pay attention to '--with-PACKAGE' options, where PACKAGE -is something like 'gnu-as' or 'x' (for the X Window System). The -'README' should mention any '--enable-' and '--with-' options that the +with an extra prefix or suffix on their names by giving `configure' the +option `--program-prefix=PREFIX' or `--program-suffix=SUFFIX'. + + Some packages pay attention to `--enable-FEATURE' options to +`configure', where FEATURE indicates an optional part of the package. +They may also pay attention to `--with-PACKAGE' options, where PACKAGE +is something like `gnu-as' or `x' (for the X Window System). The +`README' should mention any `--enable-' and `--with-' options that the package recognizes. - For packages that use the X Window System, 'configure' can usually + For packages that use the X Window System, `configure' can usually find the X include and library files automatically, but if it doesn't, -you can use the 'configure' options '--x-includes=DIR' and -'--x-libraries=DIR' to specify their locations. +you can use the `configure' options `--x-includes=DIR' and +`--x-libraries=DIR' to specify their locations. Some packages offer the ability to configure how verbose the -execution of 'make' will be. For these packages, running './configure +execution of `make' will be. For these packages, running `./configure --enable-silent-rules' sets the default to minimal output, which can be -overridden with 'make V=1'; while running './configure +overridden with `make V=1'; while running `./configure --disable-silent-rules' sets the default to verbose, which can be -overridden with 'make V=0'. +overridden with `make V=0'. Particular systems ================== - On HP-UX, the default C compiler is not ANSI C compatible. If GNU CC -is not installed, it is recommended to use the following options in + On HP-UX, the default C compiler is not ANSI C compatible. If GNU +CC is not installed, it is recommended to use the following options in order to use an ANSI C compiler: ./configure CC="cc -Ae -D_XOPEN_SOURCE=500" and if that doesn't work, install pre-built binaries of GCC for HP-UX. - HP-UX 'make' updates targets which have the same time stamps as their -prerequisites, which makes it generally unusable when shipped generated -files such as 'configure' are involved. Use GNU 'make' instead. + HP-UX `make' updates targets which have the same time stamps as +their prerequisites, which makes it generally unusable when shipped +generated files such as `configure' are involved. Use GNU `make' +instead. On OSF/1 a.k.a. Tru64, some versions of the default C compiler cannot -parse its '' header file. The option '-nodtk' can be used as a -workaround. If GNU CC is not installed, it is therefore recommended to -try +parse its `' header file. The option `-nodtk' can be used as +a workaround. If GNU CC is not installed, it is therefore recommended +to try ./configure CC="cc" @@ -240,26 +242,26 @@ and if that doesn't work, try ./configure CC="cc -nodtk" - On Solaris, don't put '/usr/ucb' early in your 'PATH'. This + On Solaris, don't put `/usr/ucb' early in your `PATH'. This directory contains several dysfunctional programs; working variants of -these programs are available in '/usr/bin'. So, if you need '/usr/ucb' -in your 'PATH', put it _after_ '/usr/bin'. +these programs are available in `/usr/bin'. So, if you need `/usr/ucb' +in your `PATH', put it _after_ `/usr/bin'. - On Haiku, software installed for all users goes in '/boot/common', -not '/usr/local'. It is recommended to use the following options: + On Haiku, software installed for all users goes in `/boot/common', +not `/usr/local'. It is recommended to use the following options: ./configure --prefix=/boot/common Specifying the System Type ========================== - There may be some features 'configure' cannot figure out + There may be some features `configure' cannot figure out automatically, but needs to determine by the type of machine the package will run on. Usually, assuming the package is built to be run on the -_same_ architectures, 'configure' can figure that out, but if it prints +_same_ architectures, `configure' can figure that out, but if it prints a message saying it cannot guess the machine type, give it the -'--build=TYPE' option. TYPE can either be a short name for the system -type, such as 'sun4', or a canonical name which has the form: +`--build=TYPE' option. TYPE can either be a short name for the system +type, such as `sun4', or a canonical name which has the form: CPU-COMPANY-SYSTEM @@ -268,101 +270,101 @@ where SYSTEM can have one of these forms: OS KERNEL-OS - See the file 'config.sub' for the possible values of each field. If -'config.sub' isn't included in this package, then this package doesn't + See the file `config.sub' for the possible values of each field. If +`config.sub' isn't included in this package, then this package doesn't need to know the machine type. If you are _building_ compiler tools for cross-compiling, you should -use the option '--target=TYPE' to select the type of system they will +use the option `--target=TYPE' to select the type of system they will produce code for. If you want to _use_ a cross compiler, that generates code for a platform different from the build platform, you should specify the "host" platform (i.e., that on which the generated programs will -eventually be run) with '--host=TYPE'. +eventually be run) with `--host=TYPE'. Sharing Defaults ================ - If you want to set default values for 'configure' scripts to share, -you can create a site shell script called 'config.site' that gives -default values for variables like 'CC', 'cache_file', and 'prefix'. -'configure' looks for 'PREFIX/share/config.site' if it exists, then -'PREFIX/etc/config.site' if it exists. Or, you can set the -'CONFIG_SITE' environment variable to the location of the site script. -A warning: not all 'configure' scripts look for a site script. + If you want to set default values for `configure' scripts to share, +you can create a site shell script called `config.site' that gives +default values for variables like `CC', `cache_file', and `prefix'. +`configure' looks for `PREFIX/share/config.site' if it exists, then +`PREFIX/etc/config.site' if it exists. Or, you can set the +`CONFIG_SITE' environment variable to the location of the site script. +A warning: not all `configure' scripts look for a site script. Defining Variables ================== Variables not defined in a site shell script can be set in the -environment passed to 'configure'. However, some packages may run +environment passed to `configure'. However, some packages may run configure again during the build, and the customized values of these variables may be lost. In order to avoid this problem, you should set -them in the 'configure' command line, using 'VAR=value'. For example: +them in the `configure' command line, using `VAR=value'. For example: ./configure CC=/usr/local2/bin/gcc -causes the specified 'gcc' to be used as the C compiler (unless it is +causes the specified `gcc' to be used as the C compiler (unless it is overridden in the site shell script). -Unfortunately, this technique does not work for 'CONFIG_SHELL' due to an -Autoconf limitation. Until the limitation is lifted, you can use this -workaround: +Unfortunately, this technique does not work for `CONFIG_SHELL' due to +an Autoconf limitation. Until the limitation is lifted, you can use +this workaround: CONFIG_SHELL=/bin/bash ./configure CONFIG_SHELL=/bin/bash -'configure' Invocation +`configure' Invocation ====================== - 'configure' recognizes the following options to control how it + `configure' recognizes the following options to control how it operates. -'--help' -'-h' - Print a summary of all of the options to 'configure', and exit. +`--help' +`-h' + Print a summary of all of the options to `configure', and exit. -'--help=short' -'--help=recursive' +`--help=short' +`--help=recursive' Print a summary of the options unique to this package's - 'configure', and exit. The 'short' variant lists options used only - in the top level, while the 'recursive' variant lists options also - present in any nested packages. + `configure', and exit. The `short' variant lists options used + only in the top level, while the `recursive' variant lists options + also present in any nested packages. -'--version' -'-V' - Print the version of Autoconf used to generate the 'configure' +`--version' +`-V' + Print the version of Autoconf used to generate the `configure' script, and exit. -'--cache-file=FILE' +`--cache-file=FILE' Enable the cache: use and save the results of the tests in FILE, - traditionally 'config.cache'. FILE defaults to '/dev/null' to + traditionally `config.cache'. FILE defaults to `/dev/null' to disable caching. -'--config-cache' -'-C' - Alias for '--cache-file=config.cache'. +`--config-cache' +`-C' + Alias for `--cache-file=config.cache'. -'--quiet' -'--silent' -'-q' +`--quiet' +`--silent' +`-q' Do not print messages saying which checks are being made. To - suppress all normal output, redirect it to '/dev/null' (any error + suppress all normal output, redirect it to `/dev/null' (any error messages will still be shown). -'--srcdir=DIR' +`--srcdir=DIR' Look for the package's source code in directory DIR. Usually - 'configure' can determine that directory automatically. + `configure' can determine that directory automatically. -'--prefix=DIR' - Use DIR as the installation prefix. *note Installation Names:: for - more details, including other options available for fine-tuning the - installation locations. +`--prefix=DIR' + Use DIR as the installation prefix. *note Installation Names:: + for more details, including other options available for fine-tuning + the installation locations. -'--no-create' -'-n' +`--no-create' +`-n' Run the configure checks, but stop before creating any output files. -'configure' also accepts some other, not widely useful, options. Run -'configure --help' for more details. +`configure' also accepts some other, not widely useful, options. Run +`configure --help' for more details. diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 08999077..7371770f 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -137,7 +137,8 @@ def __init__(self, param_file="ctem.in", isnew=True): self.ctemin = os.path.join(currentdir, 'ctem.in') os.rename(self.permfile, self.permin) os.rename(self.tempfile, self.ctemin) - self.run() + util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) + self.compute_one_interval() os.replace(self.permin, self.ctemin) self.user['ctemfile'] = self.permfile diff --git a/src/Makefile.am b/src/Makefile.am index 3a193de4..90ac9873 100755 --- a/src/Makefile.am +++ b/src/Makefile.am @@ -10,7 +10,7 @@ AM_FCFLAGS = $(IPRODUCTION) #ifort debug flags #gfortran optimized flags -AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace +#AM_FCFLAGS = -O3 -fopenmp -ffree-form -g -fbounds-check -fbacktrace #gfortran debug flags #AM_FCFLAGS = -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form -fsanitize-address-use-after-scope -fstack-check -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all From 46ece71b72ec8c3ed353f8349d4122c2cae5d259 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Mar 2022 17:15:39 -0500 Subject: [PATCH 05/25] Updated quasimc with new temp file system --- examples/quasimc-global/ctem.in | 19 +++++------ examples/quasimc-global/misc/.gitignore | 1 - examples/quasimc-global/surf/.gitignore | 3 -- python/ctem/ctem/driver.py | 45 +++++++++---------------- python/ctem/ctem/util.py | 5 +-- src/globals/module_globals.f90 | 7 ++-- src/io/io_input.f90 | 11 +++--- src/main/CTEM.f90 | 4 +-- 8 files changed, 36 insertions(+), 59 deletions(-) mode change 100755 => 100644 examples/quasimc-global/ctem.in delete mode 100644 examples/quasimc-global/misc/.gitignore delete mode 100644 examples/quasimc-global/surf/.gitignore diff --git a/examples/quasimc-global/ctem.in b/examples/quasimc-global/ctem.in old mode 100755 new mode 100644 index ca50a7fa..c5585032 --- a/examples/quasimc-global/ctem.in +++ b/examples/quasimc-global/ctem.in @@ -2,35 +2,35 @@ ! 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 15000.0 ! Diameter of test impactor (m) +testflag F ! 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.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 T ! MC run constrained by non-MC 'real' craters given in a list +quasimc T ! MC run constrained by non-MC 'real' craters given in a list realcraterlist craterlist.in ! list of 'real' craters for Quasi-MC runs ! IDL driver in uts -interval 6.280 -numintervals 100 ! Total number of intervals (total time = interval * numintervals) <--when runtype is 'single' +interval 6.280 +numintervals 100 ! 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 +popupconsole F ! Pop up console window every output interval 1 ! saveshaded F ! Output shaded relief images saverego F ! 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 (large file size) +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 intervals - ! statistical: surface is reset between intervals + ! 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 @@ -52,7 +52,6 @@ ybar_r 0.00e6 ! Regolith strength (Pa) 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) diff --git a/examples/quasimc-global/misc/.gitignore b/examples/quasimc-global/misc/.gitignore deleted file mode 100644 index 773a6df9..00000000 --- a/examples/quasimc-global/misc/.gitignore +++ /dev/null @@ -1 +0,0 @@ -*.dat diff --git a/examples/quasimc-global/surf/.gitignore b/examples/quasimc-global/surf/.gitignore deleted file mode 100644 index a17efa9a..00000000 --- a/examples/quasimc-global/surf/.gitignore +++ /dev/null @@ -1,3 +0,0 @@ -*.jpg -*.png - diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 7371770f..52174cc8 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -40,16 +40,16 @@ def __init__(self, param_file="ctem.in", isnew=True): 'shadedminh': 0.0, 'shadedmaxh': 0.0, 'workingdir': currentdir, - 'ctemfile': os.path.join(currentdir, param_file), + 'ctemfile': param_file, 'impfile': None, 'sfdcompare': None, - 'sfdfile': None, 'quasimc': None, + 'sfdfile' : None, 'realcraterlist': None, } self.user = util.read_user_input(self.user) - + # This is containsall files generated by the main Fortran CTEM program plus the ctem.dat file that # communicates the state of the simulation from the Python driver to the Fortran program self.output_filenames = { @@ -71,8 +71,11 @@ def __init__(self, param_file="ctem.in", isnew=True): 'ejmin' : 'ejecta_table_min.dat', 'testprof' : 'testprofile.dat', 'craterscale' : 'craterscale.dat', - 'craterlist' : 'craterlist.dat' + 'craterlist' : 'craterlist.dat', + 'sfdfile' : 'production.dat' } + if self.user['sfdfile'] is not None: # Override the default sfdfile name if the user supplies an alternative + self.output_filenames['sfdfile'] = self.user['sfdfile'] for k, v in self.output_filenames.items(): self.output_filenames[k] = os.path.join(currentdir, v) @@ -120,29 +123,18 @@ def __init__(self, param_file="ctem.in", isnew=True): #Read list of real craters print("quasi-MC mode is ON") + print("Generating the crater scaling data in CTEM") rclist = util.read_formatted_ascii(self.user['realcraterlist'], skip_lines = 0) - self.tempfile = os.path.join(currentdir, 'temp.in') + tempfile = os.path.join(currentdir, 'temp.in') # Generate craterlist.dat - with open(self.user['ctemfile']) as perm, open(self.tempfile, "w") as temp: - for line in perm: - temp.write(line) + shutil.copy2(self.user['ctemfile'], tempfile ) #Write a temporary input file to generate the necessary quasimc files - self.permfile = self.user['ctemfile'] - self.temp = util.write_temp_input(self.tempfile) - self.user['ctemfile'] = self.tempfile - self.user = util.read_user_input(self.user) - self.permin = os.path.join(currentdir, 'perm.in') - self.ctemin = os.path.join(currentdir, 'ctem.in') - os.rename(self.permfile, self.permin) - os.rename(self.tempfile, self.ctemin) + util.write_temp_input(tempfile) util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr) - self.compute_one_interval() - - os.replace(self.permin, self.ctemin) - self.user['ctemfile'] = self.permfile - self.user = util.read_user_input(self.user) + self.compute_one_interval(ctemin=tempfile) + os.remove(tempfile) #Interpolate craterscale.dat to get impactor sizes from crater sizes given df = pandas.read_csv(self.output_filenames['craterscale'], sep='\s+') @@ -183,7 +175,7 @@ def scale_production(self): self.production[:, 1] = self.production[:, 1] * area * self.user['interval'] # Write the scaled production function to file - util.write_production(self.user, self.production) + util.write_production(self.output_filenames['sfdfile'], self.production) return @@ -202,7 +194,7 @@ def run(self): self.redirect_outputs(['dat'], 'misc') #Execute Fortran code - self.compute_one_interval() + self.compute_one_interval(self.user['ctemfile']) # Read in output files self.read_output() @@ -352,13 +344,6 @@ def cleanup(self): os.remove(filename) except OSError as error: print(error) - # Delete scaled production file - prodfile = self.user['sfdfile'] - print(f"Deleting file {prodfile}") - try: - os.remove(prodfile) - except OSError as error: - print(error) return diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index b20bc08a..5d2c1700 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -243,8 +243,6 @@ def read_user_input(user): print('Invalid value for or missing variable GRIDSIZE in ' + inputfile) if (user['seed'] == 0): print('Invalid value for or missing variable SEED in ' + inputfile) - if (user['sfdfile'] is None): - print('Invalid value for or missing variable SFDFILE in ' + inputfile) if (user['impfile'] is None): print('Invalid value for or missing variable IMPFILE in ' + inputfile) if (user['popupconsole'] is None): @@ -336,8 +334,7 @@ def write_datfile(user, filename, seedarr): return -def write_production(user, production): - filename = user['sfdfile'] +def write_production(filename, production): np.savetxt(filename, production, fmt='%1.8e', delimiter=' ') return diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 6157dc9d..f29b7600 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -270,7 +270,7 @@ module module_globals integer(I4B) :: pbarpos character(len=PBARSIZE) :: pbarchar -! Grid array file names +! Default file names character(*),parameter :: DIAMFILE = 'surface_diam.dat' character(*),parameter :: TIMEFILE = 'surface_time.dat' character(*),parameter :: EJCOVFILE = 'surface_ejc.dat' @@ -284,7 +284,6 @@ module module_globals character(*),parameter :: POROFILE = 'porosity_porosity.dat' character(*),parameter :: DEPTHFILE = 'porosity_depth.dat' character(*),parameter :: POROTHICKFILE = 'porosity_thickness.dat' -!character(*),parameter :: THICKFILE = 'surface_crustal_thickness.dat' character(*),parameter :: POSFILE = 'surface_pos.dat' character(*),parameter :: TDISTFILE = 'tdistribution.dat' character(*),parameter :: TLISTFILE = 'tcumulative.dat' @@ -292,7 +291,9 @@ module module_globals character(*),parameter :: OLISTFILE = 'ocumulative.dat' character(*),parameter :: PDISTFILE = 'pdistribution.dat' character(*),parameter :: CRTSCLFILE = 'craterscale.dat' -character(*),parameter :: DATFILE = 'ctem.dat' +character(*),parameter :: SFDFILE = 'production.dat' +character(*),parameter :: USERFILE = 'ctem.in' +character(*),parameter :: DATFILE = 'ctem.dat' ! character(*),parameter :: MASSFILE = 'impactmass.dat' character(*),parameter :: RCFILE = 'craterlist.dat' !not sure if this is where this line should go, but putting it here for now... diff --git a/src/io/io_input.f90 b/src/io/io_input.f90 index 57995ed7..9c56298e 100644 --- a/src/io/io_input.f90 +++ b/src/io/io_input.f90 @@ -31,7 +31,7 @@ subroutine io_input(infile,user) integer(I4B), parameter :: LUN = 7 integer(I4B) :: ierr, ilength, ifirst, ilast,i character(STRMAX) :: line, token - integer(I4B), parameter :: numrequired=18 + integer(I4B), parameter :: numrequired=17 character(STRMAX),dimension(numrequired),parameter :: requiredvar = (/"GRIDSIZE ",& "NUMLAYERS", & "PIX ", & @@ -48,7 +48,6 @@ subroutine io_input(infile,user) "TRHO_B ", & "MAT ", & "PRHO ", & - "SFDFILE ", & "VELFILE "/) integer(I4B), parameter :: seismic_numrequired=5 @@ -94,6 +93,7 @@ subroutine io_input(infile,user) user%dorealistic = .false. user%doquasimc = .false. user%ejecta_truncation = 10.0_DP + write(user%sfdfile,*) trim(adjustl(SFDFILE)) open(unit=LUN,file=infile,status="old",iostat=ierr) if (ierr /= 0) then @@ -212,14 +212,13 @@ subroutine io_input(infile,user) ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) - read(token, *) user%sfdfile + read(token, *) user%velfile ismissing(17)=.false. - case (trim(requiredvar(18))) + case ("SFDFILE") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) token = line(ifirst:ilast) - read(token,*) user%velfile - ismissing(18)=.false. + read(token, *) user%sfdfile case ("DEPLIMIT") ifirst = ilast + 1 call io_get_token(line, ilength, ifirst, ilast, ierr) diff --git a/src/main/CTEM.f90 b/src/main/CTEM.f90 index db31e955..da0ea979 100644 --- a/src/main/CTEM.f90 +++ b/src/main/CTEM.f90 @@ -69,9 +69,9 @@ program CTEM call io_splash() narg = command_argument_count() if (narg /= 0) then - call get_command_argument(1, infile, status = ierr) + call get_command_argument(1, infile, status = ierr) ! Use first argument as the user input file name else - infile="ctem.in" + infile = USERFILE ! No arguments, so use the default file name for the user inputs end if write(*,*) 'Reading input file ',trim(adjustl(infile)) call io_input(infile,user) From f5636b972f823ebc828c31373bebad199e83ab9d Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Mar 2022 19:31:57 -0500 Subject: [PATCH 06/25] Removed obsolete files --- examples/quasimc-global/craterproduction.py | 216 -------------- examples/quasimc-global/ctem_driver.pro | 233 --------------- examples/quasimc-global/ctem_driver.py | 256 ---------------- examples/quasimc-global/ctem_image_dem.pro | 41 --- .../ctem_image_presentation.pro | 143 --------- .../quasimc-global/ctem_image_regolith.pro | 33 -- .../ctem_image_shaded_relief.pro | 52 ---- .../quasimc-global/ctem_io_read_input.pro | 131 -------- examples/quasimc-global/ctem_io_read_old.pro | 45 --- examples/quasimc-global/ctem_io_readers.py | 146 --------- examples/quasimc-global/ctem_io_writers.py | 282 ------------------ .../quasimc-global/ctem_window_display.pro | 249 ---------------- 12 files changed, 1827 deletions(-) delete mode 100644 examples/quasimc-global/craterproduction.py delete mode 100755 examples/quasimc-global/ctem_driver.pro delete mode 100755 examples/quasimc-global/ctem_driver.py delete mode 100755 examples/quasimc-global/ctem_image_dem.pro delete mode 100755 examples/quasimc-global/ctem_image_presentation.pro delete mode 100755 examples/quasimc-global/ctem_image_regolith.pro delete mode 100755 examples/quasimc-global/ctem_image_shaded_relief.pro delete mode 100755 examples/quasimc-global/ctem_io_read_input.pro delete mode 100755 examples/quasimc-global/ctem_io_read_old.pro delete mode 100644 examples/quasimc-global/ctem_io_readers.py delete mode 100644 examples/quasimc-global/ctem_io_writers.py delete mode 100755 examples/quasimc-global/ctem_window_display.pro diff --git a/examples/quasimc-global/craterproduction.py b/examples/quasimc-global/craterproduction.py deleted file mode 100644 index 4d81b687..00000000 --- a/examples/quasimc-global/craterproduction.py +++ /dev/null @@ -1,216 +0,0 @@ -import numpy as np -from scipy import optimize - -# The Neukum production function for the Moon and Mars -#Lunar PF from: Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86 -#Mars PF from: Ivanov (2001) SSR v. 96 pp. 87-104 -sfd_coef = { - "NPF_Moon" : [-3.0876,-3.557528,+0.781027,+1.021521,-0.156012,-0.444058,+0.019977,+0.086850,-0.005874,-0.006809,+8.25e-4, +5.54e-5], - "NPF_Mars" : [-3.384, -3.197, +1.257, +0.7915, -0.4861, -0.3630, +0.1016, +6.756e-2,-1.181e-2,-4.753e-3,+6.233e-4,+5.805e-5] - } -sfd_range = { - "NPF_Moon" : [0.01,1000], - "NPF_Mars" : [0.015,362] -} -#Exponential time constant (Ga) -tau = 6.93 -Nexp = 5.44e-14 - -time_coef = { - "NPF_Moon" : Nexp, - "NPF_Mars" : Nexp * 10**(sfd_coef.get("NPF_Mars")[0]) / 10**(sfd_coef.get("NPF_Moon")[0]) -} - -def N1(T, pfmodel): - """ - Return the N(1) value as a function of time for a particular production function model - - Parameters - ---------- - T : numpy array - Time in units of Ga - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - N1 : numpy array - The number of craters per square kilometer greater than 1 km in diameter - """ - T_valid_range = np.where((T <= 4.5) & (T >= 0.0), T, np.nan) - return time_coef.get(pfmodel) * (np.exp(tau * T_valid_range) - 1.0) + 10 ** (sfd_coef.get(pfmodel)[0]) * T_valid_range - -def CSFD(Dkm,planet): - """ - Return the cumulative size-frequency distribution for a particular production function model - - Parameters - ---------- - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - CSFD : numpy array - The number of craters per square kilometer greater than Dkm in diameter at T=1 Ga - """ - logCSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(sfd_coef.get(planet))) - return 10**(logCSFD) - -def DSFD(Dkm,planet): - """ - Return the differential size-frequency distribution for a particular production function model - - Parameters - ---------- - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - DSFD : numpy array - The differential number of craters (dN/dD) per square kilometer greater than Dkm in diameter at T = 1 Ga - """ - dcoef = sfd_coef.get(planet)[1:] - logDSFD = sum(co * np.log10(Dkm) ** i for i, co in enumerate(dcoef)) - return 10**(logDSFD) * CSFD(Dkm,planet) / Dkm - -def Tscale(T,planet): - """ - Return the number density of craters at time T relative to time T = 1 Ga - - Parameters - ---------- - T : numpy array - Time in units of Ga - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - Tscale : numpy array - N1(T) / CSFD(Dkm = 1.0) - """ - return N1(T,planet) / CSFD(1.0,planet) - -def pf_csfd(T,Dkm,planet): - """ - Return the cumulative size-frequency distribution for a particular production function model as a function of Time - - Parameters - ---------- - T : numpy array - Time in units of Ga - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - pf_csfd : numpy array - The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T - """ - D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]),Dkm,np.nan) - return CSFD(D_valid_range,planet) * Tscale(T,planet) - -def pf_dsfd(T,Dkm,planet): - """ - Return the differential size-frequency distribution for a particular production function model as a function of Time - - Parameters - ---------- - T : numpy array - Time in units of Ga - Dkm : numpy array - Diameters in units of km - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - pf_dsfd : numpy array - The cumulative number of craters per square kilometer greater than Dkm in diameter at time T T - """ - D_valid_range = np.where((Dkm >= sfd_range.get(planet)[0]) & (Dkm <= sfd_range.get(planet)[1]), Dkm, np.nan) - return DSFD(D_valid_range, planet) * Tscale(T, planet) - -def T_from_scale(TS,planet): - """ - Return the time in Ga for the given number density of craters relative to that at 1 Ga. - This is the inverse of Tscale - - Parameters - ---------- - TS : numpy array - number density of craters relative to that at 1 Ga - pfmodel : string - The production function model to use - Currently options are 'NPF_Moon' or 'NPF_Mars' - - Returns - ------- - T_from_scale : numpy array - The time in Ga - """ - def func(S): - return Tscale(S,planet) - TS - return optimize.fsolve(func, np.full_like(TS,4.4),xtol=1e-10) - - -if __name__ == "__main__": - import matplotlib.pyplot as plt - import matplotlib.ticker as ticker - print("Tests go here") - print(f"T = 1 Ga, N(1) = {pf_csfd(1.0,1.00,'NPF_Moon')}") - print(f"T = 4.2 Ga, N(1) = {pf_csfd(4.2,1.00,'NPF_Moon')}") - print("Tscale test: Should return all 1s") - Ttest = np.logspace(-4,np.log10(4.4),num=100) - Tres = T_from_scale(Tscale(Ttest,'NPF_Mars'),'NPF_Mars') - print(Ttest / Tres) - #for i,t in enumerate(Ttest): - # print(t,Tscale(t,'Moon'),Tres[i]) - - CSFDfig = plt.figure(1, figsize=(8, 7)) - ax = {'NPF_Moon': CSFDfig.add_subplot(121), - 'NPF_Mars': CSFDfig.add_subplot(122)} - - tvals = [0.01,1.0,4.0] - x_min = 1e-3 - x_max = 1e3 - y_min = 1e-9 - y_max = 1e3 - Dvals = np.logspace(np.log10(x_min), np.log10(x_max)) - for key in ax: - ax[key].title.set_text(key) - ax[key].set_xscale('log') - ax[key].set_yscale('log') - ax[key].set_ylabel('$\mathregular{N_{>D} (km^{-2})}$') - ax[key].set_xlabel('Diameter (km)') - ax[key].set_xlim(x_min, x_max) - ax[key].set_ylim(y_min, y_max) - ax[key].yaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20)) - ax[key].yaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100)) - ax[key].xaxis.set_major_locator(ticker.LogLocator(base=10.0, numticks=20)) - ax[key].xaxis.set_minor_locator(ticker.LogLocator(base=10.0, subs=np.arange(2,10), numticks=100)) - ax[key].grid(True,which="minor",ls="-",lw=0.5,zorder=5) - ax[key].grid(True,which="major",ls="-",lw=1,zorder=10) - for t in tvals: - prod = pf_csfd(t,Dvals,key) - ax[key].plot(Dvals, prod, '-', color='black', linewidth=1.0, zorder=50) - labeli = 15 - ax[key].text(Dvals[labeli],prod[labeli],f"{t:.2f} Ga", ha="left", va="top",rotation=-72) - - plt.tick_params(axis='y', which='minor') - plt.tight_layout() - plt.show() diff --git a/examples/quasimc-global/ctem_driver.pro b/examples/quasimc-global/ctem_driver.pro deleted file mode 100755 index 060d1b8d..00000000 --- a/examples/quasimc-global/ctem_driver.pro +++ /dev/null @@ -1,233 +0,0 @@ -pro ctem_driver - -;------------------------------------------------------------------------- -; Jim Richardson, Arecibo Observatory -; David Minton, Purdue University Dept. of Earth, Atmospheric, & Planetary Sciences -; July 2014 -; Cratered Terrain Evolution Model display module -; -; Inputs are read in from the ctem.in file -; -;------------------------------------------------------------------------- -;------------- Initial Setup ---------------- -;------------------------------------------------------------------------- -Compile_Opt DEFINT32 -!EXCEPT=2 - -; ----------- input file ----------------------- -infilename = 'ctem.in' -DATFILE='ctem.dat' - -; ---------- reading input files ---------- -seedarr = lon64arr(100) -seedn = 1 -totalimpacts = long64(0) -ncount = long64(0) -curyear = 0.d0 -restart = "F" -fracdone = 1.0d0 -masstot = 0.d0 - -ctem_io_read_input,infilename,interval,numintervals,gridsize,pix,seed,numlayers,sfdfile,impfile,maxcrat,ph1,shadedmaxhdefault,shadedminhdefault,shadedminh,shadedmaxh,restart,runtype,popupconsole,saveshaded,saverego,savepres,savetruelist - -seedarr(0) = seed -area = (gridsize * pix)^2 - -;read input data -pnum = file_lines(impfile) -production = dblarr(2,pnum) -productionfunction = dblarr(2,pnum) -openr,LUN,impfile,/GET_LUN -readf,LUN,productionfunction -free_lun,LUN - -;create impactor production population -production(0,*) = productionfunction(0,*) -production(1,*) = productionfunction(1,*)*area*interval - -;write out corrected production population -openw,1,sfdfile -printf,1,production -close,1 -free_lun,1 - -;set up cratering surface grid and display-only grid -surface_dem = dblarr(gridsize,gridsize) -regolith = dblarr(gridsize,gridsize) - -;set up temporary distribution bins -distl = 1 -pdistl = 1 -odist = dblarr(6,distl) -tdist = dblarr(6,distl) -pdist = dblarr(6,pdistl) -pdisttotal = dblarr(6,pdistl) - -datformat = "(I17,1X,I12,1X,E19.12,1X,A1,1X,F9.6,1X,E19.12)" - - -if strmatch(restart,'F',/fold_case) then begin ; Start with a clean slate - print, 'Starting a new run' - curyear = 0.0d0 - totalimpacts = 0 - masstot = 0.d0 - fracdone = 1.0d0 - - if strmatch(runtype,'statistical',/fold_case) then begin - ncount = 1 - openw,LUN,DATFILE,/GET_LUN - printf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - for n=0,seedn-1 do begin - printf,LUN,seedarr(n),format='(I12)' - endfor - free_lun,LUN - endif else begin - ncount = 0 - endelse - - surface_dem(*,*) = 0.0d0 - regolith(*,*) = 0.0d0 - - file_delete, 'tdistribution.dat',/allow_nonexistent - -endif else begin ; continue an old run - print, 'Continuing a previous run' - - ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass - - ;read in constants file - openr,LUN,DATFILE,/GET_LUN - readf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - seedn = 0 - while ~ eof(LUN) do begin - readf,LUN,iseed - seedarr(seedn) = long64(iseed) - seedn=seedn+1 - endwhile - -endelse - -openw,FRACDONEFILE,'fracdone.dat',/GET_LUN -openw,REGODEPTHFILE,'regolithdepth.dat',/GET_LUN - -;------------------------------------------------------------------------- -; ---------- begin loops ---------- -;------------------------------------------------------------------------- -print, 'Beginning loops' - -;define number of loop iterations and begin -;numintervals=ceil(endyear/interval)-ncount -while (ncount le numintervals) do begin - - - ; ---------- creating crater population ---------- - if (ncount gt 0) then begin - - fnum = string(ncount,format='(I6.6)') - if (file_test('misc',/DIRECTORY) eq 0) then begin - file_mkdir,'misc' - endif - ; save a copy of the ctem.dat file - fname = 'misc/ctem_' + fnum + '.dat' - file_copy, 'ctem.dat', fname, /OVERWRITE - - print, ncount, ' Calling FORTRAN routine' - ;call fortran program to create & count craters - spawn, './CTEM',/noshell - - ; ---------- reading FORTRAN output ---------- - print, ncount, ' Reading FORTRAN output' - ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass - - ;read in constants file - openr,LUN,DATFILE,/GET_LUN - readf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - seedn = 0 - while ~ eof(LUN) do begin - readf,LUN,iseed - seedarr(seedn) = long64(iseed) - seedn = seedn + 1 - endwhile - free_lun,LUN - - curyear = curyear + fracdone * interval - masstot = masstot + mass - printf,FRACDONEFILE,fracdone,curyear - flush,FRACDONEFILE - - printf,REGODEPTHFILE,curyear,mean(regolith),max(regolith),min(regolith) - flush,REGODEPTHFILE - - ;save a copy of the binned observed crater distribution - if (file_test('dist',/DIRECTORY) eq 0) then begin - file_mkdir,'dist' - endif - fname = 'dist/odist_' + fnum + '.dat' - file_copy, 'odistribution.dat', fname, /OVERWRITE - - ; save a copy of the cumulative observed crater distribution - fname = 'dist/ocum_' + fnum + '.dat' - file_copy, 'ocumulative.dat', fname, /OVERWRITE - - ;save a copy of the binned true distribution - fname = 'dist/tdist_' + fnum + '.dat' - file_copy, 'tdistribution.dat', fname, /OVERWRITE - - ;save a copy of the binned idealized production function - fname = 'dist/pdist_' + fnum + '.dat' - file_copy, 'pdistribution.dat', fname, /OVERWRITE - - ; save a copy of the cumulative true crater distribution if the user requests it - if strmatch(savetruelist,'T',/fold_case) then begin - fname = 'dist/tcum_' + fnum + '.dat' - file_copy, 'tcumulative.dat', fname, /OVERWRITE - endif - - ; save a copy of the impacted mass - fname = 'misc/mass_' + fnum + '.dat' - file_copy, 'impactmass.dat', fname, /OVERWRITE - - - endif - - ; Get the accumulated production function - pdisttotal = pdist - pdisttotal(3:5,*) = pdist(3:5,*) * curyear / interval - - ; ---------- displaying results ---------- - print, ncount, ' Displaying results' - - ctem_image_dem,ncount,gridsize,pix,surface_dem,surface_dem_image - if strmatch(saverego,'T',/fold_case) then ctem_image_regolith,ncount,gridsize,pix,regolith,regolith_image - if strmatch(saveshaded,'T',/fold_case) then begin - if (shadedminhdefault eq 1) then shadedminh = min(surface_dem) - if (shadedmaxhdefault eq 1) then shadedmaxh = max(surface_dem) - ctem_image_shaded_relief,ncount,gridsize,pix,surface_dem,surface_dem,shadedminh,shadedmaxh,'shaded',shaded_image - endif - if strmatch(savepres,'T',/fold_case) then ctem_image_presentation,ncount,gridsize,pix,curyear,odist,pdisttotal,tdist,ph1,surface_dem_image - ctem_window_display,ncount,totalimpacts,gridsize,pix,curyear,masstot,odist,pdisttotal,tdist,ph1,surface_dem,regolith,surface_dem_image,popupconsole - - ncount = ncount + 1 - - ;write out the current data file - if (strmatch(runtype,'statistical',/fold_case)) || (ncount eq 1) then begin - restart = 'F' - curyear = 0.0d0 - totalimpacts = 0 - masstot = 0.d0 - file_delete, 'tdistribution.dat',/allow_nonexistent - endif else begin - restart = 'T' - endelse - - openw,LUN,DATFILE,/GET_LUN - printf,LUN,totalimpacts,ncount,curyear,restart,fracdone,masstot,format=datformat - for n=0,seedn-1 do begin - printf,LUN,seedarr(n),format='(I12)' - endfor - free_lun,LUN -endwhile -free_lun,FRACDONEFILE -free_lun,REGODEPTHFILE - -end diff --git a/examples/quasimc-global/ctem_driver.py b/examples/quasimc-global/ctem_driver.py deleted file mode 100755 index 6fb56613..00000000 --- a/examples/quasimc-global/ctem_driver.py +++ /dev/null @@ -1,256 +0,0 @@ -#!/usr/bin/env python3 -# -#Cratered Terrain Evolution Model driver -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 - -#Import general purpose modules - -import numpy -import os -import subprocess -import shutil -import pandas -from scipy.interpolate import interp1d - -#Import CTEM modules -import ctem_io_readers -import ctem_io_writers -import craterproduction #craterproduction had to be cp'd to example dir - -#Create and initialize data dictionaries for parameters and options from CTEM.in -notset = '-NOTSET-' -currentdir = os.getcwd() + os.sep - -parameters={'restart': notset, - 'runtype': notset, - 'popupconsole': notset, - 'saveshaded': notset, - 'saverego': notset, - 'savepres': notset, - 'savetruelist': notset, - 'seedn': 1, - 'totalimpacts': 0, - 'ncount': 0, - 'curyear': 0.0, - 'fracdone': 1.0, - 'masstot': 0.0, - 'interval': 0.0, - 'numintervals': 0, - 'pix': -1.0, - 'gridsize': -1, - 'seed': 0, - 'maxcrat': 1.0, - 'shadedminhdefault': 1, - 'shadedmaxhdefault': 1, - 'shadedminh': 0.0, - 'shadedmaxh': 0.0, - 'workingdir': currentdir, - 'ctemfile': 'ctem.in', - 'datfile': 'ctem.dat', - 'impfile': notset, - 'sfdcompare': notset, - 'sfdfile': notset, - 'quasimc': notset, - 'realcraterlist': notset} - -#Read ctem.in to initialize parameter values based on user input -ctem_io_readers.read_ctemin(parameters,notset) - -#Read sfdcompare file -sfdfile = parameters['workingdir'] + parameters['sfdcompare'] -ph1 = ctem_io_readers.read_formatted_ascii(sfdfile, skip_lines = 0) - -#Set up data arrays -seedarr = numpy.zeros(100, dtype = numpy.int) -seedarr[0] = parameters['seed'] -odist = numpy.zeros([1, 6]) -pdist = numpy.zeros([1, 6]) -tdist = numpy.zeros([1, 6]) -surface_dem = numpy.zeros([parameters['gridsize'], parameters['gridsize']], dtype = numpy.float) -regolith = numpy.zeros([parameters['gridsize'], parameters['gridsize']], dtype =numpy.float) - -#Read production function file -impfile = parameters['workingdir'] + parameters['impfile'] -prodfunction = ctem_io_readers.read_formatted_ascii(impfile, skip_lines = 0) - -if (parameters['quasimc'] == 'T'): - - #Read list of real craters - print("quasi-MC mode is ON") - craterlistfile = parameters['workingdir'] + parameters['realcraterlist'] - rclist = ctem_io_readers.read_formatted_ascii(craterlistfile, skip_lines = 0) - - #Interpolate craterscale.dat to get impactor sizes from crater sizes given - df = pandas.read_csv('craterscale.dat', sep='\s+') - df['log(Dc)'] = numpy.log(df['Dcrat(m)']) - df['log(Di)'] = numpy.log(df['#Dimp(m)']) - xnew = df['log(Dc)'].values - ynew = df['log(Di)'].values - interp = interp1d(xnew, ynew, fill_value='extrapolate') - rclist[:,0] = numpy.exp(interp(numpy.log(rclist[:,0]))) - - #Convert age in Ga to "interval time" - rclist[:,5] = (parameters['interval'] * parameters['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon') - rclist = rclist[rclist[:,5].argsort()] - - #Export to dat file for Fortran use - ctem_io_writers.write_realcraters(parameters, rclist) - -#Create impactor production population -area = (parameters['gridsize'] * parameters['pix'])**2 -production = numpy.copy(prodfunction) -production[:,1] = production[:,1] * area * parameters['interval'] - -#Write corrected production function to file -ctem_io_writers.write_production(parameters, production) - -#Starting new or old run? -if (parameters['restart'].upper() == 'F'): - print('Starting a new run') - - if (parameters['runtype'].upper() == 'STATISTICAL'): - parameters['ncount'] = 1 - - #Write ctem.dat file - ctem_io_writers.write_ctemdat(parameters, seedarr) - - else: - parameters['ncount'] = 0 - - #Delete tdistribution file, if it exists - tdist_file = parameters['workingdir'] + 'tdistribution.dat' - if os.path.isfile(tdist_file): - os.remove(tdist_file) - -else: - print('Continuing a previous run') - - #Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + 'surface_dem.dat' - surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - ejecta_file = parameters['workingdir'] + 'surface_ejc.dat' - regolith = ctem_io_readers.read_unformatted_binary(ejecta_file, parameters['gridsize']) - - #Read odistribution, tdistribution, and pdistribution files - ofile = parameters['workingdir'] + 'odistribution.dat' - odist = ctem_io_readers.read_formatted_ascii(ofile, skip_lines = 1) - tfile = parameters['workingdir'] + 'tdistribution.dat' - tdist = ctem_io_readers.read_formatted_ascii(tfile, skip_lines = 1) - pfile = parameters['workingdir'] + 'pdistribution.dat' - pdist = ctem_io_readers.read_formatted_ascii(pfile, skip_lines = 1) - - #Read impact mass from file - massfile = parameters['workingdir'] + 'impactmass.dat' - impact_mass = ctem_io_readers.read_impact_mass(massfile) - - #Read ctem.dat file - ctem_io_readers.read_ctemdat(parameters, seedarr) - -#Open fracdonefile and regodepthfile for writing -filename = parameters['workingdir'] + 'fracdone.dat' -fp_frac = open(filename,'w') -filename = parameters['workingdir'] + 'regolithdepth.dat' -fp_reg = open(filename,'w') - -#Begin CTEM processing loops -print('Beginning loops') - -ctem_io_writers.create_dir_structure(parameters) - -while (parameters['ncount'] <= parameters['numintervals']): - - #Create crater population - if (parameters['ncount'] > 0): - - #Move ctem.dat - forig = parameters['workingdir'] + 'ctem.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "ctem%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - #Create crater population and display CTEM progress on screen - print(parameters['ncount'], ' Calling FORTRAN routine') - with subprocess.Popen([parameters['workingdir']+'CTEM'], stdout=subprocess.PIPE, bufsize=1,universal_newlines=True) as p: - for line in p.stdout: - print(line, end='') - - - #Optional: do not pipe CTEM progress to the screen - #subprocess.check_output([parameters['workingdir']+'CTEM']) - - #Read Fortran output - print(parameters['ncount'], ' Reading Fortran output') - - #Read surface dem(shaded relief) and ejecta data files - dem_file = parameters['workingdir'] + 'surface_dem.dat' - surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize']) - ejecta_file = parameters['workingdir'] + 'surface_ejc.dat' - regolith = ctem_io_readers.read_unformatted_binary(ejecta_file, parameters['gridsize']) - - #Read odistribution, tdistribution, and pdistribution files - ofile = parameters['workingdir'] + 'odistribution.dat' - odist = ctem_io_readers.read_formatted_ascii(ofile, skip_lines = 1) - tfile = parameters['workingdir'] + 'tdistribution.dat' - tdist = ctem_io_readers.read_formatted_ascii(tfile, skip_lines = 1) - pfile = parameters['workingdir'] + 'pdistribution.dat' - pdist = ctem_io_readers.read_formatted_ascii(pfile, skip_lines = 1) - - #Read impact mass from file - massfile = parameters['workingdir'] + 'impactmass.dat' - impact_mass = ctem_io_readers.read_impact_mass(massfile) - - #Read ctem.dat file - ctem_io_readers.read_ctemdat(parameters, seedarr) - - #Update parameters: mass, curyear, regolith properties - parameters['masstot'] = parameters['masstot'] + impact_mass - - parameters['curyear'] = parameters['curyear'] + parameters['fracdone'] * parameters['interval'] - template = "%(fracdone)9.6f %(curyear)19.12E\n" - fp_frac.write(template % parameters) - - reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (parameters['curyear'], - numpy.mean(regolith), numpy.amax(regolith), numpy.amin(regolith)) - fp_reg.write(reg_text) - - #Save copy of crater distribution files - ctem_io_writers.copy_dists(parameters) - - #Display results - print(parameters['ncount'], ' Displaying results') - - #Write surface dem, surface ejecta, shaded relief, and rplot data - ctem_io_writers.image_dem(parameters, surface_dem) - if (parameters['saverego'].upper() == 'T'): - ctem_io_writers.image_regolith(parameters, regolith) - if (parameters['saveshaded'].upper() == 'T'): - ctem_io_writers.image_shaded_relief(parameters, surface_dem) - if (parameters['savepres'].upper() == 'T'): - ctem_io_writers.create_rplot(parameters,odist,pdist,tdist,ph1) - - #Update ncount - parameters['ncount'] = parameters['ncount'] + 1 - - if ((parameters['runtype'].upper() == 'STATISTICAL') or (parameters['ncount'] == 1)): - parameters['restart'] = 'F' - parameters['curyear'] = 0.0 - parameters['totalimpacts'] = 0 - parameters['masstot'] = 0.0 - - #Delete tdistribution file, if it exists - tdist_file = parameters['workingdir'] + 'tdistribution.dat' - if os.path.isfile(tdist_file): - os.remove(tdist_file) - - else: - parameters['restart'] = 'T' - - #Write ctem.dat file - ctem_io_writers.write_ctemdat(parameters, seedarr) - -#Close updateable fracdonefile and regodepthfile files -fp_frac.close() -fp_reg.close() diff --git a/examples/quasimc-global/ctem_image_dem.pro b/examples/quasimc-global/ctem_image_dem.pro deleted file mode 100755 index ad6fbd7e..00000000 --- a/examples/quasimc-global/ctem_image_dem.pro +++ /dev/null @@ -1,41 +0,0 @@ -pro ctem_image_dem,ncount,gridsize,pix,surface_dem,surface_dem_image -; Generates the shaded surface DEM image and saves it as a jpeg output image in the 'surf' directory -; outputs surface_dem_arr, which may be scaled to use as a console image - -; This code is for reading in the x-y positions from the cumulative distribution and separating out into layers -; so that the tallied craters can be drawn as circles -;!PATH = Expand_Path('+/home/campus/daminton/coyote/') + ':' + !PATH - -Compile_Opt DEFINT32 -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[gridsize,gridsize],Set_Pixel_Depth=24, Decomposed=0 -TVLCT, red, green, blue, /GET - -sun = 20.d0 -radsun = sun * !dtor - -surface_dem_arr = dblarr(gridsize,gridsize) - -surface_dem_arr=surface_dem-shift(surface_dem,0,1) -surface_dem_arr = (0.5d0*!dpi) + atan(surface_dem_arr,pix) ; Get average slope -surface_dem_arr = abs(surface_dem_arr - radsun < 0.5d0*!dpi) ; incident angle -surface_dem_arr = 254.0d0 * cos(surface_dem_arr) ; shaded relief surface - -tv, surface_dem_arr, 0, 0, xsize=gridsize, ysize=gridsize, /device -print,'plotting' - -surface_dem_image = TVRD(True=1) -Set_Plot, thisDevice - -; save surface display -if (file_test('surf',/DIRECTORY) eq 0) then begin - file_mkdir,'surf' -endif -fnum = string(ncount,format='(I6.6)') -fname = 'surf/surf' + fnum + '.jpg' - -write_jpeg, fname, surface_dem_image, true=1, quality=100 - -end diff --git a/examples/quasimc-global/ctem_image_presentation.pro b/examples/quasimc-global/ctem_image_presentation.pro deleted file mode 100755 index 05a7d65a..00000000 --- a/examples/quasimc-global/ctem_image_presentation.pro +++ /dev/null @@ -1,143 +0,0 @@ -pro ctem_image_presentation,ncount,gridsize,pix,curyear,odist,pdist,tdist,ph1,map -; Generates the simplified version of the console display for use in presentations -; Plots the R-plot on the left and the image map on the right - -; presentation graphics -presresx = 1024 -presresy = 0.6*presresx -area = (gridsize * pix)^2 - -minx = (pix / 3.0d0) * 1d-3 -maxx = 3 * pix * gridsize * 1d-3 - -;strings for display -dum = string(format='(E10.4)', curyear) -parts = strsplit(dum, 'E', /extract) -fcuryear = strcompress(string(format='(F6.4,"x10!U",i,"!N")', parts[0], parts[1])) -time = 'Time = ' -timeunit = ' yr' - -if (file_test('presentation',/DIRECTORY) eq 0) then begin - file_mkdir,'presentation' -endif - -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[presresx,presresy],Set_Pixel_Depth=24, Decomposed=0 -loadct, 0 -TVLCT, red, green, blue, /GET - -;geometric saturation -geom = dblarr(2,2) -geomem = dblarr(2,2) -geomep = dblarr(2,2) -geomel = dblarr(2,2) -geom(0,0) = minx -geom(0,1) = maxx -geom(1,*) = 3.12636d0 -geomem(0,0) = minx -geomem(0,1) = maxx -geomem(1,*) = 0.156318d0 -geomep(0,0) = minx -geomep(0,1) = maxx -geomep(1,*) = 0.312636d0 -geomel(0,0) = minx -geomel(0,1) = maxx -geomel(1,*) = 0.0312636d0 - -; Remove zeros -nz = where(odist(5,*) ne 0.0,count) -if (count gt 0) then begin - odistnz = dblarr(6,count) - odistnz(*,*)= odist(*,nz) -endif else begin - odistnz = odist -endelse - -nz = where(tdist(5,*) ne 0.0,count) -if (count gt 0) then begin - tdistnz = dblarr(6,count) - tdistnz(*,*)= tdist(*,nz) -endif else begin - tdistnz = tdist -endelse - -nz = where(pdist(5,*) ne 0.0,count) -if (count gt 0) then begin - pdistnz = dblarr(6,count) - pdistnz(*,*)= pdist(*,nz) -endif else begin - pdistnz = pdist -endelse - -; create r-plot array containing exactly 1 crater per bin -area = (gridsize * pix * 1d-3)^2 -plo = 1 -while (sqrt(2.0d0)^plo gt minx) do begin - plo = plo - 1 -endwhile -phi = plo + 1 -while (sqrt(2.0d0)^phi lt maxx) do begin - phi = phi + 1 -endwhile -n = phi - plo -sdist = dblarr(6,n + 1) -p = plo -for i=0,n do begin - sdist(0,i) = sqrt(2.d0)^p - sdist(1,i) = sqrt(2.d0)^(p+1) - sdist(2,i) = sqrt(sdist(0,i) * sdist(1,i)) - sdist(3,i) = 1.0d0 - sdist(5,i) = (sdist(2,i))^3 / (area * (sdist(1,i) - sdist(0,i))) - p = p + 1 -endfor -sdist(4,*) = reverse(total(sdist(3,*),/cumulative),2) - - -plot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=0, color=255, thick=2.0, $ - TITLE='Crater Distribution R-Plot', charsize=1.2, $ - XTITLE='Crater Diameter (km)', $ - XRANGE=[minx,maxx], /XLOG, XSTYLE=1, $ - YTITLE='R Value', $ - YRANGE=[5.0e-4,5.0e0], /YLOG, YSTYLE=1, $ - /device, pos = [0.12,0.12,0.495,0.495]*presresx - -Dfac = sqrt(sqrt(2.0d0)) - -;display observed crater distribution -oplot, tdistnz(2,*)*1.0d-3, tdistnz(5,*), line=0, color=255, thick=1.0 -;oplot, geom(0,*), geom(1,*), line=2, color=255, thick=1.0 -oplot, geomem(0,*), geomem(1,*), line=1, color=255, thick=1.0 -oplot, geomep(0,*), geomep(1,*), line=1, color=255, thick=1.0 -;oplot, geomel(0,*), geomel(1,*), line=1, color=255, thick=1.0 -oplot, pdistnz(2,*)*1.0d-3, pdistnz(5,*), line=3, color=255, thick=1.0 -oplot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=1, color=255, thick=1.0 -oplot, sdist(0,*), sdist(5,*), line=1, color=255, thick=1.0 -oplot, ph1(0,*)*Dfac*1.0d-3, ph1(1,*), psym=1, color=255, thick=1.0 -;oplot, ph2(0,*), ph2(1,*), psym=4, color=255, thick=1.0 -;oplot, ph3(0,*), ph3(1,*), psym=5, color=255, thick=1.0 - -;draw box around the main surface display & fill -displaysize = floor(0.45*presresx) ; size of displayed surface -surfxpos = floor(0.520*presresx) -surfypos = floor(0.12*presresy) -xbox = [surfxpos - 1,surfxpos - 1,surfxpos + displaysize,surfxpos + displaysize,surfxpos - 1] -ybox = [surfypos - 1,surfypos + displaysize,surfypos + displaysize,surfypos - 1,surfypos - 1] -plotS, xbox, ybox, /device, color=255 - -mapscaled = congrid(map,3,displaysize,displaysize,/interp) - -;mapscaled=congrid(map,displaysize,displaysize) ; scale image -tv, mapscaled, surfxpos, surfypos, xsize=displaysize, ysize=displaysize, true=1, /device - -xyouts, 0.310*presresx, 0.04*presresy, time + fcuryear + timeunit, color=255, charsize=2*presresy/720., /device - -snapshot = TVRD(True=1) -Set_Plot, thisDevice -fnum = string(ncount,format='(I6.6)') -fname = 'presentation/presentation' + fnum + '.jpg' -Write_JPEG, fname, snapshot, True=1, Quality=90 - - -end diff --git a/examples/quasimc-global/ctem_image_regolith.pro b/examples/quasimc-global/ctem_image_regolith.pro deleted file mode 100755 index 8d7f07a7..00000000 --- a/examples/quasimc-global/ctem_image_regolith.pro +++ /dev/null @@ -1,33 +0,0 @@ -pro ctem_image_regolith,ncount,gridsize,pix,regolith,regolith_image -; Generates the regolith depth map image and saves it as a jpeg output image in the 'rego' directory -; outputs dregolith, which may be scaled to use as a console image -Compile_Opt DEFINT32 -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[gridsize,gridsize],Set_Pixel_Depth=24, Decomposed=0 -loadct, 39 -TVLCT, red, green, blue, /GET - -minref = pix * 1.0d-4 -regolith_scaled = dblarr(gridsize,gridsize) -maxreg = max(regolith) -minreg = min(regolith) -if minreg lt minref then minreg = minref -if maxreg lt minref then maxreg = minref + 1.0d30 -regolith_scaled = regolith > minreg -regolith_scaled = 254.0d0 * ((alog(regolith_scaled) - alog(minreg)) / (alog(maxreg) - alog(minreg))) - -; save regolith display -tv, regolith_scaled, 0, 0, xsize=gridsize, ysize=gridsize, /device -regolith_image = TVRD(True=1) -Set_Plot, thisDevice - -if (file_test('rego',/DIRECTORY) eq 0) then begin - file_mkdir,'rego' -endif -fnum = string(ncount,format='(I6.6)') -fname = 'rego/rego' + fnum + '.jpg' -write_jpeg, fname, regolith_image, true=1, quality=100 - -end diff --git a/examples/quasimc-global/ctem_image_shaded_relief.pro b/examples/quasimc-global/ctem_image_shaded_relief.pro deleted file mode 100755 index 73986f84..00000000 --- a/examples/quasimc-global/ctem_image_shaded_relief.pro +++ /dev/null @@ -1,52 +0,0 @@ -pro ctem_image_shaded_relief,ncount,gridsize,pix,surface_dem,height_vals,minh,maxh,dirname,shaded_image -; Generates the shaded depth map image and saves it as a jpeg output image in the 'rego' directory -; outputs dshaded, which may be scaled to use as a console image -; Uses the array height_vals for the color -Compile_Opt DEFINT32 -thisDevice = !D.Name -Set_Plot, 'Z' -Erase -Device, Set_Resolution=[gridsize,gridsize],Set_Pixel_Depth=24, Decomposed=0 -loadct, 33 -TVLCT, red, green, blue, /GET - -light=[[1,1,1],[0,0,0],[-1,-1,-1]] - -; convolution of the dem with the 3x3 matrix - -shaded=bytscl(convol(float(surface_dem),float(light))) -if max(shaded) eq 0 then shaded=255 - -; scale the dem to the height -if (minh gt maxh) then begin - tmp = minh - minh = maxh - maxh = tmp -endif - -demscaled = ((height_vals - minh) > 0.0) < (maxh-minh) -if ((maxh - minh) eq 0.0) then begin - demscaled = demscaled * 0.0 -endif else begin - demscaled = demscaled/(maxh-minh)*255.0 -endelse - -tv, demscaled, 0, 0, xsize=gridsize, ysize=gridsize, /device -shaded_image = TVRD(True=1) -shadedscl =float(shaded)/255.0 -shaded_imagearr = dblarr(3,gridsize,gridsize) -shaded_imagearr(0,*,*) = float(shaded_image(0,*,*)) * shadedscl(*,*) -shaded_imagearr(1,*,*) = float(shaded_image(1,*,*)) * shadedscl(*,*) -shaded_imagearr(2,*,*) = float(shaded_image(2,*,*)) * shadedscl(*,*) -shaded_image=round(shaded_imagearr) -Set_Plot, thisDevice - -if (file_test(dirname,/DIRECTORY) eq 0) then begin - file_mkdir,dirname -endif - -fnum = string(ncount,format='(I6.6)') -fname = dirname + '/' + dirname + fnum + '.jpg' -write_jpeg,fname,shaded_image,true=1,quality=90 - -end diff --git a/examples/quasimc-global/ctem_io_read_input.pro b/examples/quasimc-global/ctem_io_read_input.pro deleted file mode 100755 index ea135668..00000000 --- a/examples/quasimc-global/ctem_io_read_input.pro +++ /dev/null @@ -1,131 +0,0 @@ -pro ctem_io_read_input,infilename,interval,numintervals,gridsize,pix,seed,numlayers,sfdfile,impfile,maxcrat,ph1,shadedmaxhdefault,shadedminhdefault,shadedminh,shadedmaxh,restart,runtype,popupconsole,saveshaded,saverego,savepres,savetruelist -Compile_Opt DEFINT32 -print, 'Reading input file' -openr,infile,infilename, /GET_LUN -line="" -comment="!" -interval = 0.d0 -numintervals = 0 -pix=-1.0d0 -gridsize=-1 -seed = 0 -maxcrat = 1.0d0 -shadedmaxhdefault = 1 -shadedminhdefault = 1 -shadedminh = 0.d0 -shademaxnh = 0.d0 - - -; Set required strings to unset value -notset="-----NOTSET----" -sfdfile = notset -impfile = notset -sfdcompare = notset -restart = notset -runtype = notset -popupconsole = notset -saveshaded = notset -saverego = notset -savepres = notset -savetruelist = notset -while (not EOF(infile)) do begin - readf,infile,line - if (~strcmp(line,comment,1)) then begin - substrings = strsplit(line,' ',/extract) - if strmatch(substrings(0),'pix',/fold_case) then reads,substrings(1),pix - if strmatch(substrings(0),'gridsize',/fold_case) then reads,substrings(1),gridsize - if strmatch(substrings(0),'seed',/fold_case) then reads,substrings(1),seed - if strmatch(substrings(0),'sfdfile',/fold_case) then reads,substrings(1),sfdfile - if strmatch(substrings(0),'impfile',/fold_case) then reads,substrings(1),impfile - if strmatch(substrings(0),'maxcrat',/fold_case) then reads,substrings(1),maxcrat - if strmatch(substrings(0),'sfdcompare',/fold_case) then reads,substrings(1),sfdcompare - if strmatch(substrings(0),'interval',/fold_case) then reads,substrings(1),interval - if strmatch(substrings(0),'numintervals',/fold_case) then reads,substrings(1),numintervals - if strmatch(substrings(0),'popupconsole',/fold_case) then reads,substrings(1),popupconsole - if strmatch(substrings(0),'saveshaded',/fold_case) then reads,substrings(1),saveshaded - if strmatch(substrings(0),'saverego',/fold_case) then reads,substrings(1),saverego - if strmatch(substrings(0),'savepres',/fold_case) then reads,substrings(1),savepres - if strmatch(substrings(0),'savetruelist',/fold_case) then reads,substrings(1),savetruelist - if strmatch(substrings(0),'runtype',/fold_case) then reads,substrings(1),runtype - if strmatch(substrings(0),'restart',/fold_case) then reads,substrings(1),restart - if strmatch(substrings(0),'shadedminh',/fold_case) then begin - reads,substrings(1),shadedminh - shadedminhdefault = 0 - endif - if strmatch(substrings(0),'shadedmaxh',/fold_case) then begin - reads,substrings(1),shadedmaxh - shadedmaxhdefault = 0 - endif - end -end -if interval le 0.0d0 then begin - print,'Invalid value for or missing variable INTERVAL in ' + infilename - stop -end -if numintervals le 0 then begin - print,'Invalid value for or missing variable NUMINTERVALS in ' + infilename - stop -end -if pix le 0.0d0 then begin - print,'Invalid value for or missing variable PIX in ' + infilename - stop -end -if gridsize le 0 then begin - print,'Invalid value for or missing variable GRIDSIZE in ' + infilename - stop -end -if seed eq 0 then begin - print,'Invalid value for or missing variable SEED in ' + infilename - stop -end -if strmatch(sfdfile,notset) then begin - print,'Invalid value for or missing variable SFDFILE in ' + infilename - stop -end -if strmatch(impfile,notset) then begin - print,'Invalid value for or missing variable IMPFILE in ' + infilename - stop -end -if strmatch(popupconsole,notset) then begin - print,'Invalid value for or missing variable POPUPCONSOLE in ' + infilename - stop -end -if strmatch(saveshaded,notset) then begin - print,'Invalid value for or missing variable SAVESHADED in ' + infilename - stop -end -if strmatch(saverego,notset) then begin - print,'Invalid value for or missing variable SAVEREGO in ' + infilename - stop -end -if strmatch(savepres,notset) then begin - print,'Invalid value for or missing variable SAVEPRES in ' + infilename - stop -end -if strmatch(savetruelist,notset) then begin - print,'Invalid value for or missing variable SAVETRUELIST in ' + infilename - stop -end -if strmatch(runtype,notset) then begin - print,'Invalid value for or missing variable RUNTYPE in ' + infilename - stop -end -if strmatch(restart,notset) then begin - print,'Invalid value for or missing variable RESTART in ' + infilename - stop -end - -free_lun,infile - -ph1 = dblarr(3,1) -if ~strmatch(sfdcompare,notset) then begin - cnum = file_lines(sfdcompare) - ph1 = dblarr(3,cnum) - openr,COMP,sfdcompare,/GET_LUN - readf,COMP,ph1 - close,COMP - free_lun,COMP -end -free_lun,infile - -end diff --git a/examples/quasimc-global/ctem_io_read_old.pro b/examples/quasimc-global/ctem_io_read_old.pro deleted file mode 100755 index a9013687..00000000 --- a/examples/quasimc-global/ctem_io_read_old.pro +++ /dev/null @@ -1,45 +0,0 @@ -pro ctem_io_read_old,gridsize,surface_dem,regolith,odist,tdist,pdist,mass -Compile_Opt DEFINT32 -openr,LUN,'surface_ejc.dat',/GET_LUN -readu,LUN,regolith -free_lun,LUN - -openr,LUN,'surface_dem.dat',/GET_LUN -readu,LUN,surface_dem -free_lun,LUN - -distl = file_lines('odistribution.dat') - 1 -pdistl = file_lines('pdistribution.dat') - 1 - -;read in observed cumulative distribution -odist = dblarr(6,distl) -line = "temp" -openr,LUN,'odistribution.dat',/GET_LUN -readf,LUN,line -readf,LUN,odist -close,LUN -free_lun,LUN - -;read in true cumulative distribution -tdist = dblarr(6,distl) -openr,LUN,'tdistribution.dat',/GET_LUN -readf,LUN,line -readf,LUN,tdist -close,LUN -free_lun,LUN - -;read in production function -pdist = dblarr(6,pdistl) -openr,LUN,'pdistribution.dat',/GET_LUN -readf,LUN,line -readf,LUN,pdist -close,LUN -free_lun,LUN - -;read in accumulated mass from file -openr,LUN,'impactmass.dat',/GET_LUN -readf,LUN,mass -close,LUN -free_lun,LUN - -end diff --git a/examples/quasimc-global/ctem_io_readers.py b/examples/quasimc-global/ctem_io_readers.py deleted file mode 100644 index 0881ad1a..00000000 --- a/examples/quasimc-global/ctem_io_readers.py +++ /dev/null @@ -1,146 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model file reading utilities -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 -# -#Known issues for operation with CTEM -#1) ctem.in has 1.0d0 value for maxcrat which is not readable by Python - -import numpy - -def read_ctemin(parameters,notset): - #Read and parse ctem.in file - inputfile = parameters['workingdir'] + parameters['ctemfile'] - - #Read ctem.in file - print('Reading input file '+ parameters['ctemfile']) - fp = open(inputfile,'r') - lines = fp.readlines() - fp.close() - - #Process file text - for line in lines: - fields = line.split() - if len(fields) > 0: - if ('pix' == fields[0].lower()): parameters['pix']=float(fields[1]) - if ('gridsize' == fields[0].lower()): parameters['gridsize']=int(fields[1]) - if ('seed' == fields[0].lower()): parameters['seed']=int(fields[1]) - if ('sfdfile' == fields[0].lower()): parameters['sfdfile']=fields[1] - if ('impfile' == fields[0].lower()): parameters['impfile']=fields[1] - if ('maxcrat' == fields[0].lower()): parameters['maxcrat']=float(fields[1]) - if ('sfdcompare' == fields[0].lower()): parameters['sfdcompare']=fields[1] - if ('interval' == fields[0].lower()): parameters['interval']=float(fields[1]) - if ('numintervals' == fields[0].lower()): parameters['numintervals']=int(fields[1]) - if ('popupconsole' == fields[0].lower()): parameters['popupconsole']=fields[1] - if ('saveshaded' == fields[0].lower()): parameters['saveshaded']=fields[1] - if ('saverego' == fields[0].lower()): parameters['saverego']=fields[1] - if ('savepres' == fields[0].lower()): parameters['savepres']=fields[1] - if ('savetruelist' == fields[0].lower()): parameters['savetruelist']=fields[1] - if ('runtype' == fields[0].lower()): parameters['runtype']=fields[1] - if ('restart' == fields[0].lower()): parameters['restart']=fields[1] - if ('quasimc' == fields[0].lower()): parameters['quasimc']=fields[1] - if ('realcraterlist' == fields[0].lower()): parameters['realcraterlist']=fields[1] - if ('shadedminh' == fields[0].lower()): - parameters['shadedminh'] = float(fields[1]) - parameters['shadedminhdefault'] = 0 - if ('shadedmaxh' == fields[0].lower()): - parameters['shadedmaxh'] = float(fields[1]) - parameters['shadedmaxhdefault'] = 0 - - #Test values for further processing - if (parameters['interval'] <= 0.0): - print('Invalid value for or missing variable INTERVAL in '+ inputfile) - if (parameters['numintervals'] <= 0): - print('Invalid value for or missing variable NUMINTERVALS in '+ inputfile) - if (parameters['pix'] <= 0.0): - print('Invalid value for or missing variable PIX in '+ inputfile) - if (parameters['gridsize'] <= 0): - print('Invalid value for or missing variable GRIDSIZE in '+ inputfile) - if (parameters['seed'] == 0): - print('Invalid value for or missing variable SEED in '+ inputfile) - if (parameters['sfdfile'] == notset): - print('Invalid value for or missing variable SFDFILE in '+ inputfile) - if (parameters['impfile'] == notset): - print('Invalid value for or missing variable IMPFILE in '+ inputfile) - if (parameters['popupconsole'] == notset): - print('Invalid value for or missing variable POPUPCONSOLE in '+ inputfile) - if (parameters['saveshaded'] == notset): - print('Invalid value for or missing variable SAVESHADED in '+ inputfile) - if (parameters['saverego'] == notset): - print('Invalid value for or missing variable SAVEREGO in '+ inputfile) - if (parameters['savepres'] == notset): - print('Invalid value for or missing variable SAVEPRES in '+ inputfile) - if (parameters['savetruelist'] == notset): - print('Invalid value for or missing variable SAVETRUELIST in '+ inputfile) - if (parameters['runtype'] == notset): - print('Invalid value for or missing variable RUNTYPE in '+ inputfile) - if (parameters['restart'] == notset): - print('Invalid value for or missing variable RESTART in '+ inputfile) - - return - -def read_formatted_ascii(filename, skip_lines): - #Generalized ascii text reader - #For use with sfdcompare, production, odist, tdist, pdist data files - data = numpy.genfromtxt(filename, skip_header = skip_lines) - return data - -def read_unformatted_binary(filename, gridsize): - #Read unformatted binary files created by Fortran - #For use with surface ejecta and surface dem data files - dt = numpy.float - data = numpy.fromfile(filename, dtype = dt) - data.shape = (gridsize,gridsize) - - return data - -def read_ctemdat(parameters, seedarr): - #Read and parse ctem.dat file - datfile = parameters['workingdir'] + 'ctem.dat' - - #Read ctem.dat file - print('Reading input file '+ parameters['datfile']) - fp = open(datfile,'r') - lines = fp.readlines() - fp.close() - - #Parse file lines and update parameter fields - fields = lines[0].split() - if len(fields) > 0: - parameters['totalimpacts'] = float(fields[0]) - parameters['ncount'] = int(fields[1]) - parameters['curyear'] = float(fields[2]) - parameters['restart'] = fields[3] - parameters['fracdone'] = float(fields[4]) - parameters['masstot'] = float(fields[5]) - - #Parse remainder of file to build seed array - nlines = len(lines) - index = 1 - while (index < nlines): - fields = lines[index].split() - seedarr[index - 1] = float(fields[0]) - index += 1 - - parameters['seedn'] = index - 1 - - return - -def read_impact_mass(filename): - #Read impact mass file - - fp=open(filename,'r') - line=fp.readlines() - fp.close() - - fields = line[0].split() - if (len(fields) > 0): - mass = float(fields[0]) - else: - mass = 0 - - return mass \ No newline at end of file diff --git a/examples/quasimc-global/ctem_io_writers.py b/examples/quasimc-global/ctem_io_writers.py deleted file mode 100644 index a5e72b26..00000000 --- a/examples/quasimc-global/ctem_io_writers.py +++ /dev/null @@ -1,282 +0,0 @@ -#!/usr/local/bin/python -# -#Cratered Terrain Evolution Model file and graphical writing utilities -# -#Original IDL design: Jim Richardson, Arecibo Observatory -#Revised IDL design: David Minton, Purdue University -#Re-engineered design and Python implementation: Matthew Route, Purdue University -#August 2016 -# - -import matplotlib -from matplotlib import pyplot -import numpy -import os -import shutil -import scipy -from scipy import signal - -#Set pixel scaling common for image writing, at 1 pixel/ array element -dpi = 72.0 - -#Write production function to file production.dat -#This file format does not exactly match that generated from IDL. Does it work? -def write_production(parameters, production): - filename = parameters['workingdir'] + parameters['sfdfile'] - numpy.savetxt(filename, production, fmt='%1.8e', delimiter=' ') - - return - -def create_dir_structure(parameters): - #Create directories for various output files if they do not already exist - directories=['dist','misc','rego','rplot','surf','shaded'] - - for directory in directories: - dir_test = parameters['workingdir'] + directory - if not os.path.isdir(dir_test): - os.makedirs(dir_test) - - return - -def write_ctemdat(parameters, seedarr): - #Write various parameters and random number seeds into ctem.dat file - filename = parameters['workingdir'] + parameters['datfile'] - fp = open(filename,'w') - - template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n" - fp.write(template % parameters) - - #Write random number seeds to the file - for index in range(parameters['seedn']): - fp.write("%12d\n" % seedarr[index]) - - fp.close() - - return - -def copy_dists(parameters): - #Save copies of distribution files - - orig_list = ['odistribution', 'ocumulative', 'pdistribution', 'tdistribution'] - dest_list = ['odist', 'ocum', 'pdist', 'tdist'] - - for index in range(len(orig_list)): - forig = parameters['workingdir'] + orig_list[index] + '.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + dest_list[index] + "_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - forig = parameters['workingdir'] + 'impactmass.dat' - fdest = parameters['workingdir'] + 'misc' + os.sep + "mass_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - if (parameters['savetruelist'].upper() == 'T'): - forig = parameters['workingdir'] + 'tcumulative.dat' - fdest = parameters['workingdir'] + 'dist' + os.sep + "tcum_%06d.dat" % parameters['ncount'] - shutil.copy2(forig, fdest) - - return - -#Possible references -#http://nbviewer.jupyter.org/github/ThomasLecocq/geophysique.be/blob/master/2014-02-25%20Shaded%20Relief%20Map%20in%20Python.ipynb - -def image_dem(parameters, surface_dem): - - #Create surface dem map - solar_angle = 20.0 - dem_map = numpy.copy(surface_dem) - numpy.roll(surface_dem, 1, 0) - dem_map = (0.5 * numpy.pi) + numpy.arctan2(dem_map, parameters['pix']) - dem_map = dem_map - numpy.radians(solar_angle) * (0.5 * numpy.pi) - numpy.place(dem_map, dem_map > (0.5 * numpy.pi), 0.5 *numpy.pi) - dem_map = numpy.absolute(dem_map) - dem_map = 254.0 * numpy.cos(dem_map) - - #Save image to file - filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - fig.figimage(dem_map, cmap = matplotlib.cm.gray, origin = 'lower') - matplotlib.pyplot.savefig(filename) - - return - -def image_regolith(parameters, regolith): - - #Create scaled regolith image - minref = parameters['pix'] * 1.0e-4 - maxreg = numpy.amax(regolith) - minreg = numpy.amin(regolith) - if (minreg < minref): minreg = minref - if (maxreg < minref): maxreg = (minref + 1.0e3) - regolith_scaled = numpy.copy(regolith) - numpy.place(regolith_scaled, regolith_scaled < minref, minref) - regolith_scaled = 254.0 * ((numpy.log(regolith_scaled) - numpy.log(minreg)) / (numpy.log(maxreg) - numpy.log(minreg))) - - #Save image to file - filename = parameters['workingdir'] + 'rego' + os.sep + "rego%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - fig.figimage(regolith_scaled, cmap = matplotlib.cm.nipy_spectral, origin = 'lower') - matplotlib.pyplot.savefig(filename) - - return - -def image_shaded_relief(parameters, surface_dem): - #The color scale and appearance of this do not quite match the IDL version - - #Create image by convolving DEM with 3x3 illumination matrix - light = numpy.array([[1.0, 1.0, 1.0], [0.0, 0.0, 0.0], [-1.0, -1.0, -1.0]]) - convolved_map = scipy.signal.convolve2d(surface_dem, light, mode = 'same') - - #Adjust output to resemble IDL (north slopes illuminated, south in shadow) - convolved_map = convolved_map * -1.0 - convolved_map[0,:]=0.0 - convolved_map[-1,:] =0.0 - convolved_map[:,0]=0.0 - convolved_map[:,-1]=0.0 - - #If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data - if (parameters['shadedminhdefault'] == 1): shadedminh = numpy.amin(surface_dem) - if (parameters['shadedmaxhdefault'] == 1): shadedmaxh = numpy.amax(surface_dem) - - #If min and max appear to be reversed, then fix them - if (parameters['shadedminh'] > parameters['shadedmaxh']): - temp = parameters['shadedminh'] - shadedminh = parameters['shadedmaxh'] - shadedmaxh = temp - else: - shadedminh = parameters['shadedminh'] - shadedmaxh = parameters['shadedmaxh'] - - #If dynamic range is valid, construct a shaded DEM - dynamic_range = shadedmaxh - shadedminh - if (dynamic_range != 0): - dem_scaled = numpy.copy(surface_dem) - shadedminh - numpy.place(dem_scaled, dem_scaled < 0.0, 0.0) - numpy.place(dem_scaled, dem_scaled > dynamic_range, dynamic_range) - dem_scaled = dem_scaled / dynamic_range - else: - dem_scaled = numpy.copy(surface_dem) * 0.0 - - #Generate shaded depth map with surface_dem color scaling (RGBA) - shaded = scipy.misc.bytescale(convolved_map) - if numpy.amax(shaded) == 0: shaded=255 - shadedscl = shaded / 255.0 - - #shaded_imagearr = matplotlib.cm.jet(dem_scaled) - #print dem_scaled[0:4,0:4] - #shaded_imagearr[:,:,0] = shaded_imagearr[:,:,0] * shadedscl - #shaded_imagearr[:,:,1] = shaded_imagearr[:,:,1] * shadedscl - #shaded_imagearr[:,:,2] = shaded_imagearr[:,:,2] * shadedscl - #shaded_imagearr[:,:,3] = shaded_imagearr[:,:,3] * shadedscl - - #Delivers nearly proper coloring, but no shaded relief - shaded_imagearr = dem_scaled * shadedscl - shaded_imagearr = matplotlib.cm.jet(shaded_imagearr) - shaded_imagearr = numpy.around(shaded_imagearr, decimals = 1) - - #Save image to file - filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - fig.figimage(shaded_imagearr, cmap = matplotlib.cm.jet, origin = 'lower') - matplotlib.pyplot.savefig(filename) - return - -def create_rplot(parameters,odist,pdist,tdist,ph1): - #Parameters: empirical saturation limit and dfrac - satlimit = 3.12636 - dfrac = 2**(1./4) * 1.0e-3 - - #Calculate geometric saturation - minx = (parameters['pix'] / 3.0) * 1.0e-3 - maxx = 3 * parameters['pix'] * parameters['gridsize'] * 1.0e-3 - geomem = numpy.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]]) - geomep = numpy.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]]) - - #Create distribution arrays without zeros for plotting on log scale - idx = numpy.nonzero(odist[:,5]) - odistnz = odist[idx] - odistnz = odistnz[:,[2,5]] - - idx = numpy.nonzero(pdist[:,5]) - pdistnz = pdist[idx] - pdistnz = pdistnz[:,[2,5]] - - idx = numpy.nonzero(tdist[:,5]) - tdistnz = tdist[idx] - tdistnz = tdistnz[:,[2,5]] - - #Correct pdist - pdistnz[:,1] = pdistnz[:,1] * parameters['curyear'] / parameters['interval'] - - #Create sdist bin factors, which contain one crater per bin - area = (parameters['gridsize'] * parameters['pix'] * 1.0e-3)**2. - plo = 1 - sq2 = 2**(1./2) - while (sq2**plo > minx): - plo = plo - 1 - phi = plo + 1 - while (sq2**phi < maxx): - phi = phi + 1 - n = phi - plo + 1 - sdist = numpy.zeros([n , 2]) - p = plo - for index in range(n): - sdist[index, 0] = sq2**p - sdist[index, 1] = sq2**(2.0*p + 1.5) / (area * (sq2 - 1)) - p = p + 1 - - #Create time label - tlabel = "%5.4e" % parameters['curyear'] - tlabel = tlabel.split('e') - texp = str(int(tlabel[1])) - timelabel = 'Time = '+ r'${}$ x 10$^{}$'.format(tlabel[0], texp) + ' yrs' - - #Save image to file - filename = parameters['workingdir'] + 'rplot' + os.sep + "rplot%06d.png" % parameters['ncount'] - height = parameters['gridsize'] / dpi - width = height - fig = matplotlib.pyplot.figure(figsize = (width, height), dpi = dpi) - - #Alter background color to be black, and change axis colors accordingly - matplotlib.pyplot.style.use('dark_background') - matplotlib.pyplot.rcParams['axes.prop_cycle'] - - #Plot data - matplotlib.pyplot.plot(odistnz[:,0]*1.0e-3, odistnz[:,1], linewidth=3.0, color = 'blue') - matplotlib.pyplot.plot(pdistnz[:,0]*1.0e-3, pdistnz[:,1], linewidth=2.0, linestyle='dashdot', color = 'white') - matplotlib.pyplot.plot(tdistnz[:,0]*1.0e-3, tdistnz[:,1], linewidth=2.0, color = 'red') - - matplotlib.pyplot.plot(geomem[:,0], geomem[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - matplotlib.pyplot.plot(geomep[:,0], geomep[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - matplotlib.pyplot.plot(sdist[:,0], sdist[:,1], linewidth=2.0, linestyle =':', color = 'yellow') - - matplotlib.pyplot.plot(ph1[:,0] * dfrac, ph1[:,1], 'wo') - - #Create plot labels - matplotlib.pyplot.title('Crater Distribution R-Plot',fontsize=22) - matplotlib.pyplot.xlim([minx, maxx]) - matplotlib.pyplot.xscale('log') - matplotlib.pyplot.xlabel('Crater Diameter (km)',fontsize=18) - matplotlib.pyplot.ylim([5.0e-4, 5.0]) - matplotlib.pyplot.yscale('log') - matplotlib.pyplot.ylabel('R Value', fontsize=18) - matplotlib.pyplot.text(1.0e-2, 1.0, timelabel, fontsize=18) - - matplotlib.pyplot.tick_params(axis='both', which='major', labelsize=14) - matplotlib.pyplot.tick_params(axis='both', which='minor', labelsize=12) - - matplotlib.pyplot.savefig(filename) - - return - -def write_realcraters(parameters, realcraters): - """writes file of real craters for use in quasi-MC runs""" - - filename = parameters['workingdir'] + 'craterlist.dat' - numpy.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t') - - return \ No newline at end of file diff --git a/examples/quasimc-global/ctem_window_display.pro b/examples/quasimc-global/ctem_window_display.pro deleted file mode 100755 index d6edf5b9..00000000 --- a/examples/quasimc-global/ctem_window_display.pro +++ /dev/null @@ -1,249 +0,0 @@ -pro ctem_window_display,ncount,totalimpacts,gridsize,pix,curyear,masstot,odist,pdist,tdist,ph1,surface_dem,regolith,map,popupconsole -; Produces the console window display. The array 'map' is used to generate the image on the right-hand side -Compile_Opt DEFINT32 - -; console output graphics resolution -minref = pix * 1.0d-4 -conresx = 1280 -;conresy = 0.758*conresx -profsize = 0.00 -conresy = (0.60 + profsize) * conresx -sun = 20.d0 ; sun angle for display - -minx = (pix / 3.0d0) * 1d-3 -maxx = 3 * pix * gridsize * 1d-3 -charfac = 1.6 ; character size multiplier - -if strmatch(popupconsole,'T',/fold_case) then begin - device, decomposed=0, retain=2 - thisDevice = !D.NAME - SET_PLOT, 'Z' - TVLCT, red, green, blue, /GET - SET_PLOT, thisDevice - !P.MULTI = [0, 1, 2] - Window, 0, xsize=conresx, ysize=conresy -endif else begin - SET_PLOT, 'Z', /COPY - device, set_resolution=[conresx,conresy],decomposed=0, Z_Buffer=0 - TVLCT, red, green, blue, /GET - !P.MULTI = [0, 1, 2] - erase -endelse - -;set up first color scale -loadct, 0 - -;strings for display -time = 'Time = ' -timeunit = ' yr' -timp = 'Total impacts (dotdash) = ' -tcrt = 'Total craters (line) = ' -ocrt = 'Countable craters (bold) = ' -epwr = 'Mean regolith depth = ' -mtot = 'Impacted mass = ' -c1lab = 'Min. Elevation' -c2lab = 'Mean Elevation' -c3lab = 'Max. Elevation' - -maxelev = max(surface_dem) -minelev = min(surface_dem) -medelev = mean(surface_dem) -c1 = string(minelev, format = '(F9.1)') + ' m' -c2 = string(medelev, format = '(F9.1)') + ' m' -c3 = string(maxelev, format = '(F9.1)') + ' m' -tslp = mean(regolith) - -;geometric saturation -geom = dblarr(2,2) -geomem = dblarr(2,2) -geomep = dblarr(2,2) -geomel = dblarr(2,2) -geom(0,0) = minx -geom(0,1) = maxx -geom(1,*) = 3.12636d0 -geomem(0,0) = minx -geomem(0,1) = maxx -geomem(1,*) = 0.156318d0 -geomep(0,0) = minx -geomep(0,1) = maxx -geomep(1,*) = 0.312636d0 -geomel(0,0) = minx -geomel(0,1) = maxx -geomel(1,*) = 0.0312636d0 - -; Remove zeros -nz = where(odist(5,*) ne 0.0,count) -if (count gt 0) then begin - odistnz = dblarr(6,count) - odistnz(*,*)= odist(*,nz) -endif else begin - odistnz = odist -endelse - -nz = where(tdist(5,*) ne 0.0,count) -if (count gt 0) then begin - tdistnz = dblarr(6,count) - tdistnz(*,*)= tdist(*,nz) -endif else begin - tdistnz = tdist -endelse - -nz = where(pdist(5,*) ne 0.0,count) -if (count gt 0) then begin - pdistnz = dblarr(6,count) - pdistnz(*,*)= pdist(*,nz) - -endif else begin - pdistnz = pdist -endelse - -; create r-plot array containing exactly 1 crater per bin -area = (gridsize * pix * 1d-3)^2 -plo = 1 -while (sqrt(2.0d0)^plo gt minx) do begin - plo = plo - 1 -endwhile -phi = plo + 1 -while (sqrt(2.0d0)^phi lt maxx) do begin - phi = phi + 1 -endwhile -n = phi - plo -sdist = dblarr(6,n + 1) -p = plo -for i=0,n do begin - sdist(0,i) = sqrt(2.d0)^p - sdist(1,i) = sqrt(2.d0)^(p+1) - sdist(2,i) = sqrt(sdist(0,i) * sdist(1,i)) - sdist(3,i) = 1.0d0 - sdist(5,i) = (sdist(2,i))^3 / (area * (sdist(1,i) - sdist(0,i))) - p = p + 1 -endfor -sdist(4,*) = reverse(total(sdist(3,*),/cumulative),2) - -plot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=0, color=255, thick=2.0, $ - TITLE='Crater Distribution R-Plot', charsize=charfac*conresy/720, $ - XTITLE='Crater Diameter (km)', $ - XRANGE=[minx,maxx], /XLOG, XSTYLE=1, $ - YTITLE='R Value', $ - YRANGE=[5.0e-4,5.0e0], /YLOG, YSTYLE=1, $ - /device, pos = [0.085*conresx,(0.25 + profsize)*conresy,0.44*conresx,0.85*conresy] - -Dfac = sqrt(sqrt(2.0d0)) - -;display observed crater distribution -oplot, tdistnz(2,*)*1.0d-3, tdistnz(5,*), line=0, color=255, thick=1.0 -oplot, geom(0,*), geom(1,*), line=2, color=255, thick=1.0 -oplot, geomem(0,*), geomem(1,*), line=1, color=255, thick=1.0 -oplot, geomep(0,*), geomep(1,*), line=1, color=255, thick=1.0 -oplot, geomel(0,*), geomel(1,*), line=1, color=255, thick=1.0 -oplot, pdistnz(2,*)*1.0d-3, pdistnz(5,*), line=3, color=255, thick=1.0 -oplot, odistnz(2,*)*1.0d-3, odistnz(5,*), line=1, color=255, thick=1.0 -oplot, sdist(0,*), sdist(5,*), line=1, color=255, thick=1.0 -oplot, ph1(0,*)*Dfac*1.0d-3, ph1(1,*), psym=1, color=255, thick=1.0 -;oplot, ph2(0,*), ph2(1,*), psym=4, color=255, thick=1.0 -;oplot, ph3(0,*), ph3(1,*), psym=5, color=255, thick=1.0 - -; set up profile display -surfpos = dblarr(gridsize) -surfpro = dblarr(gridsize) -surfhpro = dblarr(gridsize) -surfrpro = dblarr(gridsize) -for k = 0,gridsize-1 do begin - surfpos(k) = (-0.5d0*(gridsize-1) + double(k)) * pix -endfor -surfpro(*) = surface_dem(*,gridsize/2-1) -surfhpro(*) = regolith(*,gridsize/2-1) -surfrpro(*) = surfpro(*) - surfhpro(*) - -; set up profile -;maxelevp = max(surfpro) -;minelevp = min(surfpro) -;medelevp = mean(surfpro) -;vertrange = 1.5 * abs(maxelevp-minelevp) -;vertadd = 0.5d0 * (vertrange - abs(maxelevp-minelevp)) -;horzrange = vertrange * 5.86666666667d0 -;if (horzrange gt (pix*gridsize)) then horzrange = pix*gridsize -;plot, surfpos(*)/1.0d3, surfpro(*)/1.0d3, line=0, color=255, thick=1.0, $ -; TITLE='Matrix Cross-Section', charsize=1.0, $ -; XTITLE='Horizontal Position (km)', $ -; XRANGE=[(-1.0d0*horzrange/2.0d3),(horzrange/2.0d3)], XSTYLE=1, $ -; YTITLE='Elevation (km)', $ -; YRANGE=[((minelevp - vertadd)/1.0d3),((maxelevp + vertadd)/1.0d3)], YSTYLE=1, $ -; /device, pos = [0.063*conresx,0.049*conresy,0.99*conresx,0.26*conresy] -;oplot, surfpos(*)*1.0d-3, surfrpro(*)*1.0d-3, line=0, color=255, thick=1.0 - -;display text -dum = string(format='(E10.4)', curyear) -parts = strsplit(dum, 'E', /extract) -fcuryear = strcompress(string(format='(F6.4,"x10!U",i,"!N")', parts[0], parts[1])) + ' yr' -ftimp = string(totalimpacts, format = '(I13)') -ftcnt = string(tdist(4,0), format = '(I13)') -focnt = string(odist(4,0), format = '(I13)') -ftslp = string(tslp, format = '(F13.3)') + ' m' - - -dum = string(format='(E10.4)',masstot) -parts = strsplit(dum, 'E', /extract) -fmtot = strcompress(string(format='(F6.4,"x10!U",i,"!N")', parts[0], parts[1])) + ' kg' - -; Display information text below the R-plot -texttop = 0.15 + profsize -textspc = 0.035 -textleft = 0.05*conresx -textbreak = 0.25*conresx -xyouts, textleft, (texttop - 0*textspc)*conresy, time, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 1*textspc)*conresy, timp, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 2*textspc)*conresy, tcrt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 3*textspc)*conresy, ocrt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft, (texttop - 4*textspc)*conresy, mtot, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 0*textspc)*conresy, fcuryear, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 1*textspc)*conresy, ftimp, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 2*textspc)*conresy, ftcnt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 3*textspc)*conresy, focnt, color=255, charsize=charfac*conresy/720., /device -xyouts, textleft + textbreak, (texttop - 4*textspc)*conresy, fmtot, color=255, charsize=charfac*conresy/720., /device - -; Display information text above the R-plot -xyouts, 0.0368*conresx, 0.972*conresy, c1lab, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.195*conresx, 0.972*conresy, c2lab, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.353*conresx, 0.972*conresy, c3lab, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.0368*conresx, 0.944*conresy, c1, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.195*conresx, 0.944*conresy, c2, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.353*conresx, 0.944*conresy, c3, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.100*conresx, 0.917*conresy, epwr, color=255, charsize=charfac*conresy/720., /device -xyouts, 0.278*conresx, 0.917*conresy, ftslp, color=255, charsize=charfac*conresy/720., /device - -;print to screen -print, ncount, ' time = ', curyear -;print, ncount, ' tot impacts = ', tdist -;print, ncount, ' tot craters = ', totalimpacts -;print, ncount, ' obs craters = ', ocrcnt -;print, ncount, ' regolith = ', tslp -;print, ncount, ' prod R = ', mean(prval(1,*)) - -;draw box around the main surface display & fill -displaysize = floor(0.53*conresx) ; size of displayed surface -surfxpos = floor(0.462*conresx) -surfypos = floor((0.05 + profsize)*conresy) -xbox = [surfxpos - 1,surfxpos - 1,surfxpos + displaysize,surfxpos + displaysize,surfxpos - 1] -ybox = [surfypos - 1,surfypos + displaysize,surfypos + displaysize,surfypos - 1,surfypos - 1] -plotS, xbox, ybox, /device, color=255 - -; for display -mapscaled = congrid(map,3,displaysize,displaysize,/interp) - -tv, mapscaled, surfxpos, surfypos, xsize=displaysize, ysize=displaysize, true=1, /device - -; ---------- saving window ---------- - -; save console window -fnum = string(ncount,format='(I6.6)') -; print screen to JPEG file -if (file_test('console',/DIRECTORY) eq 0) then begin - file_mkdir,'console' -endif -print, ncount, ' saving window' -fname = 'console/console' + fnum + '.jpg' -dwin = tvrd(true=1) -write_jpeg, fname, dwin, true=1, quality=90 - -end From 65427e4513cdb02f19ff355a868f9b78be7aa917 Mon Sep 17 00:00:00 2001 From: David Minton Date: Wed, 9 Mar 2022 19:32:30 -0500 Subject: [PATCH 07/25] Flipped y axis of the DEM when passing it into the LightSource shade functions --- python/ctem/ctem/util.py | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index 5d2c1700..d77b84f9 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -46,7 +46,7 @@ def image_dem(user, DEM): solar_angle = 20.0 # user['solar_angle'] ls = LightSource(azdeg=azimuth, altdeg=solar_angle) - dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix) + dem_img = ls.hillshade(np.flip(DEM, axis=0), vert_exag=ve, dx=pix, dy=pix) # Generate image to put into an array height = gridsize / dpi @@ -116,7 +116,7 @@ def image_shaded_relief(user, DEM): else: shadedmaxh = user['shadedmaxh'] - dem_img = ls.shade(DEM, cmap=cmap, blend_mode=mode, fraction=1.0, + dem_img = ls.shade(np.flip(DEM, axis=0), cmap=cmap, blend_mode=mode, fraction=1.0, vert_exag=ve, dx=pix, dy=pix, vmin=shadedminh, vmax=shadedmaxh) From 29e4b76fe0a2a5ab5802de9e00081dd189362e55 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Wed, 13 Jul 2022 16:56:34 -0400 Subject: [PATCH 08/25] Fixed underflow bug --- src/regolith/regolith_depth_model.f90 | 9 ++++++++- 1 file changed, 8 insertions(+), 1 deletion(-) diff --git a/src/regolith/regolith_depth_model.f90 b/src/regolith/regolith_depth_model.f90 index 8c98bc80..8f37eff8 100644 --- a/src/regolith/regolith_depth_model.f90 +++ b/src/regolith/regolith_depth_model.f90 @@ -19,6 +19,7 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p) use module_globals use module_util + use, intrinsic :: ieee_arithmetic use module_regolith, EXCEPT_THIS_ONE => regolith_depth_model implicit none @@ -35,7 +36,9 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p) integer(I4B) :: i, j real(DP) :: h, dr, f, fmin, fmax real(DP) :: ntotsubcrat + logical :: underflow + if (ieee_support_underflow_control(psum)) call ieee_set_underflow_mode(gradual=.true.) ! Smallest crater size in sub-pixel crater regime (regolith scaling column in "nflux") ! Do it in terms of number of craters ntotsubcrat = 0._DP @@ -93,7 +96,11 @@ subroutine regolith_depth_model(user,domain,finterval,nflux,p) f = 0.5_DP * dr * ( fmin + fmax ) psum = psum + f end do - p(2,i) = 1.0_DP - exp(-1.0_DP * PI * psum * t) + if (abs(psum) > epsilon(psum)) then + p(2,i) = 1.0_DP - exp(-1.0_DP * PI * psum * t) + else + p(2,i) = 1.0_DP + end if end do return From 7b427f6d8d626b20926e3550f049ef0bee2b1692 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Thu, 14 Jul 2022 14:39:02 -0400 Subject: [PATCH 09/25] Added check to make sure that sin and cos terms were within the correct bounds --- src/regolith/regolith_melt_glass.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/regolith/regolith_melt_glass.f90 b/src/regolith/regolith_melt_glass.f90 index 339bcce1..59795081 100644 --- a/src/regolith/regolith_melt_glass.f90 +++ b/src/regolith/regolith_melt_glass.f90 @@ -128,7 +128,7 @@ subroutine regolith_melt_glass(user,crater,age,age_resolution,ebh,rm,eradc,lrad, newlayer%thickness = ebh ! default value: stream tube's volume = paraboloid shell's volume newlayer%comp = 0.0_DP rints = sqrt(rm**2 - (crater%imp/2.0)**2) - cosvints = eradi / (crater%imp + eradi) + cosvints = min(max(eradi / (crater%imp + eradi), -1.0_DP), 1.0_DP) sinvints = sqrt(1.0 - cosvints**2) xvints = eradi * ( 1.0 - cosvints) * sinvints volv1 = regolith_streamtube_volume_func(eradi,0.0_DP,xvints,deltar) From 3df361b3ab07a7ac62b778b8b3bb746a52cb8ec8 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 15 Jul 2022 10:57:20 -0400 Subject: [PATCH 10/25] Check if a crater does not generate topography. If not, skip it. --- src/crater/crater_populate.f90 | 6 ++++++ 1 file changed, 6 insertions(+) diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index e60b68ce..c3c779c3 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -240,6 +240,12 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! Place crater onto the surface call crater_emplace(user,surf,crater,domain,ejbmass) + if (abs(ejbmass) < tiny(ejbmass)) then ! Crater has no topography. Discard and move on. + ntrue = ntrue -1 + mass = mass - crater%impmass + cycle + end if + call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance ! For very steep size distributions, only a fraction of the From 2c2b65553560c799e1c9f050c27d26960fa69d5c Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 15 Jul 2022 10:58:35 -0400 Subject: [PATCH 11/25] Check if the number of impactors in the subpixel diffusion step is bigger than the I8B maximum size. If so, skip, though it may be beneficial to let the user know they may need to truncate their production population. --- src/crater/crater_subpixel_diffusion.f90 | 1 + 1 file changed, 1 insertion(+) diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90 index f3c3e27e..fe8b2b30 100644 --- a/src/crater/crater_subpixel_diffusion.f90 +++ b/src/crater/crater_subpixel_diffusion.f90 @@ -94,6 +94,7 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin) end if lambda = dN * domain%parea + if (lambda > 1.0_DP * huge(N)) cycle ! Too many impactors. ! Don't parallelize the random do j = 1,user%gridsize From 79d94be8d0c8e24ac73b72207f5cc0a1396fdd64 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Fri, 15 Jul 2022 16:09:08 -0400 Subject: [PATCH 12/25] Cleaned up file io for rego tracking. File names now conform to the standard CTEM style, removed redundant set of top files, and added rego files (age, melt, comp, and stack) to the redirected save files when doregotrack is turned on. --- python/ctem/ctem/driver.py | 7 ++ src/crater/crater_populate.f90 | 2 +- src/globals/module_globals.f90 | 3 +- src/io/io_read_regotrack.f90 | 51 ++++++++------- src/io/io_write_regotrack.f90 | 115 +++++++-------------------------- 5 files changed, 58 insertions(+), 120 deletions(-) diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py index 52174cc8..afdb2aa6 100644 --- a/python/ctem/ctem/driver.py +++ b/python/ctem/ctem/driver.py @@ -59,6 +59,11 @@ def __init__(self, param_file="ctem.in", isnew=True): 'ejc': 'surface_ejc.dat', 'pos': 'surface_pos.dat', 'time': 'surface_time.dat', + 'stack': 'surface_stacknum.dat', + 'rego': 'surface_rego.dat', + 'melt': 'surface_melt.dat', + 'comp': 'surface_comp.dat', + 'age' : 'surface_age.dat', 'ocum': 'ocumulative.dat', 'odist': 'odistribution.dat', 'pdist': 'pdistribution.dat', @@ -309,6 +314,8 @@ def process_output(self): if (self.user['savetruelist'].upper() == 'T'): self.redirect_outputs(['tcum'], 'dist') self.redirect_outputs(['impmass'], 'misc') + if (self.user['saverego'].upper() == 'T') : + self.redirect_outputs(['stack','rego','age','melt','comp'], 'rego') diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index c3c779c3..50ac5760 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -89,7 +89,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! doregotrack & age simulation test real(DP) :: melt, clock, age, thick real(SP),dimension(user%gridsize, user%gridsize) :: agetop - real(SP),dimension(60) :: agetot + real(SP),dimension(MAXAGEBINS) :: agetot type(regolisttype),pointer :: current => null() real(DP) :: age_resolution diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index f29b7600..6da3fb25 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -63,6 +63,7 @@ module module_globals real(DP),parameter :: PF = 5.0e9 ! The shock pressure exceeding the hungoit elastic limit of common geological materials in Pa. (5 GPa) real(DP),parameter :: RAD_GP = 20.0_DP ! The maximum radial position of producing impact glass spherules within a transient crater (unit of crater radii, crater%rad) + type regodatatype real(SP),dimension(MAXAGEBINS) :: age real(DP) :: thickness @@ -275,7 +276,7 @@ module module_globals character(*),parameter :: TIMEFILE = 'surface_time.dat' character(*),parameter :: EJCOVFILE = 'surface_ejc.dat' character(*),parameter :: DEMFILE = 'surface_dem.dat' -character(*),parameter :: REGOFILE = 'surface_regotop.dat' +character(*),parameter :: REGOFILE = 'surface_rego.dat' character(*),parameter :: MELTFILE = 'surface_melt.dat' character(*),parameter :: COMPFILE = 'surface_comp.dat' character(*),parameter :: STACKNUMFILE = 'surface_stacknum.dat' diff --git a/src/io/io_read_regotrack.f90 b/src/io/io_read_regotrack.f90 index ed93300c..3acdbb14 100644 --- a/src/io/io_read_regotrack.f90 +++ b/src/io/io_read_regotrack.f90 @@ -28,14 +28,14 @@ subroutine io_read_regotrack(user,surf) ! Internals integer(I4B), parameter :: LUN=7 - integer(I4B), parameter :: LUM=8 - integer(I4B), parameter :: LUP=9 - integer(I4B), parameter :: LUC=10 - integer(I4B), parameter :: LUA=11 + integer(I4B), parameter :: FMELT = 10 + integer(I4B), parameter :: FREGO = 11 + integer(I4B), parameter :: FCOMP = 12 + integer(I4B), parameter :: FAGE = 13 real(DP),dimension(user%gridsize,user%gridsize) :: regotop,melt,comp - real(SP),dimension(user%gridsize,user%gridsize,60) :: age + real(SP),dimension(user%gridsize,user%gridsize,MAXAGEBINS) :: age integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num - real(DP), dimension(:),allocatable :: regotopi,melti,compi,agei + real(DP), dimension(:), allocatable :: regotopi,melti,compi,agei type(regodatatype) :: newsurfi integer(I4B) :: ioerr,i,j,k,q,itmp integer(kind=8) :: recsize @@ -46,34 +46,34 @@ subroutine io_read_regotrack(user,surf) ! Open a file for obtaining the number of stacks that is stored in each linked list ioerr = 0 recsize = sizeof(itmp) * user%gridsize * user%gridsize - open(LUP,file=STACKNUMFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr) + open(LUN,file=STACKNUMFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(STACKNUMFILE)) stop end if - read(LUP,rec=1) stacks_num - close(LUP) + read(LUN,rec=1) stacks_num + close(LUN) ! Open files for regolith thickness/melt fraction - open(LUN,file=REGOFILE,status='old',form='unformatted',iostat=ioerr) + open(FREGO,file=REGOFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(REGOFILE)) stop end if - open(LUC,file=COMPFILE,status='old',form='unformatted',iostat=ioerr) + open(FCOMP,file=COMPFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(COMPFILE)) stop end if - open(LUM,file=MELTFILE,status='old',form='unformatted',iostat=ioerr) + open(FMELT,file=MELTFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(MELTFILE)) stop end if - open(LUA,file=AGEFILE,status='old',form='unformatted',iostat=ioerr) + open(FAGE,file=AGEFILE,status='old',form='unformatted',iostat=ioerr) if (ioerr/=0) then write(*,*) 'Error! Cannot read file ',trim(adjustl(MELTFILE)) stop @@ -89,18 +89,18 @@ subroutine io_read_regotrack(user,surf) allocate(regotopi(stacks_num(i,j))) allocate(compi(stacks_num(i,j))) allocate(melti(stacks_num(i,j))) - allocate(agei(60 * stacks_num(i,j))) + allocate(agei(MAXAGEBINS * stacks_num(i,j))) do k=1,stacks_num(i,j) - read(LUN) regotop(i,j) + read(FREGO) regotop(i,j) regotopi(k) = regotop(i,j) - read(LUC) comp(i,j) + read(FCOMP) comp(i,j) compi(k) = comp(i,j) - read(LUM) melt(i,j) + read(FMELT) melt(i,j) melti(k) = melt(i,j) - read(LUA) age(i,j,:) - do q=1,60 - agei(60*k - (60-q)) = age(i,j,q) + read(FAGE) age(i,j,:) + do q=1,MAXAGEBINS + agei(MAXAGEBINS*k - (MAXAGEBINS-q)) = age(i,j,q) end do end do @@ -108,8 +108,8 @@ subroutine io_read_regotrack(user,surf) newsurfi%thickness = regotopi(k) newsurfi%comp = compi(k) newsurfi%meltfrac = melti(k) - do q=1,60 - newsurfi%age(q) = agei(60*k-(60-q)) + do q=1,MAXAGEBINS + newsurfi%age(q) = agei(MAXAGEBINS*k-(MAXAGEBINS-q)) end do call util_push(surf(i,j)%regolayer,newsurfi) end do @@ -118,8 +118,9 @@ subroutine io_read_regotrack(user,surf) end do end do - close(LUN) - close(LUC) - close(LUM) + close(FMELT) + close(FREGO) + close(FCOMP) + close(FAGE) return end subroutine io_read_regotrack diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index a6f3c2e1..5f522321 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -27,14 +27,12 @@ subroutine io_write_regotrack(user,surf) ! Regotrack Internals integer(I4B) :: i,j,k - integer(I4B), parameter :: LUN=7 - integer(I4B), parameter :: LUM=8 - integer(I4B), parameter :: LUC=9 - integer(I4B), parameter :: LUA=10 + integer(I4B), parameter :: LUN = 7 + integer(I4B), parameter :: FMELT = 10 + integer(I4B), parameter :: FREGO = 11 + integer(I4B), parameter :: FCOMP = 12 + integer(I4B), parameter :: FAGE = 13 type(regolisttype),pointer :: current => null() - real(DP),dimension(user%gridsize,user%gridsize) :: regotop,comp,melt - real(SP),dimension(user%gridsize,user%gridsize) :: agetop - real(SP),dimension(user%gridsize,user%gridsize,60) :: age integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num integer(kind=8) :: recsize real(DP) :: dtmp @@ -42,104 +40,35 @@ subroutine io_write_regotrack(user,surf) integer(I4B) :: itmp real(DP),dimension(user%gridsize,user%gridsize) :: comptop, rego real(DP),dimension(:),allocatable :: marehisto - real(DP) :: mare, z - ! Mixing - !real(DP),parameter :: zmix = 0.0_DP - !real(DP) :: z, zmare - - ! Output multiple "comphisto" files - character(len=255) :: fname - character(len=255), parameter :: clockfile = 'tic-toc.dat' - integer(I4B) :: tictoc - logical :: exist - ! Executable code - ! Output mulitple "comphisto" files - inquire(file=clockfile, exist=exist) - if (exist) then - open(LUN,file=clockfile,status='old') - read(LUN,*) tictoc - else - write(*,*) clockfile,' is missing!' - end if - tictoc = tictoc + 1 - close(LUN) - open(LUN,file=clockfile,status='replace') - write(LUN,*) tictoc - close(LUN) - - write(fname,'(a,i4.4)') 'surface_melt',tictoc - open(LUN,file=fname,status='replace',form='unformatted') - write(fname,'(a,i4.4)') 'surface_rego',tictoc - open(LUM,file=fname,status='replace',form='unformatted') - write(fname,'(a,i4.4)') 'surface_comp',tictoc - open(LUC,file=fname,status='replace',form='unformatted') - write(fname,'(a,i4.4)') 'surface_age',tictoc - open(LUA,file=fname,status='replace',form='unformatted') + open(FMELT,file=MELTFILE,status='replace',form='unformatted') + open(FREGO,file=REGOFILE,status='replace',form='unformatted') + open(FCOMP,file=COMPFILE,status='replace',form='unformatted') + open(FAGE,file=AGEFILE,status='replace',form='unformatted') + stacks_num(:,:) = 0 do j=1,user%gridsize do i=1,user%gridsize - stacks_num(i,j) = 0 current => surf(i,j)%regolayer - comptop(i,j) = current%regodata%comp - rego(i,j) = current%regodata%thickness - agetop(i,j) = current%regodata%age(1) do - if (.not. associated(current)) exit - stacks_num(i,j) = stacks_num(i,j) + 1 - regotop(i,j) = current%regodata%thickness - comp(i,j) = current%regodata%comp - melt(i,j) = current%regodata%meltfrac - age(i,j,:)= current%regodata%age(:) - write(LUM) regotop(i,j) - write(LUC) comp(i,j) - write(LUN) melt(i,j) - write(LUA) age(i,j,:) - current => current%next + if (.not. associated(current)) exit ! We've reached the bottom of the linked list + stacks_num(i,j) = stacks_num(i,j) + 1 + write(FMELT) current%regodata%meltfrac + write(FREGO) current%regodata%comp + write(FCOMP) current%regodata%meltfrac + write(FAGE) current%regodata%age(:) + current => current%next end do end do end do - close(LUN) - close(LUM) - close(LUC) - close(LUA) - - recsize = sizeof(dtmp) * user%gridsize * user%gridsize - open(LUN,file='agetop.dat',status='replace',form='unformatted',recl=recsize,access='direct') - write(LUN,rec=1) agetop - close(LUN) - - recsize = sizeof(dtmp) * user%gridsize * user%gridsize - open(LUN,file='comptop.dat',status='replace',form='unformatted',recl=recsize,access='direct') - write(LUN,rec=1) comptop - close(LUN) - - allocate(marehisto(user%gridsize)) - !open(LUN,file='comphisto',status='replace') ! Output comphisto one time - write(fname,'(a,i4.4)') 'comphisto',tictoc - open(LUN,file=fname,status='replace') - - do i=1,user%gridsize - marehisto(i) = 0.0_DP - mare = 0.0_DP - do j=1,user%gridsize - mare = mare + comptop(i,j) - end do - marehisto(i) = mare/real(user%gridsize) - write(LUN,*) real(i-user%gridsize/2)*user%pix/1000.0,marehisto(i)*100.0 - end do - close(LUN) - deallocate(marehisto) - - recsize = sizeof(dtmp) * user%gridsize * user%gridsize - open(LUN,file='regotop.dat',status='replace',form='unformatted',recl=recsize,access='direct') - write(LUN,rec=1) rego - close(LUN) + close(FMELT) + close(FREGO) + close(FCOMP) + close(FAGE) recsize = sizeof(itmp) * user%gridsize * user%gridsize - write(fname,'(a,i4.4)') 'surface_stacknum',tictoc - open(LUN,file=fname,status='replace',form='unformatted',recl=recsize,access='direct') + open(LUN,file=STACKNUMFILE,status='replace',form='unformatted',recl=recsize,access='direct') write(LUN,rec=1) stacks_num close(LUN) From 758791e3cc10df377091e5384d8ce04a4000d0eb Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 18 Jul 2022 13:09:00 -0400 Subject: [PATCH 13/25] Fixed bug in which the regolith composition data was being written to the regolith thickness file --- src/io/io_write_regotrack.f90 | 4 ++-- 1 file changed, 2 insertions(+), 2 deletions(-) diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index 5f522321..d76950e5 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -55,8 +55,8 @@ subroutine io_write_regotrack(user,surf) if (.not. associated(current)) exit ! We've reached the bottom of the linked list stacks_num(i,j) = stacks_num(i,j) + 1 write(FMELT) current%regodata%meltfrac - write(FREGO) current%regodata%comp - write(FCOMP) current%regodata%meltfrac + write(FREGO) current%regodata%thickness + write(FCOMP) current%regodata%comp write(FAGE) current%regodata%age(:) current => current%next end do From 5be9befa6682650f58055366b649beb7eb9381b9 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 18 Jul 2022 13:09:32 -0400 Subject: [PATCH 14/25] Reduced the size of the regolith buffer to prevent underflow. --- src/util/util_init_list.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/util_init_list.f90 b/src/util/util_init_list.f90 index 672f9c11..46249dc0 100644 --- a/src/util/util_init_list.f90 +++ b/src/util/util_init_list.f90 @@ -73,7 +73,7 @@ subroutine util_init_list(regolayer,initstat) if (allocstat == 0) then initstat = .true. nullify(regolayer%next) - regolayer%regodata%thickness = VBIG ! This generates a buffer layer that the model should never reach if the run is structured properly + regolayer%regodata%thickness = VBIG / 1e6_DP ! This generates a buffer layer that the model should never reach if the run is structured properly regolayer%regodata%comp = 0.0_DP regolayer%regodata%meltfrac = 0.0_DP regolayer%regodata%porosity = 0.0_DP From 53d23834d9663fa5c07e0f33daeaa8630380c273 Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Mon, 18 Jul 2022 13:12:16 -0400 Subject: [PATCH 15/25] Changed RAD_GP to 1 per Ya-Huei's request --- src/globals/module_globals.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/globals/module_globals.f90 b/src/globals/module_globals.f90 index 6da3fb25..8013b051 100644 --- a/src/globals/module_globals.f90 +++ b/src/globals/module_globals.f90 @@ -61,7 +61,7 @@ module module_globals real(DP),parameter :: RCONT = 2.25267_DP ! Coefficient of continuous ejecta size power law from Moore et al. (1974) - scaled from km to m real(DP),parameter :: EXPCONT = 1.006_DP ! Exponent of continuous ejecta size power law from Moore et al. (1974) real(DP),parameter :: PF = 5.0e9 ! The shock pressure exceeding the hungoit elastic limit of common geological materials in Pa. (5 GPa) -real(DP),parameter :: RAD_GP = 20.0_DP ! The maximum radial position of producing impact glass spherules within a transient crater (unit of crater radii, crater%rad) +real(DP),parameter :: RAD_GP = 1.0_DP ! The maximum radial position of producing impact glass spherules within a transient crater (unit of crater radii, crater%rad) type regodatatype From 8dd3d25b81867a0a57d57d5605dbdce356f399b1 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 18 Jul 2022 13:14:32 -0400 Subject: [PATCH 16/25] Reduced the buffer size even more for underflow exception --- src/util/util_init_list.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/util/util_init_list.f90 b/src/util/util_init_list.f90 index 46249dc0..f1e50754 100644 --- a/src/util/util_init_list.f90 +++ b/src/util/util_init_list.f90 @@ -73,7 +73,7 @@ subroutine util_init_list(regolayer,initstat) if (allocstat == 0) then initstat = .true. nullify(regolayer%next) - regolayer%regodata%thickness = VBIG / 1e6_DP ! This generates a buffer layer that the model should never reach if the run is structured properly + regolayer%regodata%thickness = sqrt(VBIG) ! This generates a buffer layer that the model should never reach if the run is structured properly regolayer%regodata%comp = 0.0_DP regolayer%regodata%meltfrac = 0.0_DP regolayer%regodata%porosity = 0.0_DP From 62318b1de9a994c8e7ec1df36c057108822bef47 Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 18 Jul 2022 13:24:01 -0400 Subject: [PATCH 17/25] Added capability of reading in linked list binary data (not the age bin data yet, though) --- python/ctem/ctem/util.py | 26 ++++++++++++++++++++++++-- 1 file changed, 24 insertions(+), 2 deletions(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index d77b84f9..c9ab281f 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -6,6 +6,7 @@ import matplotlib.pyplot as plt import re from tempfile import mkstemp +from scipy.io import FortranFile # Set pixel scaling common for image writing, at 1 pixel/ array element dpi = 72.0 @@ -263,16 +264,37 @@ def read_user_input(user): return user -def read_unformatted_binary(filename, gridsize): +def read_unformatted_binary(filename, gridsize, kind='DP'): # Read unformatted binary files created by Fortran # For use with surface ejecta and surface dem data files - dt = np.float + if kind == 'DP': + dt = np.dtype('f8') + elif kind == 'SP': + dt = np.dtype('f4') + elif kind == 'I4B': + dt = np.dtype(' Date: Mon, 18 Jul 2022 13:35:59 -0400 Subject: [PATCH 18/25] Fixed array reading issue with the linked list reader --- python/ctem/ctem/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index c9ab281f..3347de17 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -291,7 +291,7 @@ def read_linked_list_binary(filename, stackname, gridsize): for s in np.arange(stack[j, i]): d = f.read_reals(np.float64) datastack.append(d) - data[j, i] = d + data[j, i] = np.asarray(datastack).T return data From c33b2e9e245ad471dcb76b5c06aae6533b3bfb74 Mon Sep 17 00:00:00 2001 From: David A Minton Date: Mon, 18 Jul 2022 14:28:25 -0400 Subject: [PATCH 19/25] Changed the regotracking write tool to write out whole arrays of stacks rather than going layer-by-layer --- src/io/io_write_regotrack.f90 | 30 +++++++++++++++++++++++++----- 1 file changed, 25 insertions(+), 5 deletions(-) diff --git a/src/io/io_write_regotrack.f90 b/src/io/io_write_regotrack.f90 index d76950e5..69b265ce 100644 --- a/src/io/io_write_regotrack.f90 +++ b/src/io/io_write_regotrack.f90 @@ -34,10 +34,12 @@ subroutine io_write_regotrack(user,surf) integer(I4B), parameter :: FAGE = 13 type(regolisttype),pointer :: current => null() integer(I4B),dimension(user%gridsize,user%gridsize) :: stacks_num + real(DP),dimension(:),allocatable :: meltfrac, thickness, comp + real(SP),dimension(:,:),allocatable :: age integer(kind=8) :: recsize real(DP) :: dtmp real(SP) :: stmp - integer(I4B) :: itmp + integer(I4B) :: itmp, N real(DP),dimension(user%gridsize,user%gridsize) :: comptop, rego real(DP),dimension(:),allocatable :: marehisto @@ -47,6 +49,7 @@ subroutine io_write_regotrack(user,surf) open(FCOMP,file=COMPFILE,status='replace',form='unformatted') open(FAGE,file=AGEFILE,status='replace',form='unformatted') + ! First pass to get stack numbers stacks_num(:,:) = 0 do j=1,user%gridsize do i=1,user%gridsize @@ -54,14 +57,31 @@ subroutine io_write_regotrack(user,surf) do if (.not. associated(current)) exit ! We've reached the bottom of the linked list stacks_num(i,j) = stacks_num(i,j) + 1 - write(FMELT) current%regodata%meltfrac - write(FREGO) current%regodata%thickness - write(FCOMP) current%regodata%comp - write(FAGE) current%regodata%age(:) current => current%next end do end do end do + + ! Second pass to get data and save it + do j=1,user%gridsize + do i=1,user%gridsize + current => surf(i,j)%regolayer + N = stacks_num(i,j) + allocate(meltfrac(N),thickness(N),comp(N),age(MAXAGEBINS,N)) + do k=1,N + meltfrac(k) = current%regodata%meltfrac + thickness(k) = current%regodata%thickness + comp(k) = current%regodata%comp + age(:,k) = current%regodata%age(:) + current => current%next + end do + write(FMELT) meltfrac(:) + write(FREGO) thickness(:) + write(FCOMP) comp(:) + write(FAGE) age(:,:) + deallocate(meltfrac,thickness,comp,age) + end do + end do close(FMELT) close(FREGO) close(FCOMP) From 6cde6346b1f4f3e776798b8efa369400d199e231 Mon Sep 17 00:00:00 2001 From: David Minton Date: Mon, 18 Jul 2022 14:34:11 -0400 Subject: [PATCH 20/25] Updated utility to read in linked list data directly without needing the stack file --- python/ctem/ctem/util.py | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index 3347de17..167d0fc5 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -281,17 +281,12 @@ def read_unformatted_binary(filename, gridsize, kind='DP'): return data -def read_linked_list_binary(filename, stackname, gridsize): - stack = read_unformatted_binary(stackname,gridsize,kind='I4B') +def read_linked_list_binary(filename, gridsize): data = np.empty((gridsize,gridsize),dtype="object") with FortranFile(filename, 'r') as f: for i in np.arange(gridsize): for j in np.arange(gridsize): - datastack = [] - for s in np.arange(stack[j, i]): - d = f.read_reals(np.float64) - datastack.append(d) - data[j, i] = np.asarray(datastack).T + data[j, i] = f.read_reals(np.float64) return data From 03ecbb083c6d6822875f3c7c3e756cff76ae746c Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 19 Jul 2022 11:32:17 -0400 Subject: [PATCH 21/25] Changed incorrect array index order when reading linked list binary data --- python/ctem/ctem/util.py | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index 167d0fc5..070ed3f3 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -286,7 +286,7 @@ def read_linked_list_binary(filename, gridsize): with FortranFile(filename, 'r') as f: for i in np.arange(gridsize): for j in np.arange(gridsize): - data[j, i] = f.read_reals(np.float64) + data[i, j] = f.read_reals(np.float64) return data From 654054d62d11823f64fa905283cb076e0a8f8a4b Mon Sep 17 00:00:00 2001 From: Austin Blevins Date: Tue, 19 Jul 2022 11:35:38 -0400 Subject: [PATCH 22/25] Pushed Ya-Huei's change to crater_populate that ensured ages weren't in a single bin --- src/crater/crater_populate.f90 | 9 ++------- 1 file changed, 2 insertions(+), 7 deletions(-) diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90 index 50ac5760..7ed11620 100644 --- a/src/crater/crater_populate.f90 +++ b/src/crater/crater_populate.f90 @@ -89,7 +89,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! doregotrack & age simulation test real(DP) :: melt, clock, age, thick real(SP),dimension(user%gridsize, user%gridsize) :: agetop - real(SP),dimension(MAXAGEBINS) :: agetot + real(SP),dimension(60) :: agetot type(regolisttype),pointer :: current => null() real(DP) :: age_resolution @@ -240,12 +240,6 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt ! Place crater onto the surface call crater_emplace(user,surf,crater,domain,ejbmass) - if (abs(ejbmass) < tiny(ejbmass)) then ! Crater has no topography. Discard and move on. - ntrue = ntrue -1 - mass = mass - crater%impmass - cycle - end if - call ejecta_distance_estimate(user,crater,domain,crater%ejdis) ! Fast but imprecise estimate of the total ejecta distance ! For very steep size distributions, only a fraction of the @@ -319,6 +313,7 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt call crater_superdomain(user,surf,age,age_resolution,prod,nflux,domain,finterval) call regolith_depth_model(user,domain,finterval,nflux,p) call regolith_subcrater_mix(user,surf,domain,nflux,finterval,p) + age = age - finterval * user%interval end if ! Do periodic subpixel processes on the whole grid From e65999901ec20117cb067a8c01480f38d8b2b57e Mon Sep 17 00:00:00 2001 From: David Minton Date: Tue, 19 Jul 2022 11:39:58 -0400 Subject: [PATCH 23/25] Added optional kind parameter to linked list binary reader in order to select the appropriate CTEM data type --- python/ctem/ctem/util.py | 12 ++++++++++-- 1 file changed, 10 insertions(+), 2 deletions(-) diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py index 070ed3f3..da9d0b8b 100644 --- a/python/ctem/ctem/util.py +++ b/python/ctem/ctem/util.py @@ -281,12 +281,20 @@ def read_unformatted_binary(filename, gridsize, kind='DP'): return data -def read_linked_list_binary(filename, gridsize): +def read_linked_list_binary(filename, gridsize, kind='DP'): + if kind == 'DP': + dt = np.dtype('f8') + elif kind == 'SP': + dt = np.dtype('f4') + elif kind == 'I4B': + dt = np.dtype(' Date: Sun, 24 Jul 2022 22:45:43 -0400 Subject: [PATCH 24/25] fixed a bug in error handling --- src/io/io_read_surf.f90 | 2 +- 1 file changed, 1 insertion(+), 1 deletion(-) diff --git a/src/io/io_read_surf.f90 b/src/io/io_read_surf.f90 index af611eea..4b52660d 100644 --- a/src/io/io_read_surf.f90 +++ b/src/io/io_read_surf.f90 @@ -69,7 +69,7 @@ subroutine io_read_surf(user,surf) open(LUN,file=TIMEFILE,status='old',form='unformatted',recl=recsize,access='direct',iostat=ioerr) if (ioerr/=0) then - write(*,*) 'Error! Cannot read file ',trim(adjustl(DIAMFILE)) + write(*,*) 'Error! Cannot read file ',trim(adjustl(TIMEFILE)) stop end if do i=1,user%numlayers From 75e57322db1f9b4480558337e0f37810b5912df5 Mon Sep 17 00:00:00 2001 From: David Minton Date: Fri, 23 Sep 2022 16:08:19 -0400 Subject: [PATCH 25/25] Disabled everything except the rim code in the realistic morphology function --- src/crater/crater_realistic_topography.f90 | 124 ++++++++++----------- 1 file changed, 62 insertions(+), 62 deletions(-) diff --git a/src/crater/crater_realistic_topography.f90 b/src/crater/crater_realistic_topography.f90 index 115601ff..07a1cf5e 100644 --- a/src/crater/crater_realistic_topography.f90 +++ b/src/crater/crater_realistic_topography.f90 @@ -112,14 +112,14 @@ end subroutine crater_realistic_slope_texture select case(crater%morphtype) case("COMPLEX","PEAKRING","MULTIRING") call complex_terrace(user,surf,crater,deltaMtot) - call complex_wall_texture(user,surf,crater,domain,deltaMtot) - call complex_floor(user,surf,crater,deltaMtot) - call complex_peak(user,surf,crater,deltaMtot) + !call complex_wall_texture(user,surf,crater,domain,deltaMtot) + !call complex_floor(user,surf,crater,deltaMtot) + !call complex_peak(user,surf,crater,deltaMtot) end select ! Retrieve the size of the ejecta dem and correct for indexing - inc = (size(ejecta_dem,1) - 1) / 2 - call ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) + !inc = (size(ejecta_dem,1) - 1) / 2 + !call ejecta_texture(user,surf,crater,deltaMtot,inc,ejecta_dem) ! Do a final pass of the slope collapse with a shallower slope than normal to smooth out all of the sharp edges @@ -448,63 +448,63 @@ subroutine complex_terrace(user,surf,crater,deltaMtot) crater%maxinc = max(crater%maxinc,inc) flr = crater%floordiam / crater%fcrat - do terrace = 1, nterraces - router = flr + (1._DP - flr) * (terrace / real(nterraces,kind=DP))**(terracefac) ! The radius of the outer edge of the terrace - rinner = flr + (1._DP - flr) * ((terrace - 1) / real(nterraces,kind=DP))**(terracefac) ! The radius of the inner edge of the terrace - - do j = -inc,inc - do i = -inc,inc - xpi = crater%xlpx + i - ypi = crater%ylpx + j - - xbar = xpi * user%pix - crater%xl - ybar = ypi * user%pix - crater%yl - - ! periodic boundary conditions - call util_periodic(xpi,ypi,user%gridsize) - newdem = surf(xpi,ypi)%dem - - r = sqrt(xbar**2 + ybar**2) / crater%frad - - ! Make scalloped terraces - znoise = scallop_width**(1._DP / (2 * scallop_p)) - xynoise = (nscallops / PI) / crater%fcrat - dnoise = util_perlin_noise(xynoise * xbar + terrace * offset * rn(1), & - xynoise * ybar + terrace * offset * rn(2)) * znoise - noise = (dnoise**2)**scallop_p - hprof = (r / router)**(-1) - - ! Make textured floor of terrace - tnoise = 0.0_DP - do octave = 1, num_oct_tfloor - xynoise = xy_size_tfloor * freq_tfloor ** (octave - 1) - znoise = noise_height_tfloor * (pers_tfloor ) ** (octave - 1) - tnoise = tnoise + util_perlin_noise(xynoise * xbar + offset * rn(1), & - xynoise * ybar + offset * rn(2))* znoise - end do - - isterrace = 1.0_DP - noise < hprof - - if (r < 1.0_DP) then - tprof = crater_profile(user,crater,rinner) * (r/rinner)**h_scallop_profile + tnoise ! This is the floor profile that replaces the old one at each terrace - else - tprof = (crater_profile(user,crater,rinner) + crater%ejrim * (rinner)**(-EJPROFILE)) + tnoise ! This is the floor profile that replaces the old one at each terrace - end if - - if (isterrace) then - newdem = min(tprof,newdem) - elchange = newdem - surf(xpi,ypi)%dem - deltaMtot = deltaMtot + elchange - surf(xpi,ypi)%dem = newdem - ! Save the minimum elevation below the floor - dfloor = crater%melev - crater%floordepth - newdem - if (dfloor > 0.0_DP) upshift = max(upshift,dfloor) - end if - - end do - end do - - end do + ! do terrace = 1, nterraces + ! router = flr + (1._DP - flr) * (terrace / real(nterraces,kind=DP))**(terracefac) ! The radius of the outer edge of the terrace + ! rinner = flr + (1._DP - flr) * ((terrace - 1) / real(nterraces,kind=DP))**(terracefac) ! The radius of the inner edge of the terrace + + ! do j = -inc,inc + ! do i = -inc,inc + ! xpi = crater%xlpx + i + ! ypi = crater%ylpx + j + + ! xbar = xpi * user%pix - crater%xl + ! ybar = ypi * user%pix - crater%yl + + ! ! periodic boundary conditions + ! call util_periodic(xpi,ypi,user%gridsize) + ! newdem = surf(xpi,ypi)%dem + + ! r = sqrt(xbar**2 + ybar**2) / crater%frad + + ! ! Make scalloped terraces + ! znoise = scallop_width**(1._DP / (2 * scallop_p)) + ! xynoise = (nscallops / PI) / crater%fcrat + ! dnoise = util_perlin_noise(xynoise * xbar + terrace * offset * rn(1), & + ! xynoise * ybar + terrace * offset * rn(2)) * znoise + ! noise = (dnoise**2)**scallop_p + ! hprof = (r / router)**(-1) + + ! ! Make textured floor of terrace + ! tnoise = 0.0_DP + ! do octave = 1, num_oct_tfloor + ! xynoise = xy_size_tfloor * freq_tfloor ** (octave - 1) + ! znoise = noise_height_tfloor * (pers_tfloor ) ** (octave - 1) + ! tnoise = tnoise + util_perlin_noise(xynoise * xbar + offset * rn(1), & + ! xynoise * ybar + offset * rn(2))* znoise + ! end do + + ! isterrace = 1.0_DP - noise < hprof + + ! if (r < 1.0_DP) then + ! tprof = crater_profile(user,crater,rinner) * (r/rinner)**h_scallop_profile + tnoise ! This is the floor profile that replaces the old one at each terrace + ! else + ! tprof = (crater_profile(user,crater,rinner) + crater%ejrim * (rinner)**(-EJPROFILE)) + tnoise ! This is the floor profile that replaces the old one at each terrace + ! end if + + ! if (isterrace) then + ! newdem = min(tprof,newdem) + ! elchange = newdem - surf(xpi,ypi)%dem + ! deltaMtot = deltaMtot + elchange + ! surf(xpi,ypi)%dem = newdem + ! ! Save the minimum elevation below the floor + ! dfloor = crater%melev - crater%floordepth - newdem + ! if (dfloor > 0.0_DP) upshift = max(upshift,dfloor) + ! end if + + ! end do + ! end do + + ! end do ! Now make the scalloped rim do j = -inc,inc