From d7630ccd3225e46f2863e755216c51aeb1324461 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Wed, 23 Feb 2022 15:19:12 -0500
Subject: [PATCH 01/19] Replaced Python code with revised module

---
 python/ctem/ctem/__init__.py           |   1 +
 python/ctem/ctem/ctem_draw_surf.py     |  69 -----
 python/ctem/ctem/ctem_driver.py        | 229 ---------------
 python/ctem/ctem/ctem_io_readers.py    | 144 ----------
 python/ctem/ctem/ctem_io_writers.py    | 264 -----------------
 python/ctem/ctem/ctem_viewer_3d.py     |  92 ------
 python/ctem/ctem/driver.py             | 262 +++++++++++++++++
 python/ctem/ctem/io.py                 | 382 +++++++++++++++++++++++++
 python/ctem/ctem/viewer3d.py           | 122 ++++++++
 python/ctem/tests/testio/ctem.in       |  70 +++++
 python/ctem/tests/testio/testio.py     |   7 +
 python/ctem/tests/viewer3d/polytest.py |  13 +-
 12 files changed, 851 insertions(+), 804 deletions(-)
 delete mode 100644 python/ctem/ctem/ctem_draw_surf.py
 delete mode 100644 python/ctem/ctem/ctem_driver.py
 delete mode 100644 python/ctem/ctem/ctem_io_readers.py
 delete mode 100644 python/ctem/ctem/ctem_io_writers.py
 delete mode 100644 python/ctem/ctem/ctem_viewer_3d.py
 create mode 100644 python/ctem/ctem/driver.py
 create mode 100644 python/ctem/ctem/io.py
 create mode 100644 python/ctem/ctem/viewer3d.py
 create mode 100755 python/ctem/tests/testio/ctem.in
 create mode 100644 python/ctem/tests/testio/testio.py

diff --git a/python/ctem/ctem/__init__.py b/python/ctem/ctem/__init__.py
index 4774c6af..165a1bb0 100644
--- a/python/ctem/ctem/__init__.py
+++ b/python/ctem/ctem/__init__.py
@@ -1,3 +1,4 @@
 # from ctem.ctem_io_readers import *
 # from ctem.ctem_io_writers import *
 # from ctem.ctem_driver import *
