Restructured python files and fixed bugs to make the code more compat…

…ible with ifx

ctem/NPF.py

@@ -0,0 +1,104 @@
+import numpy as np
+# The Neukum production function
+# Ivanov, Neukum, and Hartmann (2001) SSR v. 96 pp. 55-86
+def N1lunar(T):
+	return 5.44e-14 * (np.exp(6.93*T) - 1) + 8.38e-4 * T
+
+
+def Nlunar(D):
+    # Lunar crater SFD
+    aL00 = -3.0876
+    aL01 = -3.557528
+    aL02 =  0.781027
+    aL03 =  1.021521
+    aL04 = -0.156012
+    aL05 = -0.444058
+    aL06 =  0.019977
+    aL07 =  0.086850
+    aL08 = -0.005874
+    aL09 = -0.006809
+    aL10 =  8.25e-4
+    aL11 =  5.54e-5
+
+    return np.where((D > 0.01) & (D < 1000.0),
+        10**(aL00
+        + aL01 * np.log10(D) ** 1
+        + aL02 * np.log10(D) ** 2
+        + aL03 * np.log10(D) ** 3
+        + aL04 * np.log10(D) ** 4
+        + aL05 * np.log10(D) ** 5
+        + aL06 * np.log10(D) ** 6
+        + aL07 * np.log10(D) ** 7
+        + aL08 * np.log10(D) ** 8
+        + aL09 * np.log10(D) ** 9
+        + aL10 * np.log10(D) ** 10
+        + aL11 * np.log10(D) ** 11),float('nan'))
+
+def Rproj(D):
+    #Projectile SFD
+    aP00 =  0.0
+    aP01 = -1.375458
+    aP02 =  1.272521e-1
+    aP03 = -1.282166
+    aP04 = -3.074558e-1
+    aP05 =  4.149280e-1
+    aP06 =  1.910668e-1
+    aP07 = -4.260980e-2
+    aP08 = -3.976305e-2
+    aP09 = -3.180179e-3
+    aP10 =  2.799369e-3
+    aP11 =  6.892223e-4
+    aP12 =  2.614385e-6
+    aP13 = -1.416178e-5
+    aP14 = -1.191124e-6
+
+    return np.where((D > 1e-4) & (D < 300),
+        10**(aP00
+        + aP01 * np.log10(D)**1
+        + aP02 * np.log10(D)**2
+        + aP03 * np.log10(D)**3
+        + aP04 * np.log10(D)**4
+        + aP05 * np.log10(D)**5
+        + aP06 * np.log10(D)**6
+        + aP07 * np.log10(D)**7
+        + aP08 * np.log10(D)**8
+        + aP09 * np.log10(D)**9
+        + aP10 * np.log10(D)**10
+        + aP11 * np.log10(D)**11
+        + aP12 * np.log10(D)**12
+        + aP13 * np.log10(D)**13
+        + aP14 * np.log10(D)**14),float('nan'))
+
+def N1mars(T):
+    # Mars crater SFD
+    # Ivanov (2001) SSR v. 96 pp. 87-104
+    return 2.68e-14 * (np.exp(6.93 * T) - 1) + 4.13e-4 * T
+
+def Nmars(D):
+    aM00 = -3.384
+    aM01 = -3.197
+    aM02 = 1.257
+    aM03 = 0.7915
+    aM04 = -0.4861
+    aM05 = -0.3630
+    aM06 = 0.1016
+    aM07 = 6.756e-2
+    aM08 = -1.181e-2
+    aM09 = -4.753e-3
+    aM10 = 6.233e-4
+    aM11 = 5.805e-5
+
+    return np.where((D > 0.01) & (D < 1000.),
+        10**(aM00
+        + aM01 * np.log10(D)**1
+        + aM02 * np.log10(D)**2
+        + aM03 * np.log10(D)**3
+        + aM04 * np.log10(D)**4
+        + aM05 * np.log10(D)**5
+        + aM06 * np.log10(D)**6
+        + aM07 * np.log10(D)**7
+        + aM08 * np.log10(D)**8
+        + aM09 * np.log10(D)**9
+        + aM10 * np.log10(D)**10
+        + aM11 * np.log10(D)**11),
+        np.full_like(D, np.nan))

ctem/__init__.py

@@ -0,0 +1 @@
+from ctem.driver import *

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()