tried and failed to add quasimc routines to driver.py

examples/global-lunar-bombardment/driver.py

@@ -0,0 +1,3 @@
+import ctem
+sim = ctem.Simulation()
+sim.run()

python/ctem/ctem/craterproduction.py

@@ -0,0 +1,216 @@
+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()

python/ctem/ctem/driver.py

@@ -4,6 +4,9 @@
 import shutil
 from ctem import util
 import sys
+import pandas
+from ctem import craterproduction
+from scipy.interpolate import interp1d
 
 class Simulation:
     """
@@ -41,7 +44,9 @@ def __init__(self, param_file="ctem.in", isnew=True):
             'ctemfile': os.path.join(currentdir, param_file),
             'impfile': None,
             'sfdcompare': None,
-            'sfdfile': None
+            'sfdfile': None,
+            'quasimc': None,
+            'realcraterlist': None
         }
         self.user = util.read_user_input(self.user)
 
@@ -65,7 +70,8 @@ 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'
         }
 
         for k, v in self.output_filenames.items():
@@ -100,6 +106,33 @@ 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")
+                    craterlistfile = self.user['workingdir'] + self.user['realcraterlist']
+                    rclist = util.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)'] = 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(user, rclist)
+
+
 
                 util.write_datfile(self.user, self.output_filenames['dat'], self.seedarr)
             else:
@@ -219,6 +252,9 @@ def read_output(self):
 
         # Read ctem.dat file
         util.read_datfile(self.user, self.output_filenames['dat'], self.seedarr)
+
+        # Read craterlist.dat file
+        self.realcraterlist = util.read_formatted_ascii(self.output_filenames['craterlist'], skip_lines=1)
 
     def process_output(self):
         """

python/ctem/ctem/util.py

@@ -400,3 +400,11 @@ def write_production(user, production):
 
     return
 
+
+def write_realcraters(user, realcraters):
+    """Writes file of real craters for use in quasi-MC runs"""
+
+    filename = user['workingdir'] + user['craterlist']
+    np.savetxt(filename, realcraters, fmt='%1.8e', delimiter='\t')
+
+    return

src/Makefile.am

@@ -10,7 +10,7 @@ AM_FCFLAGS =  $(IPRODUCTION)
 #ifort debug flags
 
 #gfortran optimized flags
-#AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
+AM_FCFLAGS =  -O3 -fopenmp -ffree-form  -g -fbounds-check -fbacktrace
 #gfortran debug flags
 #AM_FCFLAGS =  -O0 -g -fopenmp -fbounds-check -Wall -Warray-bounds -Warray-temporaries -Wimplicit-interface -ffree-form -fsanitize-address-use-after-scope -fstack-check  -fsanitize=bounds-strict -fsanitize=undefined -fsanitize=signed-integer-overflow -fsanitize=object-size -fstack-protector-all