+from ctem.driver import *
\ No newline at end of file
diff --git a/python/ctem/ctem/ctem_draw_surf.py b/python/ctem/ctem/ctem_draw_surf.py
deleted file mode 100644
index 74bf81e1..00000000
--- a/python/ctem/ctem/ctem_draw_surf.py
+++ /dev/null
@@ -1,69 +0,0 @@
-#!/usr/local/bin/python
-#
-#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 CTEM modules
-import ctem_io_readers
-import ctem_io_writers
-
-#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}
-
-#Read ctem.in to initialize parameter values based on user input
-ctem_io_readers.read_ctemin(parameters,notset)
-ctem_io_writers.create_dir_structure(parameters)
-
-dem_file = parameters['workingdir'] + 'surface_dem.dat'
-surface_dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize'])
-
-#Set up data arrays
-seedarr = numpy.zeros(100, dtype = numpy.int)
-seedarr[0] = parameters['seed']
-
-#Read ctem.dat file
-ctem_io_readers.read_ctemdat(parameters, seedarr)
-
-#Write surface dem
-ctem_io_writers.image_dem(parameters, surface_dem)
\ No newline at end of file
diff --git a/python/ctem/ctem/ctem_driver.py b/python/ctem/ctem/ctem_driver.py
deleted file mode 100644
index 29057053..00000000
--- a/python/ctem/ctem/ctem_driver.py
+++ /dev/null
@@ -1,229 +0,0 @@
-#!/usr/local/bin/python
-#
-#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 CTEM modules
-import ctem_io_readers
-import ctem_io_writers
-
-#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}
-
-#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)
-
-#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, 
-                              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/python/ctem/ctem/ctem_io_readers.py b/python/ctem/ctem/ctem_io_readers.py
deleted file mode 100644
index 84cf7e06..00000000
--- a/python/ctem/ctem/ctem_io_readers.py
+++ /dev/null
@@ -1,144 +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 ('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/python/ctem/ctem/ctem_io_writers.py b/python/ctem/ctem/ctem_io_writers.py
deleted file mode 100644
index 7bf21ce4..00000000
--- a/python/ctem/ctem/ctem_io_writers.py
+++ /dev/null
@@ -1,264 +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
-from matplotlib.colors import LightSource
-import matplotlib.cm as cm
-
-#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, DEM):
-    dpi = 300.0  # 72.0
-    pix = parameters['pix']
-    gridsize = parameters['gridsize']
-    ve = 1.0
-    azimuth = 300.0 #parameters['azimuth']
-    solar_angle = 20.0 #parameters['solar_angle']
-
-    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
-    dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix)
-
-    # Generate image to put into an array
-    height = gridsize / dpi
-    width = gridsize / dpi
-    fig = matplotlib.pyplot.figure(figsize=(width, height), dpi=dpi)
-    ax = matplotlib.pyplot.axes([0, 0, 1, 1])
-    ax.imshow(dem_img, interpolation="nearest", cmap='gray',vmin=0.0,vmax=1.0)
-    matplotlib.pyplot.axis('off')
-    # Save image to file
-    filename = parameters['workingdir'] + 'surf' + os.sep + "surf%06d.png" % parameters['ncount']
-    matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0)
-
-    return
-
-def image_shaded_relief(parameters, DEM):
-    dpi = 300.0 #72.0
-    pix = parameters['pix']
-    gridsize = parameters['gridsize']
-    ve = 1.0
-    mode = 'overlay'
-    azimuth = 300.0 #parameters['azimuth']
-    solar_angle = 20.0  #parameters['solar_angle']
-
-    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
-    cmap = cm.cividis
-
-    # If min and max appear to be reversed, then fix them
-    if (parameters['shadedminh'] > parameters['shadedmaxh']):
-        temp = parameters['shadedminh']
-        parameters['shadedminh'] = parameters['shadedmaxh']
-        parameters['shadedmaxh'] = temp
-    else:
-        parameters['shadedminh'] = parameters['shadedminh']
-        parameters['shadedmaxh'] = parameters['shadedmaxh']
-
-    # If no shadedmin/max parameters are read in from ctem.dat, determine the values from the data
-    if (parameters['shadedminhdefault'] == 1):
-        shadedminh = numpy.amin(DEM)
-    else:
-        shadedminh = parameters['shadedminh']
-    if (parameters['shadedmaxhdefault'] == 1):
-        shadedmaxh = numpy.amax(DEM)
-    else:
-        shadedmaxh = parameters['shadedmaxh']
-
-    dem_img = ls.shade(DEM, cmap=cmap,blend_mode=mode, fraction=1.0,
-                       vert_exag=ve, dx=pix, dy=pix,
-                       vmin=shadedminh, vmax=shadedmaxh)
-
-    #Generate image to put into an array
-    height = gridsize / dpi
-    width = gridsize / dpi
-    fig = matplotlib.pyplot.figure(figsize=(width, height), dpi=dpi)
-    ax = matplotlib.pyplot.axes([0, 0, 1, 1])
-    ax.imshow(dem_img,interpolation="nearest",vmin=0.0,vmax=1.0)
-    matplotlib.pyplot.axis('off')
-    # Save image to file
-    filename = parameters['workingdir'] + 'shaded' + os.sep + "shaded%06d.png" % parameters['ncount']
-    matplotlib.pyplot.savefig(filename,dpi=dpi,bbox_inches=0)
-    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 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
\ No newline at end of file
diff --git a/python/ctem/ctem/ctem_viewer_3d.py b/python/ctem/ctem/ctem_viewer_3d.py
deleted file mode 100644
index be967d10..00000000
--- a/python/ctem/ctem/ctem_viewer_3d.py
+++ /dev/null
@@ -1,92 +0,0 @@
-import numpy as np
-import os
-import open3d
-import ctem_io_readers
-import os
-
-
-class polysurface:
-    """A model of a self-gravitating small body"""
-    def __init__(self):
-        self.polyfile = ""
-        #used for Open3d module
-        self.mesh = open3d.geometry.TriangleMesh()
-        return
-    
-    def ctem2mesh(self, surface_file):
-    
-        # 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}
-
-        # Read ctem.in to initialize parameter values based on user input
-        ctem_io_readers.read_ctemin(parameters, notset)
-        # Read surface dem(shaded relief) and ejecta data files
-        dem_file = parameters['workingdir'] + surface_file
-        self.dem = ctem_io_readers.read_unformatted_binary(dem_file, parameters['gridsize'])
-        # Build mesh grid
-        s = parameters['gridsize']
-        pix = parameters['pix']
-        ygrid, xgrid = np.mgrid[-s / 2:s / 2, -s / 2:s / 2] * pix
-    
-        xvals = xgrid.flatten()
-        yvals = ygrid.flatten()
-        zvals = self.dem.flatten()
-        verts = np.array((xvals, yvals, zvals)).T
-        faces = np.empty([2 * (s - 1) ** 2, 3], dtype=np.int64)
-        for j in range(s - 1):
-            for i in range(s - 1):
-                idx = (s - 1) * j + i
-                # faces[idx,:] = [ j*s+i, j*s+i+1, (j+1)*s+i+1, (j+1)*s+i ]
-                faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i]
-                idx += (s - 1) ** 2
-                faces[idx, :] = [j * s + i + 1, (j + 1) * s + i + 1, (j + 1) * s + i]
-        self.mesh.vertices = open3d.utility.Vector3dVector(verts)
-        self.mesh.triangles = open3d.utility.Vector3iVector(faces)
-        self.mesh.compute_vertex_normals()
-        
-    
-if __name__ == '__main__':
-    import matplotlib.pyplot as plt
-
-    surf = polysurface()
-    surf.ctem2mesh(surface_file="surface_dem.dat")
-
-    zmax = np.amax(surf.dem)
-    zmin = np.amin(surf.dem)
-    cnorm = (surf.dem.flatten() - zmin) / (zmax - zmin)
-    cval = plt.cm.terrain(cnorm)[:,:3]
-    surf.mesh.vertex_colors = open3d.utility.Vector3dVector(cval)
-    open3d.visualization.draw_geometries([surf.mesh])
-    
-
diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
new file mode 100644
index 00000000..7dd8dfef
--- /dev/null
+++ b/python/ctem/ctem/driver.py
@@ -0,0 +1,262 @@
+import numpy as np
+import os
+import subprocess
+import shutil
+from ctem import io
+
+class Simulation:
+    """
+    This is a class that defines the basic CTEM simulation object. It initializes a dictionary of user and reads
+    it in from a file (default name is ctem.in). It also creates the directory structure needed to store simulation
+    files if necessary.
+    """
+    def __init__(self, param_file="ctem.in", isnew=True):
+        currentdir = os.getcwd()
+        self.user = {
+            'restart': None,
+            'runtype': None,
+            'popupconsole': None,
+            'saveshaded': None,
+            'saverego': None,
+            'savepres': None,
+            'savetruelist': None,
+            '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': os.path.join(currentdir, param_file),
+            'impfile': None,
+            'sfdcompare': None,
+            'sfdfile': None
+        }
+        self.user = io.read_user_input(self.user)
+        
+        self.output_filenames = {
+            'dat': 'ctem.dat',
+            'dem': 'surface_dem.dat',
+            'diam': 'surface_diam.dat',
+            'ejc': 'surface_ejc.dat',
+            'pos': 'surface_pos.dat',
+            'time': 'surface_time.dat',
+            'ocum': 'ocumulative.dat',
+            'odist': 'odistribution.dat',
+            'pdist': 'pdistribution.dat',
+            'tcum' : 'tcumulative.dat',
+            'tdist' : 'tdistribution.dat',
+            'impmass' : 'impactmass.dat',
+            'fracdone' : 'fracdone.dat',
+            'regodepth' : 'regolithdepth.dat',
+        }
+        
+        for k, v in self.output_filenames.items():
+            self.output_filenames[k] = os.path.join(currentdir, v)
+        
+        # Set up data arrays
+        self.seedarr = np.zeros(100, dtype=int)
+        self.seedarr[0] = self.user['seed']
+        self.odist = np.zeros([1, 6])
+        self.pdist = np.zeros([1, 6])
+        self.tdist = np.zeros([1, 6])
+        self.surface_dem = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
+        self.regolith = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
+        
+        if self.user['sfdcompare'] is not None:
+            # Read sfdcompare file
+            sfdfile = os.path.join(self.user['workingdir'], self.user['sfdcompare'])
+            self.ph1 = io.read_formatted_ascii(sfdfile, skip_lines=0)
+           
+
+        # Starting new or old run?
+        if (self.user['restart'].upper() == 'F' and isnew):
+            print('Starting a new run')
+
+            io.create_dir_structure(self.user)
+            # Delete any old output files
+            for k, v in self.output_filenames.items():
+                if os.path.isfile(v):
+                    os.remove(v)
+
+            # Scale the production function to the simulation domain
+            self.scale_production()
+            
+            if (self.user['runtype'].upper() == 'STATISTICAL'):
+                self.user['ncount'] = 1
+
+                # Write ctem.dat file
+                io.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+            else:
+                self.user['ncount'] = 0
+        else:
+            print('Continuing a previous run')
+            self.process_interval(isnew=False)
+
+        return
+    
+    def scale_production(self):
+        """
+        Scales the production function to the simulation domain and saves the scaled production function to a file that
+        will be read in by the main Fortran program.
+        """
+        
+        # Read production function file
+        impfile = os.path.join(self.user['workingdir'], self.user['impfile'])
+        prodfunction = io.read_formatted_ascii(impfile, skip_lines=0)
+
+        # Create impactor production population
+        area = (self.user['gridsize'] * self.user['pix']) ** 2
+        self.production = np.copy(prodfunction)
+        self.production[:, 1] = self.production[:, 1] * area * self.user['interval']
+
+        # Write the scaled production function to file
+        io.write_production(self.user, self.production)
+        
+        return
+    
+    def run(self):
+        """
+        Runs a complete simulation over all intervals specified in the input user. This method loops over all
+        intervals and process outputs into images, and redirect output files to their apporpriate folders.
+        
+        This method replaces most of the functionality of the original ctem_driver.py
+        """
+
+        print('Beginning loops')
+        while (self.user['ncount'] <= self.user['numintervals']):
+            if (self.user['ncount'] > 0):
+                # Move ctem.dat
+                self.redirect_outputs(['dat'], 'misc')
+                
+                #Execute Fortran code
+                self.compute_one_interval()
+                
+                # Process output files
+                self.process_interval()
+
+            # Display results
+            print(self.user['ncount'], '  Displaying results')
+
+            # Write surface dem, surface ejecta, shaded relief, and rplot data
+            io.image_dem(self.user, self.surface_dem)
+            if (self.user['saverego'].upper() == 'T'):
+                io.image_regolith(self.user, self.surface_ejc)
+            if (self.user['saveshaded'].upper() == 'T'):
+                io.image_shaded_relief(self.user, self.surface_dem)
+            if (self.user['savepres'].upper() == 'T'):
+                io.create_rplot(self.user, self.odist, self.pdist, self.tdist, self.ph1)
+
+            # Update ncount
+            self.user['ncount'] = self.user['ncount'] + 1
+
+            if ((self.user['runtype'].upper() == 'STATISTICAL') or (self.user['ncount'] == 1)):
+                self.user['restart'] = 'F'
+                self.user['curyear'] = 0.0
+                self.user['totalimpacts'] = 0
+                self.user['masstot'] = 0.0
+
+                # Delete tdistribution file, if it exists
+                tdist_file = self.user['workingdir'] + 'tdistribution.dat'
+                if os.path.isfile(tdist_file):
+                    os.remove(tdist_file)
+
+            else:
+                self.user['restart'] = 'T'
+
+            # Write ctem.dat file
+            io.write_datfile(self.user, self.output_filenames['dat'],self.seedarr)
+
+        return
+    
+    def compute_one_interval(self):
+        """
+        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')
+        with subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')],
+                              stdout=subprocess.PIPE,
+                              universal_newlines=True) as p:
+            for line in p.stdout:
+                print(line, end='')
+
+        return
+    
+    def process_interval(self, isnew=True):
+        """
+        Reads in all of the files generated by the Fortran code after one interval and redirects them to appropriate
+        folders for storage after appending the interval number to the filename.
+        """
+        # Read Fortran output
+        print(self.user['ncount'], '  Reading Fortran output')
+
+        # Read surface dem(shaded relief) and ejecta data files
+        self.surface_dem = io.read_unformatted_binary(self.output_filenames['dem'], self.user['gridsize'])
+        self.surface_ejc = io.read_unformatted_binary(self.output_filenames['ejc'], self.user['gridsize'])
+
+        # Read odistribution, tdistribution, and pdistribution files
+        self.odist = io.read_formatted_ascii(self.output_filenames['odist'], skip_lines=1)
+        self.tdist = io.read_formatted_ascii(self.output_filenames['tdist'], skip_lines=1)
+        self.pdist = io.read_formatted_ascii(self.output_filenames['pdist'], skip_lines=1)
+
+        # Read impact mass from file
+        impact_mass = io.read_impact_mass(self.output_filenames['impmass'])
+
+        # Read ctem.dat file
+        io.read_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+
+        # Save copy of crater distribution files
+        # Update user: mass, curyear, regolith properties
+        self.user['masstot'] = self.user['masstot'] + impact_mass
+
+        self.user['curyear'] = self.user['curyear'] + self.user['fracdone'] * self.user['interval']
+        template = "%(fracdone)9.6f %(curyear)19.12E\n"
+        with open(self.output_filenames['fracdone'], 'w') as fp_frac:
+            fp_frac.write(template % self.user)
+
+        reg_text = "%19.12E %19.12E %19.12E %19.12E\n" % (self.user['curyear'],
+                                                          np.mean(self.surface_ejc), np.amax(self.surface_ejc),
+                                                          np.amin(self.surface_ejc))
+        with open(self.output_filenames['regodepth'], 'w') as fp_reg:
+            fp_reg.write(reg_text)
+
+        if isnew:
+            self.redirect_outputs(['odist', 'ocum', 'pdist', 'tdist'], 'dist')
+            if (self.user['savetruelist'].upper() == 'T'):
+                self.redirect_outputs(['tcum'], 'dist')
+            self.redirect_outputs(['impmass'], 'misc')
+        return
+
+    def redirect_outputs(self, filekeys, foldername):
+        """
+        Copies a set of output files from the working directory into a subfolder.
+        Takes as input the dictionaries of user parameters and output file names, the dictionary keys to the file names
+        you wish to redirect, and the redirection destination folder name. The new name will be the key name + zero padded
+        interval number.
+
+        Example:
+            calling redirect_outputs(['odist','tdist'], 'dist') when user['ncount'] is 1
+            should do the following:
+            copy 'odistribution.dat' to 'dist/odist_000001.dat'
+            copy 'tdistribution.dat' to 'dist/tdist_000001.dat'
+        """
+        
+        for k in filekeys:
+            forig = self.output_filenames[k]
+            fdest = os.path.join(self.user['workingdir'], foldername, f"{k}_{self.user['ncount']:06d}.dat")
+            shutil.copy2(forig, fdest)
+
+if __name__ == '__main__':
+    sim = Simulation()
+    sim.run()
\ No newline at end of file
diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py
new file mode 100644
index 00000000..d0d383e5
--- /dev/null
+++ b/python/ctem/ctem/io.py
@@ -0,0 +1,382 @@
+import numpy as np
+import os
+import shutil
+from matplotlib.colors import LightSource
+import matplotlib.cm as cm
+import matplotlib.pyplot as plt
+
+# Set pixel scaling common for image writing, at 1 pixel/ array element
+dpi = 72.0
+
+def create_dir_structure(user):
+    # 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 = os.path.join(user['workingdir'], directory)
+        if not os.path.isdir(dir_test):
+            os.makedirs(dir_test)
+    
+    return
+
+
+def create_rplot(user, odist, pdist, tdist, ph1):
+    # Parameters: empirical saturation limit and dfrac
+    satlimit = 3.12636
+    dfrac = 2 ** (1. / 4) * 1.0e-3
+    
+    # Calculate geometric saturation
+    minx = (user['pix'] / 3.0) * 1.0e-3
+    maxx = 3 * user['pix'] * user['gridsize'] * 1.0e-3
+    geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]])
+    geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]])
+    
+    # Create distribution arrays without zeros for plotting on log scale
+    idx = np.nonzero(odist[:, 5])
+    odistnz = odist[idx]
+    odistnz = odistnz[:, [2, 5]]
+    
+    idx = np.nonzero(pdist[:, 5])
+    pdistnz = pdist[idx]
+    pdistnz = pdistnz[:, [2, 5]]
+    
+    idx = np.nonzero(tdist[:, 5])
+    tdistnz = tdist[idx]
+    tdistnz = tdistnz[:, [2, 5]]
+    
+    # Correct pdist
+    pdistnz[:, 1] = pdistnz[:, 1] * user['curyear'] / user['interval']
+    
+    # Create sdist bin factors, which contain one crater per bin
+    area = (user['gridsize'] * user['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 = np.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" % user['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 = os.path.join(user['workingdir'],'rplot',"rplot%06d.png" % user['ncount'])
+    height = user['gridsize'] / dpi
+    width = height
+    fig = plt.figure(figsize=(width, height), dpi=dpi)
+    
+    # Alter background color to be black, and change axis colors accordingly
+    plt.style.use('dark_background')
+    plt.rcParams['axes.prop_cycle']
+    
+    # Plot data
+    plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue')
+    plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white')
+    plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red')
+    
+    plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow')
+    plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow')
+    plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow')
+    
+    plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo')
+    
+    # Create plot labels
+    plt.title('Crater Distribution R-Plot', fontsize=22)
+    plt.xlim([minx, maxx])
+    plt.xscale('log')
+    plt.xlabel('Crater Diameter (km)', fontsize=18)
+    plt.ylim([5.0e-4, 5.0])
+    plt.yscale('log')
+    plt.ylabel('R Value', fontsize=18)
+    plt.text(1.0e-2, 1.0, timelabel, fontsize=18)
+    
+    plt.tick_params(axis='both', which='major', labelsize=14)
+    plt.tick_params(axis='both', which='minor', labelsize=12)
+    
+    plt.savefig(filename)
+    
+    return
+
+def image_dem(user, DEM):
+    dpi = 300.0  # 72.0
+    pix = user['pix']
+    gridsize = user['gridsize']
+    ve = 1.0
+    azimuth = 300.0  # user['azimuth']
+    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)
+    
+    # Generate image to put into an array
+    height = gridsize / dpi
+    width = gridsize / dpi
+    fig = plt.figure(figsize=(width, height), dpi=dpi)
+    ax = plt.axes([0, 0, 1, 1])
+    ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0)
+    plt.axis('off')
+    # Save image to file
+    filename = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount'])
+    plt.savefig(filename, dpi=dpi, bbox_inches=0)
+    
+    return
+
+
+def image_regolith(user, regolith):
+    # Create scaled regolith image
+    minref = user['pix'] * 1.0e-4
+    maxreg = np.amax(regolith)
+    minreg = np.amin(regolith)
+    if (minreg < minref): minreg = minref
+    if (maxreg < minref): maxreg = (minref + 1.0e3)
+    regolith_scaled = np.copy(regolith)
+    np.place(regolith_scaled, regolith_scaled < minref, minref)
+    regolith_scaled = 254.0 * (
+            (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg)))
+    
+    # Save image to file
+    filename = os.path.join(user['workingdir'],'rego', "rego%06d.png" % user['ncount'])
+    height = user['gridsize'] / dpi
+    width = height
+    fig = plt.figure(figsize=(width, height), dpi=dpi)
+    fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower')
+    plt.savefig(filename)
+    
+    return
+
+
+def image_shaded_relief(user, DEM):
+    dpi = 300.0  # 72.0
+    pix = user['pix']
+    gridsize = user['gridsize']
+    ve = 1.0
+    mode = 'overlay'
+    azimuth = 300.0  # user['azimuth']
+    solar_angle = 20.0  # user['solar_angle']
+    
+    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
+    cmap = cm.cividis
+    
+    # If min and max appear to be reversed, then fix them
+    if (user['shadedminh'] > user['shadedmaxh']):
+        temp = user['shadedminh']
+        user['shadedminh'] = user['shadedmaxh']
+        user['shadedmaxh'] = temp
+    else:
+        user['shadedminh'] = user['shadedminh']
+        user['shadedmaxh'] = user['shadedmaxh']
+    
+    # If no shadedmin/max user are read in from ctem.dat, determine the values from the data
+    if (user['shadedminhdefault'] == 1):
+        shadedminh = np.amin(DEM)
+    else:
+        shadedminh = user['shadedminh']
+    if (user['shadedmaxhdefault'] == 1):
+        shadedmaxh = np.amax(DEM)
+    else:
+        shadedmaxh = user['shadedmaxh']
+    
+    dem_img = ls.shade(DEM, cmap=cmap, blend_mode=mode, fraction=1.0,
+                       vert_exag=ve, dx=pix, dy=pix,
+                       vmin=shadedminh, vmax=shadedmaxh)
+    
+    # Generate image to put into an array
+    height = gridsize / dpi
+    width = gridsize / dpi
+    fig = plt.figure(figsize=(width, height), dpi=dpi)
+    ax = plt.axes([0, 0, 1, 1])
+    ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0)
+    plt.axis('off')
+    # Save image to file
+    filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount'])
+    plt.savefig(filename, dpi=dpi, bbox_inches=0)
+    return user
+
+
+def read_datfile(user, datfile, seedarr):
+    # Read and parse ctem.dat file
+    
+    # Read ctem.dat file
+    print('Reading input file ' + 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:
+        user['totalimpacts'] = real2float(fields[0])
+        user['ncount'] = int(fields[1])
+        user['curyear'] = real2float(fields[2])
+        user['restart'] = fields[3]
+        user['fracdone'] = real2float(fields[4])
+        user['masstot'] = real2float(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] = real2float(fields[0])
+        index += 1
+    
+    user['seedn'] = index - 1
+    
+    return
+
+
+def read_formatted_ascii(filename, skip_lines):
+    # Generalized ascii text reader
+    # For use with sfdcompare, production, odist, tdist, pdist data files
+    data = np.genfromtxt(filename, skip_header=skip_lines)
+    return data
+
+
+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 = real2float(fields[0])
+    else:
+        mass = 0
+    
+    return mass
+
+
+# Write production function to file production.dat
+# This file format does not exactly match that generated from IDL. Does it work?
+def read_user_input(user):
+    # Read and parse ctem.in file
+    inputfile = user['ctemfile']
+    
+    # Read ctem.in file
+    print('Reading input file ' + user['ctemfile'])
+    with open(inputfile, 'r') as fp:
+        lines = fp.readlines()
+    
+    # Process file text
+    for line in lines:
+        fields = line.split()
+        if len(fields) > 0:
+            if ('pix' == fields[0].lower()): user['pix'] = real2float(fields[1])
+            if ('gridsize' == fields[0].lower()): user['gridsize'] = int(fields[1])
+            if ('seed' == fields[0].lower()): user['seed'] = int(fields[1])
+            if ('sfdfile' == fields[0].lower()): user['sfdfile'] = os.path.join(user['workingdir'],fields[1])
+            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 ('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]
+            if ('saveshaded' == fields[0].lower()): user['saveshaded'] = fields[1]
+            if ('saverego' == fields[0].lower()): user['saverego'] = fields[1]
+            if ('savepres' == fields[0].lower()): user['savepres'] = fields[1]
+            if ('savetruelist' == fields[0].lower()): user['savetruelist'] = fields[1]
+            if ('runtype' == fields[0].lower()): user['runtype'] = fields[1]
+            if ('restart' == fields[0].lower()): user['restart'] = fields[1]
+            if ('shadedminh' == fields[0].lower()):
+                user['shadedminh'] = real2float(fields[1])
+                user['shadedminhdefault'] = 0
+            if ('shadedmaxh' == fields[0].lower()):
+                user['shadedmaxh'] = real2float(fields[1])
+                user['shadedmaxhdefault'] = 0
+    
+    # Test values for further processing
+    if (user['interval'] <= 0.0):
+        print('Invalid value for or missing variable INTERVAL in ' + inputfile)
+    if (user['numintervals'] <= 0):
+        print('Invalid value for or missing variable NUMINTERVALS in ' + inputfile)
+    if (user['pix'] <= 0.0):
+        print('Invalid value for or missing variable PIX in ' + inputfile)
+    if (user['gridsize'] <= 0):
+        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):
+        print('Invalid value for or missing variable POPUPCONSOLE in ' + inputfile)
+    if (user['saveshaded'] is None):
+        print('Invalid value for or missing variable SAVESHADED in ' + inputfile)
+    if (user['saverego'] is None):
+        print('Invalid value for or missing variable SAVEREGO in ' + inputfile)
+    if (user['savepres'] is None):
+        print('Invalid value for or missing variable SAVEPRES in ' + inputfile)
+    if (user['savetruelist'] is None):
+        print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile)
+    if (user['runtype'] is None):
+        print('Invalid value for or missing variable RUNTYPE in ' + inputfile)
+    if (user['restart'] is None):
+        print('Invalid value for or missing variable RESTART in ' + inputfile)
+    
+    return user
+
+
+def read_unformatted_binary(filename, gridsize):
+    # Read unformatted binary files created by Fortran
+    # For use with surface ejecta and surface dem data files
+    dt = np.float
+    data = np.fromfile(filename, dtype=dt)
+    data.shape = (gridsize, gridsize)
+    
+    return data
+
+
+def real2float(realstr):
+    """
+    Converts a Fortran-generated ASCII string of a real value into a np float type. Handles cases where double precision
+    numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E'
+
+    Parameters
+    ----------
+    realstr : string
+        Fortran-generated ASCII string of a real value.
+
+    Returns
+    -------
+     : float
+        The converted floating point value of the input string
+    """
+    return float(realstr.replace('d', 'E').replace('D', 'E'))
+
+
+def write_datfile(user, filename, seedarr):
+    # Write various user and random number seeds into ctem.dat file
+    fp = open(filename, 'w')
+    
+    template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n"
+    fp.write(template % user)
+    
+    # Write random number seeds to the file
+    for index in range(user['seedn']):
+        fp.write("%12d\n" % seedarr[index])
+    
+    fp.close()
+    
+    return
+
+
+def write_production(user, production):
+    filename = user['sfdfile']
+    np.savetxt(filename, production, fmt='%1.8e', delimiter='   ')
+    
+    return
+
diff --git a/python/ctem/ctem/viewer3d.py b/python/ctem/ctem/viewer3d.py
new file mode 100644
index 00000000..e3d66cb7
--- /dev/null
+++ b/python/ctem/ctem/viewer3d.py
@@ -0,0 +1,122 @@
+import numpy as np
+import open3d
+import ctem
+
+class Polysurface(ctem.Simulation):
+    """A model of a self-gravitating small body"""
+    def __init__(self):
+        ctem.Simulation.__init__(self, isnew=False) # Initialize the data structures, but doesn't start a new run
+        #used for Open3d module
+        self.mesh = open3d.geometry.TriangleMesh()
+        return
+    
+    def ctem2blockmesh(self):
+        # Build mesh grid
+        s = self.user['gridsize']
+        pix = self.user['pix']
+        ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix
+        y0, x0 = ygrid - pix/2, xgrid - pix/2
+        y1, x1 = ygrid - pix/2, xgrid + pix/2
+        y2, x2 = ygrid + pix/2, xgrid + pix/2
+        y3, x3 = ygrid + pix/2, xgrid - pix/2
+
+        xvals = np.append(
+                np.append(
+                np.append(x0.flatten(),
+                          x1.flatten()),
+                          x2.flatten()),
+                          x3.flatten())
+        yvals = np.append(
+                np.append(
+                np.append(y0.flatten(),
+                          y1.flatten()),
+                          y2.flatten()),
+                          y3.flatten())
+        zvals = np.append(
+                np.append(
+                np.append(self.dem.flatten(),
+                          self.dem.flatten()),
+                          self.dem.flatten()),
+                          self.dem.flatten())
+        verts = np.array((xvals, yvals, zvals)).T
+        nface_triangles = 10
+        faces = np.full([nface_triangles*s**2, 3], -1, dtype=np.int64)
+        for j in range(s):
+            for i in range(s):
+                i0 = s*j + i
+                i1 = i0 + s**2
+                i2 = i1 + s**2
+                i3 = i2 + s**2
+
+                fidx = np.arange(nface_triangles*i0,nface_triangles*(i0+1))
+                # Make the two top triangles
+                faces[fidx[0],:] = [i0, i1, i2]
+                faces[fidx[1],:] = [i0, i2, i3]
+                # Make the two west side triangles
+                if i > 0:
+                    faces[fidx[2],:] = [i0, i3, i2-1]
+                    faces[fidx[3],:] = [i0, i2-1, i1-1]
+                # Make the two south side triangles
+                if j > 0:
+                    faces[fidx[4],:] = [i1, i0, i3-s ]
+                    faces[fidx[5],:] = [i1, i3-s, i2-s]
+                # Make the two east side triangles
+                if i < (s - 1):
+                    faces[fidx[6],:] = [i2, i1, i0+1]
+                    faces[fidx[7],:] = [i2, i0+1, i3+1]
+                #Make the two north side triangles
+                if j < (s -1):
+                    faces[fidx[8],:] = [i3, i2, i1+s]
+                    faces[fidx[9],:] = [i3, i1+s, i0+s]
+        nz = faces[:,0] != -1
+        f2 = faces[nz,:]
+        self.mesh.vertices = open3d.utility.Vector3dVector(verts)
+        self.mesh.triangles = open3d.utility.Vector3iVector(f2)
+        self.mesh.compute_vertex_normals()
+        self.mesh.compute_triangle_normals()
+
+        return
+    
+    def ctem2trimesh(self):
+        # Build mesh grid
+        s = self.user['gridsize']
+        pix = self.user['pix']
+        ygrid, xgrid = np.mgrid[-s/2:s/2, -s/2:s/2] * pix
+
+        xvals = xgrid.flatten()
+        yvals = ygrid.flatten()
+        zvals = self.surface_dem.flatten()
+        verts = np.array((xvals, yvals, zvals)).T
+        faces = np.full([2 * (s-1)**2, 3], -1, dtype=np.int64)
+        for j in range(s - 1):
+            for i in range(s - 1):
+                idx = (s - 1) * j + i
+                faces[idx, :] = [j * s + i, j * s + i + 1, (j + 1) * s + i + 1]
+                idx += (s - 1) ** 2
+                faces[idx, :] = [(j + 1) * s + i + 1, (j + 1) * s + i, j * s + i ]
+        self.mesh.vertices = open3d.utility.Vector3dVector(verts)
+        self.mesh.triangles = open3d.utility.Vector3iVector(faces)
+        self.mesh.compute_vertex_normals()
+        self.mesh.compute_triangle_normals()
+
+        return
+    
+    def visualize(self):
+        vis = open3d.visualization.Visualizer()
+        vis.create_window()
+        vis.add_geometry(self.mesh)
+        opt = vis.get_render_option()
+        opt.background_color = np.asarray([0, 0, 0])
+
+        self.mesh.paint_uniform_color([0.5, 0.5, 0.5])
+        vis.run()
+        vis.destroy_window()
+
+
+if __name__ == '__main__':
+    sim = Polysurface()
+    sim.ctem2trimesh()
+    sim.visualize()
+
+    
+
diff --git a/python/ctem/tests/testio/ctem.in b/python/ctem/tests/testio/ctem.in
new file mode 100755
index 00000000..e57310dc
--- /dev/null
+++ b/python/ctem/tests/testio/ctem.in
@@ -0,0 +1,70 @@
+! CTEM Input file
+
+
+! Testing input. These are used to perform non-Monte Carlo tests.
+testflag     T                                 ! Set to T to create a single crater with user-defined impactor properties
+testimp      5.934e3 ! 93km crater            ! Diameter of test impactor (m)
+testvel      15.0e3                            ! Velocity of test crater (m/s)
+testang      90.0                              ! Impact angle of test crater (deg) - Default 90.0
+testxoffset  0.0e0                             ! x-axis offset of crater center from grid center (m) - Default 0.0
+testyoffset  0.0e0                             ! y-axis offset of crater center from grid center (m) - Default 0.0
+tallyonly    F                                 ! Tally the craters without generating any craters
+testtally    F   
+
+! IDL driver in uts
+interval     1.0
+numintervals 1                                  ! Total number of intervals
+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 
+saveshaded   T                                  ! 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)
+runtype      statistical                        !    single:      craters accumulate in successive intervals
+                                                !    statistical: surface is reset between intervals
+sfdcompare   FassettCounts.txt
+
+! CTEM required inputs
+seed         23790                    ! Random number generator seed
+gridsize     1000                          ! Size of grid in pixels
+numlayers    10                                 ! Number of perched layers
+pix          0.3468773099817e3                   ! Pixel size (m) - Copernicus DEM compariso
+mat          rock                               ! Material (rock or ice)
+! Bedrock scaling parameters
+mu_b         0.55e0                             ! Experimentally derived parameter for bedrock crater scaling law 
+kv_b         0.20e0                             ! Experimentally derived parameter for bedrock crater scaling law 
+trho_b       2250.0e0                           ! Target density (bedrock) (kg/m**3)
+ybar_b       0.0e6                               ! Bedrock strength (Pa)
+! Regolith scaling parameters
+mu_r         0.55e0                             ! Experimentally derived parameter for regolith crater scaling law 
+kv_r         0.20e0                             ! Experimentally derived parameter for regolith crater scaling law
+trho_r       2250.0e0                           ! Target density (regolith) (kg/m**3)
+ybar_r       0.00e6                              ! Regolith strength (Pa)
+! Body parameters
+gaccel       1.62e0                             ! Gravitational acceleration at target (m/s**2)
+trad         1737.35e3                          ! Target radius (m)
+prho         2500.0e0                           ! Projectile density (kg/m**3)
+sfdfile      production.dat                     ! Impactor SFD file (col 1: Dimp (m), col 2: ! impactors > D
+velfile      lunar-MBA-impactor-velocities.dat  ! Impactor velocity dist file
+
+! Seismic shaking input (required if seismic shaking is set to T, otherwise ignored)
+doseismic    F                                  ! Perform seismic shaking calculations with each impact - Default F
+
+! Optional inputF These have internally set default values that work reasonable well. Comment them out with
+deplimit        9e99                            ! Depth limit for craters (m) - Default is to ignore.
+maxcrat         1.00e0                          ! Fraction of gridsize that maximum crater can be - Default 1.0
+killatmaxcrater F                               ! Stop the run if a crater larger than the maximum is produced - Default F      
+basinimp        35.0e3                          ! Size of impactor to switch to lunar basin scaling law - Default is to ignore
+docollapse      T                               ! Do slope collapse - Default T
+dosoftening     F                               ! Do ejecta softening - Default T
+doangle         T                               ! Vary the impact angle. Set to F to have only vertical impacts - Default T
+Kd1             0.0001
+psi             2.000
+fe              2.00
+ejecta_truncation 2.0
+dorays  F
+superdomain F
+
+dorealistic T
+
diff --git a/python/ctem/tests/testio/testio.py b/python/ctem/tests/testio/testio.py
new file mode 100644
index 00000000..1212193e
--- /dev/null
+++ b/python/ctem/tests/testio/testio.py
@@ -0,0 +1,7 @@
+import ctem
+
+d = dir(ctem)
+
+sim = ctem.Simulation()
+for k, v in sim.user.items():
+    print(k, v)
\ No newline at end of file
diff --git a/python/ctem/tests/viewer3d/polytest.py b/python/ctem/tests/viewer3d/polytest.py
index 13cee521..c528374c 100644
--- a/python/ctem/tests/viewer3d/polytest.py
+++ b/python/ctem/tests/viewer3d/polytest.py
@@ -1,8 +1,9 @@
 import numpy as np
 from ctem import ctem_viewer_3d
-
-surf = ctem_viewer_3d.surf()
-surf.io_read_ply("216kleopatra-mesh.ply")
-v = np.asarray(surf.mesh)
-t = np.asarray(surf.mesh)
-surf.render()
\ No newline at end of file
+import open3d
+surf = ctem_viewer_3d.polysurface()
+#surf.io_read_ply("216kleopatra-mesh.ply")
+surf.mesh = open3d.geometry.TriangleMesh.create_box()
+v = np.asarray(surf.mesh.vertices)
+t = np.asarray(surf.mesh.triangles)
+surf.visualize()
\ No newline at end of file

From 300550f63826fd5eccd6ce94579ea387a733cc5f Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 14:07:49 -0400
Subject: [PATCH 02/19] Added matplotlib close methods after each image is
 generated so that they don't stay open and use up a bunch of memory

---
 python/ctem/ctem/io.py | 4 ++++
 1 file changed, 4 insertions(+)

diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py
index d0d383e5..144b419d 100644
--- a/python/ctem/ctem/io.py
+++ b/python/ctem/ctem/io.py
@@ -105,6 +105,7 @@ def create_rplot(user, odist, pdist, tdist, ph1):
     plt.tick_params(axis='both', which='minor', labelsize=12)
     
     plt.savefig(filename)
+    plt.close()
     
     return
 
@@ -129,6 +130,7 @@ def image_dem(user, DEM):
     # Save image to file
     filename = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount'])
     plt.savefig(filename, dpi=dpi, bbox_inches=0)
+    plt.close()
     
     return
 
@@ -201,6 +203,8 @@ def image_shaded_relief(user, DEM):
     # Save image to file
     filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount'])
     plt.savefig(filename, dpi=dpi, bbox_inches=0)
+    plt.close()
+
     return user
 
 

From 35fad7ba9f371f83c96be4727dd61443866cd899 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 14:11:19 -0400
Subject: [PATCH 03/19] Adjusted some parameters for generating the shaded and
 surf images to try to fix the 'flashing' problem that I started noticing on
 runs.

---
 python/ctem/ctem/driver.py |  6 +++---
 python/ctem/ctem/util.py   | 15 ++++++++++-----
 2 files changed, 13 insertions(+), 8 deletions(-)

diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index 401a29b1..89f05823 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -33,8 +33,8 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'gridsize': -1,
             'seed': 0,
             'maxcrat': 1.0,
-            'shadedminhdefault': 1,
-            'shadedmaxhdefault': 1,
+            'shadedminhdefault': 0,
+            'shadedmaxhdefault': 0,
             'shadedminh': 0.0,
             'shadedmaxh': 0.0,
             'workingdir': currentdir,
@@ -305,4 +305,4 @@ def cleanup(self):
 
 if __name__ == '__main__':
     sim = Simulation()
-    sim.run()
\ No newline at end of file
+    sim.run()
diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py
index 7d363f2f..488280dc 100644
--- a/python/ctem/ctem/util.py
+++ b/python/ctem/ctem/util.py
@@ -121,6 +121,7 @@ def create_rplot(user, odist, pdist, tdist, ph1):
     plt.tick_params(axis='both', which='minor', labelsize=12)
     
     plt.savefig(filename)
+    plt.close()
     
     return
 
@@ -130,7 +131,7 @@ def image_dem(user, DEM):
     gridsize = user['gridsize']
     ve = 1.0
     azimuth = 300.0  # user['azimuth']
-    solar_angle = 20.0  # user['solar_angle']
+    solar_angle = 15.0  # user['solar_angle']
     
     ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
     dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix)
@@ -145,6 +146,7 @@ def image_dem(user, DEM):
     # Save image to file
     filename = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount'])
     plt.savefig(filename, dpi=dpi, bbox_inches=0)
+    plt.close()
     
     return
 
@@ -168,6 +170,7 @@ def image_regolith(user, regolith):
     fig = plt.figure(figsize=(width, height), dpi=dpi)
     fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower')
     plt.savefig(filename)
+    plt.close()
     
     return
 
@@ -177,13 +180,13 @@ def image_shaded_relief(user, DEM):
     pix = user['pix']
     gridsize = user['gridsize']
     ve = 1.0
-    mode = 'overlay'
+    mode = 'hsv'
     azimuth = 300.0  # user['azimuth']
-    solar_angle = 20.0  # user['solar_angle']
+    solar_angle = 15.0  # user['solar_angle']
     
     ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
     cmap = cm.cividis
-    
+
     # If min and max appear to be reversed, then fix them
     if (user['shadedminh'] > user['shadedmaxh']):
         temp = user['shadedminh']
@@ -192,7 +195,7 @@ def image_shaded_relief(user, DEM):
     else:
         user['shadedminh'] = user['shadedminh']
         user['shadedmaxh'] = user['shadedmaxh']
-    
+
     # If no shadedmin/max user are read in from ctem.dat, determine the values from the data
     if (user['shadedminhdefault'] == 1):
         shadedminh = np.amin(DEM)
@@ -217,6 +220,8 @@ def image_shaded_relief(user, DEM):
     # Save image to file
     filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount'])
     plt.savefig(filename, dpi=dpi, bbox_inches=0)
+    plt.close()
+
     return user
 
 

From fa1f4b73e2c68af51e5983308ce035c73dd0b429 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 14:14:01 -0400
Subject: [PATCH 04/19] aded .idea to ignore because Pycharm

---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index f7ab64c6..1277aaf5 100644
--- a/.gitignore
+++ b/.gitignore
@@ -71,3 +71,5 @@ examples/global-lunar-bombardment/scale.ipynb
 *.png
 
 python/ctem/ctem.egg-info/
+
+.idea

From be5288ea771ce148782faec7f10c17f987484112 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 15:06:37 -0400
Subject: [PATCH 05/19] Modernized the driver on the global bombardment example

---
 .../craterproduction.py                       | 216 --------------
 .../global-lunar-bombardment/ctem_driver.py   | 260 +---------------
 .../ctem_io_readers.py                        | 146 ---------
 .../ctem_io_writers.py                        | 282 ------------------
 4 files changed, 4 insertions(+), 900 deletions(-)
 delete mode 100644 examples/global-lunar-bombardment/craterproduction.py
 delete mode 100644 examples/global-lunar-bombardment/ctem_io_readers.py
 delete mode 100644 examples/global-lunar-bombardment/ctem_io_writers.py

diff --git a/examples/global-lunar-bombardment/craterproduction.py b/examples/global-lunar-bombardment/craterproduction.py
deleted file mode 100644
index 4d81b687..00000000
--- a/examples/global-lunar-bombardment/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/global-lunar-bombardment/ctem_driver.py b/examples/global-lunar-bombardment/ctem_driver.py
index 7665387e..7903997c 100755
--- a/examples/global-lunar-bombardment/ctem_driver.py
+++ b/examples/global-lunar-bombardment/ctem_driver.py
@@ -1,256 +1,4 @@
-#!/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()
+#!/usr/bin/env python3
+import ctem
+sim = ctem.Simulation()
+sim.run()
diff --git a/examples/global-lunar-bombardment/ctem_io_readers.py b/examples/global-lunar-bombardment/ctem_io_readers.py
deleted file mode 100644
index c825fdc9..00000000
--- a/examples/global-lunar-bombardment/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/global-lunar-bombardment/ctem_io_writers.py b/examples/global-lunar-bombardment/ctem_io_writers.py
deleted file mode 100644
index e7a2192c..00000000
--- a/examples/global-lunar-bombardment/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

From 2153bd1c7392b639c89c310a681b9786a1776aa8 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 15:07:20 -0400
Subject: [PATCH 06/19] Wrote an override of the matplotlib LightSource class
 to remove the intensity scaling step which was causing the movies of the
 surface evolution to flicker.

---
 src/io/io_read_const.i90 | 54 ++++++++++++++++++++++++++++++++++++++++
 1 file changed, 54 insertions(+)
 create mode 100644 src/io/io_read_const.i90

diff --git a/src/io/io_read_const.i90 b/src/io/io_read_const.i90
new file mode 100644
index 00000000..12f34959
--- /dev/null
+++ b/src/io/io_read_const.i90
@@ -0,0 +1,54 @@
+# 1 "/Users/daminton/git/CTEM/src/io/io_read_const.f90"
+!**********************************************************************************************************************************
+!
+!  Unit Name   : io_read_const
+!  Unit Type   : subroutine
+!  Project     : CTEM
+!  Language    : Fortran 2003
+!
+!  Description : Read parameters passed between CTEM and the IDL CTEM driver routine
+!
+!  Input
+!    Arguments :
+!              :
+!
+!  Output
+!    Arguments : seed  : Current random number generator seeds
+!              : ncount : Current step count
+!              : totalimpacts  : Total number of impacts in the initial time step
+!              : nsteps : Number of steps to execute
+!              : restart : F = new run, T = restart old run
+!
+!  Notes       :
+!
+!**********************************************************************************************************************************
+
+subroutine io_read_const(totalimpacts,ncount,curyear,restart,fracdone,masstot,seedarr)
+   use module_globals
+   use module_io, EXCEPT_THIS_ONE => io_read_const
+   implicit none
+   
+! Arguments
+   integer(I8B),intent(out) :: totalimpacts
+   integer(I4B),intent(out) :: ncount
+   logical,intent(out) :: restart
+   real(DP),intent(out) :: curyear,fracdone,masstot
+   integer(I4B),dimension(:),intent(out) :: seedarr
+   
+! Internals
+   integer(I4B),parameter :: cfile=21
+   integer(I4B) :: ioerr,l
+
+   open(unit=cfile,file=DATFILE,status='old')
+   read(cfile,*) totalimpacts,ncount,curyear,restart,fracdone,masstot
+   do l=1,size(seedarr)
+      read(cfile,*,iostat=ioerr) seedarr(l)
+      if (ioerr/=0) seedarr(l)=abs(mod(((l-1)*181)*((seedarr(1)-83)*359),104729))
+   end do
+
+   close(cfile)
+
+   
+   return
+
+end subroutine io_read_const

From c6c266382eaef55eb461de52fd533d211434a453 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 15:07:20 -0400
Subject: [PATCH 07/19] Revert "Wrote an override of the matplotlib LightSource
 class to remove the intensity scaling step which was causing the movies of
 the surface evolution to flicker."

This reverts commit 2153bd1c7392b639c89c310a681b9786a1776aa8.
---
 src/io/io_read_const.i90 | 54 ----------------------------------------
 1 file changed, 54 deletions(-)
 delete mode 100644 src/io/io_read_const.i90

diff --git a/src/io/io_read_const.i90 b/src/io/io_read_const.i90
deleted file mode 100644
index 12f34959..00000000
--- a/src/io/io_read_const.i90
+++ /dev/null
@@ -1,54 +0,0 @@
-# 1 "/Users/daminton/git/CTEM/src/io/io_read_const.f90"
-!**********************************************************************************************************************************
-!
-!  Unit Name   : io_read_const
-!  Unit Type   : subroutine
-!  Project     : CTEM
-!  Language    : Fortran 2003
-!
-!  Description : Read parameters passed between CTEM and the IDL CTEM driver routine
-!
-!  Input
-!    Arguments :
-!              :
-!
-!  Output
-!    Arguments : seed  : Current random number generator seeds
-!              : ncount : Current step count
-!              : totalimpacts  : Total number of impacts in the initial time step
-!              : nsteps : Number of steps to execute
-!              : restart : F = new run, T = restart old run
-!
-!  Notes       :
-!
-!**********************************************************************************************************************************
-
-subroutine io_read_const(totalimpacts,ncount,curyear,restart,fracdone,masstot,seedarr)
-   use module_globals
-   use module_io, EXCEPT_THIS_ONE => io_read_const
-   implicit none
-   
-! Arguments
-   integer(I8B),intent(out) :: totalimpacts
-   integer(I4B),intent(out) :: ncount
-   logical,intent(out) :: restart
-   real(DP),intent(out) :: curyear,fracdone,masstot
-   integer(I4B),dimension(:),intent(out) :: seedarr
-   
-! Internals
-   integer(I4B),parameter :: cfile=21
-   integer(I4B) :: ioerr,l
-
-   open(unit=cfile,file=DATFILE,status='old')
-   read(cfile,*) totalimpacts,ncount,curyear,restart,fracdone,masstot
-   do l=1,size(seedarr)
-      read(cfile,*,iostat=ioerr) seedarr(l)
-      if (ioerr/=0) seedarr(l)=abs(mod(((l-1)*181)*((seedarr(1)-83)*359),104729))
-   end do
-
-   close(cfile)
-
-   
-   return
-
-end subroutine io_read_const

From 00f0fb6779036f00de8d32225e28682fae9da4f3 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 15:09:37 -0400
Subject: [PATCH 08/19] Auto stash before revert of "Wrote an override of the
 matplotlib LightSource class to remove the intensity scaling step which was
 causing the movies of the surface evolution to flicker."

added .i90 files to gitignore
---
 .gitignore | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/.gitignore b/.gitignore
index 1277aaf5..3688c36b 100644
--- a/.gitignore
+++ b/.gitignore
@@ -73,3 +73,5 @@ examples/global-lunar-bombardment/scale.ipynb
 python/ctem/ctem.egg-info/
 
 .idea
+
+*.i90

From 77b05cc86cf333181edffeb3917a849cca302c72 Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Fri, 21 Oct 2022 15:10:14 -0400
Subject: [PATCH 09/19] Wrote an override of the matplotlib LightSource class
 to remove the intensity scaling step which was causing the movies of the
 surface evolution to flicker.

---
 python/ctem/ctem/util.py | 69 ++++++++++++++++++++++++++++++++++++----
 1 file changed, 63 insertions(+), 6 deletions(-)

diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py
index 488280dc..c750e774 100644
--- a/python/ctem/ctem/util.py
+++ b/python/ctem/ctem/util.py
@@ -5,6 +5,62 @@
 import matplotlib.cm as cm
 import matplotlib.pyplot as plt
 
+
+class CTEMLightSource(LightSource):
+    """
+    Override one function in LightSource to prevent the contrast from being rescaled.
+    """
+    def shade_normals(self, normals, fraction=1.):
+        """
+        Calculate the illumination intensity for the normal vectors of a
+        surface using the defined azimuth and elevation for the light source.
+
+        Imagine an artificial sun placed at infinity in some azimuth and
+        elevation position illuminating our surface. The parts of the surface
+        that slope toward the sun should brighten while those sides facing away
+        should become darker.
+
+        Changes by David Minton: The matplotlib version of this rescales the intensity
+        in a way that causes the brightness level of the images to change between
+        simulation outputs. This makes movies of the surface evolution appear to flicker.
+
+        Parameters
+        ----------
+        fraction : number, optional
+            Increases or decreases the contrast of the hillshade.  Values
+            greater than one will cause intermediate values to move closer to
+            full illumination or shadow (and clipping any values that move
+            beyond 0 or 1). Note that this is not visually or mathematically
+            the same as vertical exaggeration.
+
+        Returns
+        -------
+        ndarray
+            A 2D array of illumination values between 0-1, where 0 is
+            completely in shadow and 1 is completely illuminated.
+        """
+
+        intensity = normals.dot(self.direction)
+
+        # Apply contrast stretch
+        imin, imax = intensity.min(), intensity.max()
+        intensity *= fraction
+
+        # Rescale to 0-1, keeping range before contrast stretch
+        # If constant slope, keep relative scaling (i.e. flat should be 0.5,
+        # fully occluded 0, etc.)
+        # if (imax - imin) > 1e-6:
+        #     # Strictly speaking, this is incorrect. Negative values should be
+        #     # clipped to 0 because they're fully occluded. However, rescaling
+        #     # in this manner is consistent with the previous implementation and
+        #     # visually appears better than a "hard" clip.
+        #     intensity -= imin
+        #     intensity /= (imax - imin)
+        intensity = np.clip(intensity, 0, 1)
+
+        return intensity
+
+
 # Set pixel scaling common for image writing, at 1 pixel/ array element
 dpi = 72.0
 
@@ -125,16 +181,17 @@ def create_rplot(user, odist, pdist, tdist, ph1):
     
     return
 
+
 def image_dem(user, DEM):
     dpi = 300.0  # 72.0
     pix = user['pix']
     gridsize = user['gridsize']
     ve = 1.0
     azimuth = 300.0  # user['azimuth']
-    solar_angle = 15.0  # user['solar_angle']
+    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)
+    ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle)
+    dem_img = ls.hillshade(DEM, vert_exag=ve, dx=pix, dy=pix, fraction=1)
     
     # Generate image to put into an array
     height = gridsize / dpi
@@ -180,11 +237,11 @@ def image_shaded_relief(user, DEM):
     pix = user['pix']
     gridsize = user['gridsize']
     ve = 1.0
-    mode = 'hsv'
+    mode = 'overlay'
     azimuth = 300.0  # user['azimuth']
-    solar_angle = 15.0  # user['solar_angle']
+    solar_angle = 20.0  # user['solar_angle']
     
-    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
+    ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle)
     cmap = cm.cividis
 
     # If min and max appear to be reversed, then fix them

From 2ca5a86e62bbc502f82a8fb719f95de7e21bb69c Mon Sep 17 00:00:00 2001
From: David Minton <daminton@purdue.edu>
Date: Sat, 19 Nov 2022 11:16:37 -0500
Subject: [PATCH 10/19] Fixed issue where progress bar makes a new line each
 time it updates when running via Python

---
 .gitignore                               |   2 +
 INSTALL                                  | 368 -----------------------
 python/ctem/ctem/driver.py               |  21 +-
 src/crater/crater_populate.f90           |  12 +-
 src/crater/crater_subpixel_diffusion.f90 |  12 +-
 src/io/io_updatePbar.f90                 |  22 +-
 6 files changed, 39 insertions(+), 398 deletions(-)
 delete mode 100644 INSTALL

diff --git a/.gitignore b/.gitignore
index 3688c36b..0eeffdd9 100644
--- a/.gitignore
+++ b/.gitignore
@@ -75,3 +75,5 @@ python/ctem/ctem.egg-info/
 .idea
 
 *.i90
+
+INSTALL
diff --git a/INSTALL b/INSTALL
deleted file mode 100644
index 8865734f..00000000
--- a/INSTALL
+++ /dev/null
@@ -1,368 +0,0 @@
-Installation Instructions
-*************************
-
-   Copyright (C) 1994-1996, 1999-2002, 2004-2016 Free Software
-Foundation, Inc.
-
-   Copying and distribution of this file, with or without modification,
-are permitted in any medium without royalty provided the copyright
-notice and this notice are preserved.  This file is offered as-is,
-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
-instructions specific to this package.  Some packages provide this
-'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
-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
-you can run in the future to recreate the current configuration, and a
-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.
-
-   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
-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
-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 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.
-
-     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.
-
-  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
-     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
-     privileges.
-
-  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
-     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
-     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
-     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
-     distcheck', which can by used by developers to test that all other
-     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'
-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:
-
-     ./configure CC=c99 CFLAGS=-g LIBS=-lposix
-
-   *Note Defining Variables::, for more details.
-
-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
-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.
-
-   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
-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
-this:
-
-     ./configure CC="gcc -arch i386 -arch x86_64 -arch ppc -arch ppc64" \
-                 CXX="g++ -arch i386 -arch x86_64 -arch ppc -arch ppc64" \
-                 CPP="gcc -E" CXXCPP="g++ -E"
-
-   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.
-
-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
-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
-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
-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
-both of the following shortcuts of passing variable assignments to the
-'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
-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
-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.
-
-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
-package recognizes.
-
-   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.
-
-   Some packages offer the ability to configure how verbose the
-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
---disable-silent-rules' sets the default to verbose, which can be
-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
-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.
-
-   On OSF/1 a.k.a. Tru64, some versions of the default C compiler cannot
-parse its '<wchar.h>' 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"
-
-and if that doesn't work, try
-
-     ./configure CC="cc -nodtk"
-
-   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'.
-
-   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
-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
-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:
-
-     CPU-COMPANY-SYSTEM
-
-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
-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
-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'.
-
-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.
-
-Defining Variables
-==================
-
-   Variables not defined in a site shell script can be set in the
-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:
-
-     ./configure CC=/usr/local2/bin/gcc
-
-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:
-
-     CONFIG_SHELL=/bin/bash ./configure CONFIG_SHELL=/bin/bash
-
-'configure' Invocation
-======================
-
-   '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=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.
-
-'--version'
-'-V'
-     Print the version of Autoconf used to generate the 'configure'
-     script, and exit.
-
-'--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
-     disable caching.
-
-'--config-cache'
-'-C'
-     Alias for '--cache-file=config.cache'.
-
-'--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
-     messages will still be shown).
-
-'--srcdir=DIR'
-     Look for the package's source code in directory DIR.  Usually
-     '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.
-
-'--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.
diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index 89f05823..9ac146f6 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -180,17 +180,20 @@ def compute_one_interval(self):
         # 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')],
+            with subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')],
                               stdout=subprocess.PIPE,
                               stderr=subprocess.PIPE,
-                              universal_newlines=True)
-            for line in p.stdout:
-                print(line, end='')
-            res = p.communicate()
-            if p.returncode != 0:
-                for line in res[1]:
-                    print(line, end='')
-                raise Exception ("CTEM Failure")
+                              universal_newlines=True) as p:
+                for line in p.stdout:
+                    if "%" in line:
+                        print(line.replace('\n','\r'), end='')
+                    else:
+                        print(line,end='')
+                res = p.communicate()
+                if p.returncode != 0:
+                    for line in res[1]:
+                        print(line, end='')
+                    raise Exception ("CTEM Failure")
         except:
             print("Error executing main CTEM program")
             sys.exit()
diff --git a/src/crater/crater_populate.f90 b/src/crater/crater_populate.f90
index e60b68ce..e5e8f4db 100644
--- a/src/crater/crater_populate.f90
+++ b/src/crater/crater_populate.f90
@@ -158,7 +158,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
       !if in quasiMC mode: check to see if it's time for a real crater
       if (user%doquasimc) then
          if ((user%rctime > timestamp_old) .and. (user%rctime < crater%timestamp)) then
-            write(*,*) "Real crater at ", crater%timestamp
+            write(message,*) "Real @ ", crater%timestamp
+            call io_updatePbar(message)
             user%testflag = .true.
             user%testimp = rclist(1, rccount)
             user%testvel = rclist(2, rccount)
@@ -169,8 +170,10 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
       end if
       ! generate random crater
       call crater_generate(user,crater,domain,prod,production_list,vdist,surf)
-      if (user%testflag) write(*,*) 'Dcrat = ',crater%fcrat
-      if (user%testflag) write(*,*) 'Dtrans = ',crater%rad*2
+      if (user%testflag) then
+         write(message,'("Dc=",F8.1," Dt=",F8.1)') crater%fcrat, crater%rad*2
+         call io_updatePbar(message)
+      end if
       if (crater%fcrat > domain%biggest_crater) then ! End the run if the crater is too big
          if ( user%testflag .eqv. .false. ) then
             if (user%killatmaxcrater) then 
@@ -188,7 +191,8 @@ subroutine crater_populate(user,surf,crater,domain,prod,production_list,vdist,nt
             user%testflag = .false.
             rccount = rccount + 1
             if (rccount > domain%rcnum) then
-               write(*,*) "Real crater list complete."
+               write(message,*) "Real crater list complete."
+               call io_updatePbar(message)
                user%rctime = 1e30
             else
                user%rctime = rclist(6,rccount)
diff --git a/src/crater/crater_subpixel_diffusion.f90 b/src/crater/crater_subpixel_diffusion.f90
index f3c3e27e..095b53ec 100644
--- a/src/crater/crater_subpixel_diffusion.f90
+++ b/src/crater/crater_subpixel_diffusion.f90
@@ -207,12 +207,12 @@ subroutine crater_subpixel_diffusion(user,surf,nflux,domain,finterval,kdiffin)
    kdiff(user%gridsize + 1,:) = kdiff(1,:)
    kdiff(:,user%gridsize + 1) = kdiff(:,1)  
 
-   write(*,*)
-   write(*,*) 'avgkdiff = ',sum(kdiff) / (user%gridsize + 2)**2 / finterval
-   write(*,*)
-   open(unit=55,file='avgkdiff.dat',status='unknown',position='append')
-   write(55,*) sum(kdiff) / (user%gridsize + 2)**2 / finterval
-   close(55)
+   ! write(*,*)
+   ! write(*,*) 'avgkdiff = ',sum(kdiff) / (user%gridsize + 2)**2 / finterval
+   ! write(*,*)
+   ! open(unit=55,file='avgkdiff.dat',status='unknown',position='append')
+   ! write(55,*) sum(kdiff) / (user%gridsize + 2)**2 / finterval
+   ! close(55)
 
 
    call util_diffusion_solver(user,surf,user%gridsize + 2,indarray,kdiff,cumulative_elchange,maxhits)
diff --git a/src/io/io_updatePbar.f90 b/src/io/io_updatePbar.f90
index 7ba679d9..eab51800 100644
--- a/src/io/io_updatePbar.f90
+++ b/src/io/io_updatePbar.f90
@@ -42,17 +42,17 @@ subroutine io_updatePbar(message)
 do k = 1,endpoint !- 1 
    bar(6+k:6+k)=pbarchar(k:k)
 end do
-!select case(flip)
-!case(1)
-!   bar(6+endpoint:6+endpoint)="/"
-!case(2)
-!   bar(6+endpoint:6+endpoint)="-"
-!case(3)
-!   bar(6+endpoint:6+endpoint)="\"
-!case(4)
-!   bar(6+endpoint:6+endpoint)="|"
-!end select
-!flip = flip + 1
+select case(flip)
+case(1)
+  bar(6+endpoint:6+endpoint)="/"
+case(2)
+  bar(6+endpoint:6+endpoint)="-"
+case(3)
+  bar(6+endpoint:6+endpoint)="\"
+case(4)
+  bar(6+endpoint:6+endpoint)="|"
+end select
+flip = flip + 1
 if (flip > 4) flip = 1
 ! print the progress bar.  
 write(fmtlabel,'("(A1,A",I2.2,",1X,A",I2.2,",$)")') PBARSIZE+7,MESSAGESIZE

From aa5b2d847085093222d5e413e9fe52a159a6e2ad Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Wed, 1 Feb 2023 13:46:36 -0500
Subject: [PATCH 11/19] Replaced the util.py function with the updated one from
 the allocatable_arrays branch

---
 python/ctem/ctem/util.py | 191 +++++++++++++++++----------------------
 1 file changed, 85 insertions(+), 106 deletions(-)

diff --git a/python/ctem/ctem/util.py b/python/ctem/ctem/util.py
index c750e774..dab91401 100644
--- a/python/ctem/ctem/util.py
+++ b/python/ctem/ctem/util.py
@@ -4,8 +4,12 @@
 from matplotlib.colors import LightSource
 import matplotlib.cm as cm
 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
 class CTEMLightSource(LightSource):
     """
     Override one function in LightSource to prevent the contrast from being rescaled.
