Add an example for how to generate synthetic crater distributions using craterstats?

This is maybe more of a technical question on poisson processes than a craterstats question, but could improve the software documentation.

I am trying to generate synthetic crater size distributions, generated for varying timespans (varying absolute start times and durations). To begin, I am trying to reproduce the workflow and Figure 1a in Michael et al., 2016 "Planetary surface dating from crater size-frequency distribution measurements: Poisson timing analysis".

I'm having trouble with the step where "the instantaneous cratering rate λt for craters larger than dmin is found from the CF and PF." I think this question could also be a helpful example for using the craterstats package without the cli, if you decide to take the package in the direction of supporting a pythonic api.

T = 1 # 1 billion years

# define the diameters of interest and a discretized reference between them
crater_range_interest = [50/1000, 10] # km, as set in Figure 1
x0 = np.linspace(-10, 4, 5000)
production_valid = np.where(np.logical_and(10**x0 >= production_function.range[0], 10**x0 <= production_function.range[1]))
x = x0[production_valid]

# set up random number gen and lists to accumulate craters generated
rng = np.random.default_rng()
list_d = []   # list of diameters
list_dt = []  # list of interarrival times

# set dmin and area to match Figure 1
dmin = crater_range_interest[0]
Area = 100

# set up time/counter
t = -T
it = 0

# run the simulation while time is less than 0 (or alternatively up to the number of craters in Figure 1)
while t < 0:
# while it < 2590: # i.e., the number of craters in Figure 1 data

    U = rng.uniform(0, 1)
    a0 = chronology_function.a0(np.abs(t))  # or 1 always?

    # calculate the cumulative form of production function
    prod0 = (U * production_function.C(dmin, a0))  # cumulative produced /km2/Ga for reference d
    cumulative_produced = production_function.C(10**x, a0)  # cumulative produced for each size, /km2/Ga

    # find the diameter for this event
    didx = np.argmax(prod0 >= cumulative_produced)
    d = 10**(x[didx])

    # compute lambdat from production function and chronology function
    d_production_rate = Area * (
        production_function.C(10**x[didx-1], a0) -
        production_function.C(10**x[didx], a0)) # production rate Eqn 2(ish)
    N1_production_rate = chronology_function.phi(np.abs(t)) # production rate of 1km crater at time t

    # what is the correct step here?!
    lambdat = d_production_rate / N1_production_rate  # relate the N1 density to the d production?

    # convert to deltat
    U = rng.uniform(0, 1)
    deltat = -(np.log(U) / (lambdat))  # exponential distribution of arrival times with mean lambdat

    # append the time and d to a list
    list_d.append(d)
    list_dt.append(deltat)

    t += (deltat)
    it += 1

print(f"Total time was {T+t} Ga, with {len(list_d)} craters.")

From the above code, I can generate a plot that somewhat matches Figure 1a (below), but the recovered distribution age is only ~1 Ga when the crater number is fixed to num_craters=2590, as in Figure 1a. Moreover, the recovered timescales vary substantially when using simulation durations other than 1 Ga (500 Ma, gives surface estimates of ~200 Ma).

Total time was 1.066659270408508 Ga, with 1091 craters.

So, I think I'm doing something wrong here, and it seems to have to do with determining the interarrival time via lambda_t, maybe missing a time-relevant conversion, when going from the crater density via chronology function to production rate of size d crater.

Could you provide any advice on how to simulate the crater distributions with craterstats? I'm happy to add a complete example to documentation when I sort this out.

ggmichael / craterstats

Add an example for how to generate synthetic crater distributions using craterstats? #18