Replaced the util.py function with the updated one from the allocatab…

…le_arrays branch

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