@@ -61,16 +65,14 @@ def shade_normals(self, normals, fraction=1.):
         return intensity
 
 
-# Set pixel scaling common for image writing, at 1 pixel/ array element
-dpi = 72.0
 
 # These are directories that are created by CTEM in order to store intermediate results
-directories = ['dist', 'misc', 'rego', 'rplot', 'surf', 'shaded']
 
-def create_dir_structure(user):
+def create_dir_structure(user, directories):
     """
     Create directories for various output files if they do not already exist
     """
+
     for dirname in directories:
         dirpath = os.path.join(user['workingdir'], dirname)
         if not os.path.isdir(dirpath):
@@ -79,7 +81,7 @@ def create_dir_structure(user):
     
     return
 
-def destroy_dir_structure(user):
+def destroy_dir_structure(user, directories):
     """
     Deletes directories generated by create_dir_structure
     """
@@ -91,97 +93,6 @@ def destroy_dir_structure(user):
 
     return
 
-
-def create_rplot(user, odist, pdist, tdist, ph1):
-    # Parameters: empirical saturation limit and dfrac
-    satlimit = 3.12636
-    dfrac = 2 ** (1. / 4) * 1.0e-3
-    
-    # Calculate geometric saturation
-    minx = (user['pix'] / 3.0) * 1.0e-3
-    maxx = 3 * user['pix'] * user['gridsize'] * 1.0e-3
-    geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]])
-    geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]])
-    
-    # Create distribution arrays without zeros for plotting on log scale
-    idx = np.nonzero(odist[:, 5])
-    odistnz = odist[idx]
-    odistnz = odistnz[:, [2, 5]]
-    
-    idx = np.nonzero(pdist[:, 5])
-    pdistnz = pdist[idx]
-    pdistnz = pdistnz[:, [2, 5]]
-    
-    idx = np.nonzero(tdist[:, 5])
-    tdistnz = tdist[idx]
-    tdistnz = tdistnz[:, [2, 5]]
-    
-    # Correct pdist
-    pdistnz[:, 1] = pdistnz[:, 1] * user['curyear'] / user['interval']
-    
-    # Create sdist bin factors, which contain one crater per bin
-    area = (user['gridsize'] * user['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 = np.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" % user['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 = os.path.join(user['workingdir'],'rplot',"rplot%06d.png" % user['ncount'])
-    height = user['gridsize'] / dpi
-    width = height
-    fig = plt.figure(figsize=(width, height), dpi=dpi)
-    
-    # Alter background color to be black, and change axis colors accordingly
-    plt.style.use('dark_background')
-    plt.rcParams['axes.prop_cycle']
-    
-    # Plot data
-    plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue')
-    plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white')
-    plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red')
-    
-    plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    
-    plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo')
-    
-    # Create plot labels
-    plt.title('Crater Distribution R-Plot', fontsize=22)
-    plt.xlim([minx, maxx])
-    plt.xscale('log')
-    plt.xlabel('Crater Diameter (km)', fontsize=18)
-    plt.ylim([5.0e-4, 5.0])
-    plt.yscale('log')
-    plt.ylabel('R Value', fontsize=18)
-    plt.text(1.0e-2, 1.0, timelabel, fontsize=18)
-    
-    plt.tick_params(axis='both', which='major', labelsize=14)
-    plt.tick_params(axis='both', which='minor', labelsize=12)
-    
-    plt.savefig(filename)
-    plt.close()
-    
-    return
-
-
 def image_dem(user, DEM):
     dpi = 300.0  # 72.0
     pix = user['pix']
@@ -243,7 +154,7 @@ def image_shaded_relief(user, DEM):
     
     ls = CTEMLightSource(azdeg=azimuth, altdeg=solar_angle)
     cmap = cm.cividis
-
+    
     # If min and max appear to be reversed, then fix them
     if (user['shadedminh'] > user['shadedmaxh']):
         temp = user['shadedminh']
@@ -252,7 +163,7 @@ def image_shaded_relief(user, DEM):
     else:
         user['shadedminh'] = user['shadedminh']
         user['shadedmaxh'] = user['shadedmaxh']
-
+    
     # If no shadedmin/max user are read in from ctem.dat, determine the values from the data
     if (user['shadedminhdefault'] == 1):
         shadedminh = np.amin(DEM)
@@ -278,7 +189,7 @@ def image_shaded_relief(user, DEM):
     filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount'])
     plt.savefig(filename, dpi=dpi, bbox_inches=0)
     plt.close()
-
+    
     return user
 
 
@@ -363,6 +274,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]
@@ -372,6 +284,7 @@ def read_user_input(user):
             if ('savetruelist' == fields[0].lower()): user['savetruelist'] = fields[1]
             if ('runtype' == fields[0].lower()): user['runtype'] = fields[1]
             if ('restart' == fields[0].lower()): user['restart'] = fields[1]
+            if ('quasimc' == fields[0].lower()): user['quasimc'] = fields[1]
             if ('shadedminh' == fields[0].lower()):
                 user['shadedminh'] = real2float(fields[1])
                 user['shadedminhdefault'] = 0
@@ -390,8 +303,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):
@@ -412,16 +323,40 @@ 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('<i4')
+    elif kind == 'I8B':
+        dt = np.dtypye('<i8')
     data = np.fromfile(filename, dtype=dt)
     data.shape = (gridsize, gridsize)
     
     return data
 
 
+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('<i4')
+    elif kind == 'I8B':
+        dt = np.dtypye('<i8')
+    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):
+                data[i, j] = f.read_reals(dt)
+    return data
+
+
 def real2float(realstr):
     """
     Converts a Fortran-generated ASCII string of a real value into a np float type. Handles cases where double precision
@@ -439,6 +374,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
@@ -456,9 +418,26 @@ 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
 
+
+def write_realcraters(filename, realcraters):
+    """Writes file of real craters for use in quasi-MC runs"""
+
+    np.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t')
+
+    return
+
+def write_temp_input(user, filename):
+    """Makes changes to a temporary input file for use when generating craterlist.dat for quasimc runs"""
+
+    sed('testflag', 'testflag T!', filename)
+    sed('testimp', f'testimp {user["pix"]*1e-3} !', filename) # Make a tiny test crater. We don't care about the crater itself, just that we run CTEM once to get all of the converted files
+    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 52160b5d92b47d4068de8a2e6290b3ae3350823e Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 10 Mar 2023 16:23:19 -0500
Subject: [PATCH 12/19] Pulled driver.py from the allocatable_arrays branch

---
 python/ctem/ctem/driver.py | 140 +++++++++++++++++++++++++++----------
 1 file changed, 102 insertions(+), 38 deletions(-)

diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index 9ac146f6..24f16c9c 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -4,6 +4,12 @@
 import shutil
 from ctem import util
 import sys
+import pandas
+from ctem import craterproduction
+from scipy.interpolate import interp1d
+from ctem import __file__ as _pyfile
+from pathlib import Path
+import warnings
 
 class Simulation:
     """
@@ -19,7 +25,6 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'popupconsole': None,
             'saveshaded': None,
             'saverego': None,
-            'savepres': None,
             'savetruelist': None,
             'seedn': 1,
             'totalimpacts': 0,
@@ -33,18 +38,31 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'gridsize': -1,
             'seed': 0,
             'maxcrat': 1.0,
-            'shadedminhdefault': 0,
-            'shadedmaxhdefault': 0,
+            'shadedminhdefault': 1,
+            'shadedmaxhdefault': 1,
             '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,
         }
+
+        # Get the location of the CTEM executable
+        self.ctem_executable = Path(_pyfile).parent.parent.parent.parent / "build" / "src" / "CTEM"
+        if not self.ctem_executable.exists():
+            print(f"CTEM driver not found at {self.ctem_executable}. Trying current directory.")
+            self.ctem_executable = Path(currentdir) / "CTEM"
+            if not self.ctem_executable.exists():
+                warnings.warn(f"Cannot find the CTEM driver {str(self.ctem_executable)}", stacklevel=2)
+                self.ctem_executable = None
+
+
         self.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 = {
@@ -54,6 +72,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',
@@ -65,12 +88,22 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'ejmax' : 'ejecta_table_max.dat',
             'ejmin' : 'ejecta_table_min.dat',
             'testprof' : 'testprofile.dat',
-            'craterscale' : 'craterscale.dat'
+            'craterscale' : 'craterscale.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)
-        
+
+        self.directories = ['dist', 'misc', 'surf']
+        if self.user['saveshaded'].upper() == 'T':
+            self.directories.append('shaded')
+        if self.user['saverego'].upper() == 'T':
+            self.directories.append('rego')
+            
         # Set up data arrays
         self.seedarr = np.zeros(100, dtype=int)
         self.seedarr[0] = self.user['seed']
@@ -79,6 +112,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
         self.tdist = np.zeros([1, 6])
         self.surface_dem = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
         self.surface_ejc = np.zeros([self.user['gridsize'], self.user['gridsize']], dtype=float)
+        self.ph1 = None
         
         if self.user['sfdcompare'] is not None:
             # Read sfdcompare file
@@ -92,7 +126,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
             if (self.user['restart'].upper() == 'F'):
                 print('Starting a new run')
 
-                util.create_dir_structure(self.user)
+                util.create_dir_structure(self.user, self.directories)
                 # Delete any old output files
                 for k, v in self.output_filenames.items():
                     if os.path.isfile(v):
@@ -100,7 +134,42 @@ def __init__(self, param_file="ctem.in", isnew=True):
 
                 # Scale the production function to the simulation domain
                 self.scale_production()
-                
+
+                # Setup Quasi-MC run
+
+                if (self.user['quasimc'] == 'T'):
+
+                    #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)
+                    tempfile = os.path.join(currentdir, 'temp.in')
+
+                    # Generate craterlist.dat
+                    shutil.copy2(self.user['ctemfile'], tempfile )
+                    
+                    #Write a temporary input file to generate the necessary quasimc files
+                    util.write_temp_input(self.user, tempfile)
+                    util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+                    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+')
+                    df['log(Dc)'] = np.log(df['Dcrat(m)'])
+                    df['log(Di)'] = np.log(df['#Dimp(m)'])
+                    xnew = df['log(Dc)'].values
+                    ynew = df['log(Di)'].values
+                    interp = interp1d(xnew, ynew, fill_value='extrapolate')
+                    rclist[:,0] = np.exp(interp(np.log(rclist[:,0])))
+    
+                    #Convert age in Ga to "interval time"
+                    rclist[:,5] = (self.user['interval'] * self.user['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon')
+                    rclist = rclist[rclist[:,5].argsort()]
+
+                    #Export to dat file
+                    util.write_realcraters(self.output_filenames['craterlist'], rclist)
+
                 util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
             else:
                 print('Continuing a previous run')
@@ -124,7 +193,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
     
@@ -143,12 +212,12 @@ 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()
 
-            # Process the o utput files
+            # Process the output files
             self.process_output()
 
             # Update ncount
@@ -173,27 +242,28 @@ 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:
-            with subprocess.Popen([os.path.join(self.user['workingdir'], 'CTEM')],
+            p = subprocess.Popen([os.path.join(self.user['workingdir'], self.ctem_executable), ctemin],
                               stdout=subprocess.PIPE,
                               stderr=subprocess.PIPE,
-                              universal_newlines=True) as p:
-                for line in p.stdout:
-                    if "%" in line:
-                        print(line.replace('\n','\r'), end='')
-                    else:
-                        print(line,end='')
-                res = p.communicate()
-                if p.returncode != 0:
-                    for line in res[1]:
-                        print(line, end='')
-                    raise Exception ("CTEM Failure")
+                              universal_newlines=True,
+                               shell=False)
+            for line in p.stdout:
+                if "%" in line:
+                    print(line.replace('\n','\r'), end='')
+                else:
+                    print(line,end='')
+            res = p.communicate()
+            if p.returncode != 0:
+                for line in res[1]:
+                    print(line, end='')
+                raise Exception ("CTEM Failure")
         except:
             print("Error executing main CTEM program")
             sys.exit()
@@ -222,6 +292,7 @@ def read_output(self):
 
         # Read ctem.dat file
         util.read_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+
         
     def process_output(self):
         """
@@ -231,7 +302,7 @@ def process_output(self):
         # Display results
         print(self.user['ncount'], '  Generating surface images and plots')
 
-        # Write surface dem, surface ejecta, shaded relief, and rplot data
+        # Write surface dem, surface ejecta and shaded relief data
         util.image_dem(self.user, self.surface_dem)
         if (self.user['saverego'].upper() == 'T'):
             util.image_regolith(self.user, self.surface_ejc)
@@ -239,8 +310,6 @@ def process_output(self):
             util.image_shaded_relief(self.user, self.surface_dem)
         
         if self.user['ncount'] > 0: # These aren't available yet from the initial conditions
-            if (self.user['savepres'].upper() == 'T'):
-                util.create_rplot(self.user, self.odist, self.pdist, self.tdist, self.ph1)
             
             # Save copy of crater distribution files
             # Update user: mass, curyear, regolith properties
@@ -261,6 +330,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')
             
 
 
@@ -289,23 +360,16 @@ def cleanup(self):
         """
         # This is a list of files generated by the main Fortran program
         print("Deleting all files generated by CTEM")
-        util.destroy_dir_structure(self.user)
+        util.destroy_dir_structure(self.user, self.directories)
         for key, filename in self.output_filenames.items():
             print(f"Deleting file {filename}")
             try:
                 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
 
 if __name__ == '__main__':
     sim = Simulation()
-    sim.run()
+    sim.run()
\ No newline at end of file

From b71f385f2118dd85aa33a22ec42f8e38f5e55f1e Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 10 Mar 2023 16:24:40 -0500
Subject: [PATCH 13/19] Got rid of quasimc stuff from master where it isn't
 implemented yet

---
 python/ctem/ctem/driver.py | 36 ------------------------------------
 1 file changed, 36 deletions(-)

diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index 24f16c9c..bbc9a616 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -46,7 +46,6 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'ctemfile': param_file,
             'impfile': None,
             'sfdcompare': None,
-            'quasimc': None,
             'sfdfile' : None,
             'realcraterlist': None,
         }
@@ -135,41 +134,6 @@ def __init__(self, param_file="ctem.in", isnew=True):
                 # Scale the production function to the simulation domain
                 self.scale_production()
 
-                # Setup Quasi-MC run
-
-                if (self.user['quasimc'] == 'T'):
-
-                    #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)
-                    tempfile = os.path.join(currentdir, 'temp.in')
-
-                    # Generate craterlist.dat
-                    shutil.copy2(self.user['ctemfile'], tempfile )
-                    
-                    #Write a temporary input file to generate the necessary quasimc files
-                    util.write_temp_input(self.user, tempfile)
-                    util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
-                    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+')
-                    df['log(Dc)'] = np.log(df['Dcrat(m)'])
-                    df['log(Di)'] = np.log(df['#Dimp(m)'])
-                    xnew = df['log(Dc)'].values
-                    ynew = df['log(Di)'].values
-                    interp = interp1d(xnew, ynew, fill_value='extrapolate')
-                    rclist[:,0] = np.exp(interp(np.log(rclist[:,0])))
-    
-                    #Convert age in Ga to "interval time"
-                    rclist[:,5] = (self.user['interval'] * self.user['numintervals']) - craterproduction.Tscale(rclist[:,5], 'NPF_Moon')
-                    rclist = rclist[rclist[:,5].argsort()]
-
-                    #Export to dat file
-                    util.write_realcraters(self.output_filenames['craterlist'], rclist)
-
                 util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
             else:
                 print('Continuing a previous run')

From d7c25e02b106b3c3421be9fa96411e444fbc6ed3 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 10 Mar 2023 16:27:04 -0500
Subject: [PATCH 14/19] Got rid of extra imports from other branch

---
 python/ctem/ctem/driver.py | 5 -----
 1 file changed, 5 deletions(-)

diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index bbc9a616..dd5ae9e0 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -4,12 +4,7 @@
 import shutil
 from ctem import util
 import sys
-import pandas
-from ctem import craterproduction
-from scipy.interpolate import interp1d
 from ctem import __file__ as _pyfile
-from pathlib import Path
-import warnings
 
 class Simulation:
     """

From c79d6d338affe8480cbd50ff06276ca2f1b1c668 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 10 Mar 2023 16:27:49 -0500
Subject: [PATCH 15/19] Added back stuff we actually need

---
 python/ctem/ctem/driver.py | 2 ++
 1 file changed, 2 insertions(+)

diff --git a/python/ctem/ctem/driver.py b/python/ctem/ctem/driver.py
index dd5ae9e0..63deb939 100644
--- a/python/ctem/ctem/driver.py
+++ b/python/ctem/ctem/driver.py
@@ -5,6 +5,8 @@
 from ctem import util
 import sys
 from ctem import __file__ as _pyfile
+from pathlib import Path
+import warnings
 
 class Simulation:
     """

From 2d6a990c1ffa19c100a956e6601bcaa6855b2ca5 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Fri, 10 Mar 2023 16:31:39 -0500
Subject: [PATCH 16/19] Got rid of quasimc stuff

---
 examples/morphology_test_cases/global/ctem_io_writers.py | 7 -------
 python/ctem/ctem/util.py                                 | 8 --------
 2 files changed, 15 deletions(-)

diff --git a/examples/morphology_test_cases/global/ctem_io_writers.py b/examples/morphology_test_cases/global/ctem_io_writers.py
index e7a2192c..d07d7d54 100644
--- a/examples/morphology_test_cases/global/ctem_io_writers.py
+++ b/examples/morphology_test_cases/global/ctem_io_writers.py
@@ -273,10 +273,3 @@ def create_rplot(parameters,odist,pdist,tdist,ph1):
    
     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/python/ctem/ctem/util.py b/python/ctem/ctem/util.py
index dab91401..61314f4c 100644
--- a/python/ctem/ctem/util.py
+++ b/python/ctem/ctem/util.py
@@ -274,7 +274,6 @@ 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]
@@ -424,13 +423,6 @@ def write_production(filename, production):
     return
 
 
-def write_realcraters(filename, realcraters):
-    """Writes file of real craters for use in quasi-MC runs"""
-
-    np.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t')
-
-    return
-
 def write_temp_input(user, filename):
     """Makes changes to a temporary input file for use when generating craterlist.dat for quasimc runs"""
 

From 29518f67f2ab965d8ef516ef0c3e5b1acd265bcb Mon Sep 17 00:00:00 2001
From: MintoDA1 <51412913+MintoDA1@users.noreply.github.com>
Date: Mon, 11 Sep 2023 15:20:08 -0400
Subject: [PATCH 17/19] Replaced tabs with spaces

---
 src/realistic/realistic_crater_topography.f90 | 10 +++++-----
 1 file changed, 5 insertions(+), 5 deletions(-)

diff --git a/src/realistic/realistic_crater_topography.f90 b/src/realistic/realistic_crater_topography.f90
index 81cb1dab..6768fa6b 100644
--- a/src/realistic/realistic_crater_topography.f90
+++ b/src/realistic/realistic_crater_topography.f90
@@ -164,11 +164,11 @@ subroutine make_realistic_crater(user,surf,crater,deltaMtot)
    psd_elevation%num_sine=psd_elevation%num_vertices/2
    psd_elevation%input_x_max=PI*psd_elevation%diameter_in_km
    psd_elevation%diameter_in_km_trans=20.0_DP
-   psd_elevation%bp2_x_k_b_sigma=[0.0417732075218209_DP ,	2.48787783435609e-18_DP ,	0.496231041304548_DP ,	0.00210183162260986_DP ,	0.793427517984220_DP ,	0.308862248697524_DP]
-   psd_elevation%bp2_y_k_b_sigma=[0.0593460236785858_DP ,	0.408350328939867_DP ,	1.53480069842064_DP ,	-0.00201958329120022_DP ,	1.63566246833559_DP ,	0.809757199553671_DP]
-   psd_elevation%bp3_y_k_b_sigma=[0.0561324763037810_DP ,	2.46101909781936_DP ,	0.537503435665947_DP ,	-0.00292059238533839_DP ,	3.64208047160175_DP ,	0.462627466875698_DP]
-   psd_elevation%bp4_y_k_b_sigma=[0.0159696581958859_DP ,	3.25122353414224_DP ,	0.660958440041299_DP ,	0.00638368649794204_DP ,	3.44294296810112_DP ,	0.769783029462181_DP]
-   psd_elevation%slope12_k_b_sigma=[0.0279866462789263_DP ,	2.53280320608848_DP ,	0.492438940196540_DP ,	0.000799959919780047_DP ,	3.07653693327141_DP ,	0.400693610967074_DP]
+   psd_elevation%bp2_x_k_b_sigma=[0.0417732075218209_DP , 2.48787783435609e-18_DP , 0.496231041304548_DP , 0.00210183162260986_DP , 0.793427517984220_DP , 0.308862248697524_DP]
+   psd_elevation%bp2_y_k_b_sigma=[0.0593460236785858_DP , 0.408350328939867_DP , 1.53480069842064_DP , -0.00201958329120022_DP , 1.63566246833559_DP , 0.809757199553671_DP]
+   psd_elevation%bp3_y_k_b_sigma=[0.0561324763037810_DP , 2.46101909781936_DP , 0.537503435665947_DP , -0.00292059238533839_DP , 3.64208047160175_DP , 0.462627466875698_DP]
+   psd_elevation%bp4_y_k_b_sigma=[0.0159696581958859_DP , 3.25122353414224_DP , 0.660958440041299_DP , 0.00638368649794204_DP , 3.44294296810112_DP , 0.769783029462181_DP]
+   psd_elevation%slope12_k_b_sigma=[0.0279866462789263_DP , 2.53280320608848_DP , 0.492438940196540_DP , 0.000799959919780047_DP , 3.07653693327141_DP , 0.400693610967074_DP]
    psd_elevation%misfit_k_b_sigma=0.000000_DP
    call Calculate_am_wl_phase_from_diameter(psd_elevation,amplitude_rim_elevation,wavelength_rim_elevation,phase_rim_elevation)
 

From bf43045f23b7690c88f2a0cb7fa93b0980dd3760 Mon Sep 17 00:00:00 2001
From: David A Minton <daminton@purdue.edu>
Date: Mon, 11 Sep 2023 15:26:14 -0400
Subject: [PATCH 18/19] Removed old file from pre-restructure

---
 python/ctem/ctem/io.py | 386 -----------------------------------------
 1 file changed, 386 deletions(-)
 delete mode 100644 python/ctem/ctem/io.py

diff --git a/python/ctem/ctem/io.py b/python/ctem/ctem/io.py
deleted file mode 100644
index 144b419d..00000000
--- a/python/ctem/ctem/io.py
+++ /dev/null
@@ -1,386 +0,0 @@
-import numpy as np
-import os
-import shutil
-from matplotlib.colors import LightSource
-import matplotlib.cm as cm
-import matplotlib.pyplot as plt
-
-# Set pixel scaling common for image writing, at 1 pixel/ array element
-dpi = 72.0
-
-def create_dir_structure(user):
-    # 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 = os.path.join(user['workingdir'], directory)
-        if not os.path.isdir(dir_test):
-            os.makedirs(dir_test)
-    
-    return
-
-
-def create_rplot(user, odist, pdist, tdist, ph1):
-    # Parameters: empirical saturation limit and dfrac
-    satlimit = 3.12636
-    dfrac = 2 ** (1. / 4) * 1.0e-3
-    
-    # Calculate geometric saturation
-    minx = (user['pix'] / 3.0) * 1.0e-3
-    maxx = 3 * user['pix'] * user['gridsize'] * 1.0e-3
-    geomem = np.array([[minx, satlimit / 20.0], [maxx, satlimit / 20.0]])
-    geomep = np.array([[minx, satlimit / 10.0], [maxx, satlimit / 10.0]])
-    
-    # Create distribution arrays without zeros for plotting on log scale
-    idx = np.nonzero(odist[:, 5])
-    odistnz = odist[idx]
-    odistnz = odistnz[:, [2, 5]]
-    
-    idx = np.nonzero(pdist[:, 5])
-    pdistnz = pdist[idx]
-    pdistnz = pdistnz[:, [2, 5]]
-    
-    idx = np.nonzero(tdist[:, 5])
-    tdistnz = tdist[idx]
-    tdistnz = tdistnz[:, [2, 5]]
-    
-    # Correct pdist
-    pdistnz[:, 1] = pdistnz[:, 1] * user['curyear'] / user['interval']
-    
-    # Create sdist bin factors, which contain one crater per bin
-    area = (user['gridsize'] * user['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 = np.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" % user['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 = os.path.join(user['workingdir'],'rplot',"rplot%06d.png" % user['ncount'])
-    height = user['gridsize'] / dpi
-    width = height
-    fig = plt.figure(figsize=(width, height), dpi=dpi)
-    
-    # Alter background color to be black, and change axis colors accordingly
-    plt.style.use('dark_background')
-    plt.rcParams['axes.prop_cycle']
-    
-    # Plot data
-    plt.plot(odistnz[:, 0] * 1.0e-3, odistnz[:, 1], linewidth=3.0, color='blue')
-    plt.plot(pdistnz[:, 0] * 1.0e-3, pdistnz[:, 1], linewidth=2.0, linestyle='dashdot', color='white')
-    plt.plot(tdistnz[:, 0] * 1.0e-3, tdistnz[:, 1], linewidth=2.0, color='red')
-    
-    plt.plot(geomem[:, 0], geomem[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    plt.plot(geomep[:, 0], geomep[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    plt.plot(sdist[:, 0], sdist[:, 1], linewidth=2.0, linestyle=':', color='yellow')
-    
-    plt.plot(ph1[:, 0] * dfrac, ph1[:, 1], 'wo')
-    
-    # Create plot labels
-    plt.title('Crater Distribution R-Plot', fontsize=22)
-    plt.xlim([minx, maxx])
-    plt.xscale('log')
-    plt.xlabel('Crater Diameter (km)', fontsize=18)
-    plt.ylim([5.0e-4, 5.0])
-    plt.yscale('log')
-    plt.ylabel('R Value', fontsize=18)
-    plt.text(1.0e-2, 1.0, timelabel, fontsize=18)
-    
-    plt.tick_params(axis='both', which='major', labelsize=14)
-    plt.tick_params(axis='both', which='minor', labelsize=12)
-    
-    plt.savefig(filename)
-    plt.close()
-    
-    return
-
-def image_dem(user, DEM):
-    dpi = 300.0  # 72.0
-    pix = user['pix']
-    gridsize = user['gridsize']
-    ve = 1.0
-    azimuth = 300.0  # user['azimuth']
-    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)
-    
-    # Generate image to put into an array
-    height = gridsize / dpi
-    width = gridsize / dpi
-    fig = plt.figure(figsize=(width, height), dpi=dpi)
-    ax = plt.axes([0, 0, 1, 1])
-    ax.imshow(dem_img, interpolation="nearest", cmap='gray', vmin=0.0, vmax=1.0)
-    plt.axis('off')
-    # Save image to file
-    filename = os.path.join(user['workingdir'],'surf',"surf%06d.png" % user['ncount'])
-    plt.savefig(filename, dpi=dpi, bbox_inches=0)
-    plt.close()
-    
-    return
-
-
-def image_regolith(user, regolith):
-    # Create scaled regolith image
-    minref = user['pix'] * 1.0e-4
-    maxreg = np.amax(regolith)
-    minreg = np.amin(regolith)
-    if (minreg < minref): minreg = minref
-    if (maxreg < minref): maxreg = (minref + 1.0e3)
-    regolith_scaled = np.copy(regolith)
-    np.place(regolith_scaled, regolith_scaled < minref, minref)
-    regolith_scaled = 254.0 * (
-            (np.log(regolith_scaled) - np.log(minreg)) / (np.log(maxreg) - np.log(minreg)))
-    
-    # Save image to file
-    filename = os.path.join(user['workingdir'],'rego', "rego%06d.png" % user['ncount'])
-    height = user['gridsize'] / dpi
-    width = height
-    fig = plt.figure(figsize=(width, height), dpi=dpi)
-    fig.figimage(regolith_scaled, cmap=cm.nipy_spectral, origin='lower')
-    plt.savefig(filename)
-    
-    return
-
-
-def image_shaded_relief(user, DEM):
-    dpi = 300.0  # 72.0
-    pix = user['pix']
-    gridsize = user['gridsize']
-    ve = 1.0
-    mode = 'overlay'
-    azimuth = 300.0  # user['azimuth']
-    solar_angle = 20.0  # user['solar_angle']
-    
-    ls = LightSource(azdeg=azimuth, altdeg=solar_angle)
-    cmap = cm.cividis
-    
-    # If min and max appear to be reversed, then fix them
-    if (user['shadedminh'] > user['shadedmaxh']):
-        temp = user['shadedminh']
-        user['shadedminh'] = user['shadedmaxh']
-        user['shadedmaxh'] = temp
-    else:
-        user['shadedminh'] = user['shadedminh']
-        user['shadedmaxh'] = user['shadedmaxh']
-    
-    # If no shadedmin/max user are read in from ctem.dat, determine the values from the data
-    if (user['shadedminhdefault'] == 1):
-        shadedminh = np.amin(DEM)
-    else:
-        shadedminh = user['shadedminh']
-    if (user['shadedmaxhdefault'] == 1):
-        shadedmaxh = np.amax(DEM)
-    else:
-        shadedmaxh = user['shadedmaxh']
-    
-    dem_img = ls.shade(DEM, cmap=cmap, blend_mode=mode, fraction=1.0,
-                       vert_exag=ve, dx=pix, dy=pix,
-                       vmin=shadedminh, vmax=shadedmaxh)
-    
-    # Generate image to put into an array
-    height = gridsize / dpi
-    width = gridsize / dpi
-    fig = plt.figure(figsize=(width, height), dpi=dpi)
-    ax = plt.axes([0, 0, 1, 1])
-    ax.imshow(dem_img, interpolation="nearest", vmin=0.0, vmax=1.0)
-    plt.axis('off')
-    # Save image to file
-    filename = os.path.join(user['workingdir'],'shaded',"shaded%06d.png" % user['ncount'])
-    plt.savefig(filename, dpi=dpi, bbox_inches=0)
-    plt.close()
-
-    return user
-
-
-def read_datfile(user, datfile, seedarr):
-    # Read and parse ctem.dat file
-    
-    # Read ctem.dat file
-    print('Reading input file ' + 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:
-        user['totalimpacts'] = real2float(fields[0])
-        user['ncount'] = int(fields[1])
-        user['curyear'] = real2float(fields[2])
-        user['restart'] = fields[3]
-        user['fracdone'] = real2float(fields[4])
-        user['masstot'] = real2float(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] = real2float(fields[0])
-        index += 1
-    
-    user['seedn'] = index - 1
-    
-    return
-
-
-def read_formatted_ascii(filename, skip_lines):
-    # Generalized ascii text reader
-    # For use with sfdcompare, production, odist, tdist, pdist data files
-    data = np.genfromtxt(filename, skip_header=skip_lines)
-    return data
-
-
-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 = real2float(fields[0])
-    else:
-        mass = 0
-    
-    return mass
-
-
-# Write production function to file production.dat
-# This file format does not exactly match that generated from IDL. Does it work?
-def read_user_input(user):
-    # Read and parse ctem.in file
-    inputfile = user['ctemfile']
-    
-    # Read ctem.in file
-    print('Reading input file ' + user['ctemfile'])
-    with open(inputfile, 'r') as fp:
-        lines = fp.readlines()
-    
-    # Process file text
-    for line in lines:
-        fields = line.split()
-        if len(fields) > 0:
-            if ('pix' == fields[0].lower()): user['pix'] = real2float(fields[1])
-            if ('gridsize' == fields[0].lower()): user['gridsize'] = int(fields[1])
-            if ('seed' == fields[0].lower()): user['seed'] = int(fields[1])
-            if ('sfdfile' == fields[0].lower()): user['sfdfile'] = os.path.join(user['workingdir'],fields[1])
-            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 ('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]
-            if ('saveshaded' == fields[0].lower()): user['saveshaded'] = fields[1]
-            if ('saverego' == fields[0].lower()): user['saverego'] = fields[1]
-            if ('savepres' == fields[0].lower()): user['savepres'] = fields[1]
-            if ('savetruelist' == fields[0].lower()): user['savetruelist'] = fields[1]
-            if ('runtype' == fields[0].lower()): user['runtype'] = fields[1]
-            if ('restart' == fields[0].lower()): user['restart'] = fields[1]
-            if ('shadedminh' == fields[0].lower()):
-                user['shadedminh'] = real2float(fields[1])
-                user['shadedminhdefault'] = 0
-            if ('shadedmaxh' == fields[0].lower()):
-                user['shadedmaxh'] = real2float(fields[1])
-                user['shadedmaxhdefault'] = 0
-    
-    # Test values for further processing
-    if (user['interval'] <= 0.0):
-        print('Invalid value for or missing variable INTERVAL in ' + inputfile)
-    if (user['numintervals'] <= 0):
-        print('Invalid value for or missing variable NUMINTERVALS in ' + inputfile)
-    if (user['pix'] <= 0.0):
-        print('Invalid value for or missing variable PIX in ' + inputfile)
-    if (user['gridsize'] <= 0):
-        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):
-        print('Invalid value for or missing variable POPUPCONSOLE in ' + inputfile)
-    if (user['saveshaded'] is None):
-        print('Invalid value for or missing variable SAVESHADED in ' + inputfile)
-    if (user['saverego'] is None):
-        print('Invalid value for or missing variable SAVEREGO in ' + inputfile)
-    if (user['savepres'] is None):
-        print('Invalid value for or missing variable SAVEPRES in ' + inputfile)
-    if (user['savetruelist'] is None):
-        print('Invalid value for or missing variable SAVETRUELIST in ' + inputfile)
-    if (user['runtype'] is None):
-        print('Invalid value for or missing variable RUNTYPE in ' + inputfile)
-    if (user['restart'] is None):
-        print('Invalid value for or missing variable RESTART in ' + inputfile)
-    
-    return user
-
-
-def read_unformatted_binary(filename, gridsize):
-    # Read unformatted binary files created by Fortran
-    # For use with surface ejecta and surface dem data files
-    dt = np.float
-    data = np.fromfile(filename, dtype=dt)
-    data.shape = (gridsize, gridsize)
-    
-    return data
-
-
-def real2float(realstr):
-    """
-    Converts a Fortran-generated ASCII string of a real value into a np float type. Handles cases where double precision
-    numbers in exponential notation use 'd' or 'D' instead of 'e' or 'E'
-
-    Parameters
-    ----------
-    realstr : string
-        Fortran-generated ASCII string of a real value.
-
-    Returns
-    -------
-     : float
-        The converted floating point value of the input string
-    """
-    return float(realstr.replace('d', 'E').replace('D', 'E'))
-
-
-def write_datfile(user, filename, seedarr):
-    # Write various user and random number seeds into ctem.dat file
-    fp = open(filename, 'w')
-    
-    template = "%(totalimpacts)17d %(ncount)12d %(curyear)19.12E %(restart)s %(fracdone)9.6f %(masstot)19.12E\n"
-    fp.write(template % user)
-    
-    # Write random number seeds to the file
-    for index in range(user['seedn']):
-        fp.write("%12d\n" % seedarr[index])
-    
-    fp.close()
-    
-    return
-
-
-def write_production(user, production):
-    filename = user['sfdfile']
-    np.savetxt(filename, production, fmt='%1.8e', delimiter='   ')
-    
-    return
-

From 53c4246cb7f499d94c09ffbe5d19e24ddd4a71fd Mon Sep 17 00:00:00 2001
From: MintoDA1 <51412913+MintoDA1@users.noreply.github.com>
Date: Tue, 12 Sep 2023 11:42:21 -0400
Subject: [PATCH 19/19] Changed relative path to CTEM executable

---
 .gitignore          | 3 +--
 CMakeLists.txt      | 2 +-
 ctem/ctem/driver.py | 2 +-
 3 files changed, 3 insertions(+), 4 deletions(-)

diff --git a/.gitignore b/.gitignore
index ab7aae65..733dc353 100644
--- a/.gitignore
+++ b/.gitignore
@@ -6,6 +6,5 @@
 !README.md
 !version.txt
 !src/**.f90
-!ctem/setup.py
-!ctem/ctem/*.py
+!ctem/**.py
 !cmake/Modules/*.cmake
\ No newline at end of file
diff --git a/CMakeLists.txt b/CMakeLists.txt
index fd259936..eb4fdb32 100644
--- a/CMakeLists.txt
+++ b/CMakeLists.txt
@@ -29,7 +29,7 @@ ENDIF ()
 OPTION(USE_COARRAY "Use Coarray Fortran for parallelization of test particles" OFF)
 OPTION(USE_OPENMP "Use OpenMP for parallelization" ON)
 OPTION(USE_SIMD "Use SIMD vectorization" ON)
-OPTION(BUILD_SHARED_LIBS "Build using shared libraries" OFF)
+OPTION(BUILD_SHARED_LIBS "Build using shared libraries" ON)
 
 # The following section is modified from Numpy f2py documentation
 IF(PROJECT_SOURCE_DIR STREQUAL PROJECT_BINARY_DIR)
diff --git a/ctem/ctem/driver.py b/ctem/ctem/driver.py
index 0c9bd611..0249c19d 100644
--- a/ctem/ctem/driver.py
+++ b/ctem/ctem/driver.py
@@ -52,7 +52,7 @@ def __init__(self, param_file="ctem.in", isnew=True):
         }
 
         # Get the location of the CTEM executable
-        self.ctem_executable = Path(_pyfile).parent.parent.parent.parent / "build" / "src" / "CTEM"
+        self.ctem_executable = Path(_pyfile).parent.parent.parent / "bin" / "CTEM"
         if not self.ctem_executable.exists():
             print(f"CTEM driver not found at {self.ctem_executable}. Trying current directory.")
             self.ctem_executable = Path(currentdir) / "CTEM"