Neukum production function

Question

I am trying to plot a function which looks like that: . Where, D range from 0.01 to 300km, N is number of craters with dia > D per sq.km per giga year, coefficients values are in code.

I get empty plot and obviously I am doing something very wrong which I am not able to understand. I am sharing my code.

'''

a_0 = -3.0876
a_1 = -3.557528
a_2 = 0.781027
a_3 = 1.021521
a_4 = -0.156012
a_5 = -0.444058
a_6 = 0.019977
a_7 = 0.086850
a_8 = -0.005874
a_9 = -0.006809
a_10 = 8.25*10**-4
a_11 = 5.54*10**-5
a_n = a_0 + a_1 + a_2 + a_3 + a_4 + a_5 + a_6 + a_7 + a_8 + a_9 + a_10 + a_11

for d in range(1, 200000):
    n = a_0 + np.multiply(a_1, np.log(d))
    + np.multiply(a_2, np.log(d)**2)
    + np.multiply(a_3, np.log(d)**3)
    + np.multiply(a_4, np.log(d)**4)
    + np.multiply(a_5, np.log(d)**5)
    + np.multiply(a_6, np.log(d)**6)
    + np.multiply(a_7, np.log(d)**7)
    + np.multiply(a_8, np.log(d)**8)
    + np.multiply(a_9, np.log(d)**9)
    + np.multiply(a_10, np.log(d)**10)
    + np.multiply(a_11, np.log(d)**11)
print(10**n)
plt.plot(10**n, color='green')
plt.show()

''' The graphical curve should look something like that

You calculate the value n 200000 times, but you only save the last one. Maybe also save all the n values in an array? — JohanC, Nov 18 '19 at 19:26
My advice is to reduce the number of values to something much smaller, say 20, and print out all of the values in the array before trying to plot it. That might shed some light. — Robert Dodier, Nov 18 '19 at 19:58

score 3 · Accepted Answer · answered Nov 18 '19 at 20:48

I can not understand your code, but if you want to plot the function you wrote you can try this:

import numpy as np
import matplotlib.pyplot as plt
a_0 = -3.0876
a_1 = -3.557528
a_2 = 0.781027
a_3 = 1.021521
a_4 = -0.156012
a_5 = -0.444058
a_6 = 0.019977
a_7 = 0.086850
a_8 = -0.005874
a_9 = -0.006809
a_10 = 8.25*10**(-4)
a_11 = 5.54*10**(-5)

N_POINTS = 100
a_coeff = np.array([a_0, a_1, a_2, a_3, a_4, a_5, a_6, a_7, a_8, a_9, a_10, a_11])
distance = np.logspace(np.log10(0.01), np.log10(200), N_POINTS)

exponents = np.arange(12)
distance_matrix = distance[:, np.newaxis]*np.ones([N_POINTS, 12])
N = 10**np.sum(a_coeff * (np.log10(distance_matrix)**exponents), axis=1)

fig, ax = plt.subplots(1, figsize=(3, 8))
ax.scatter(distance, N)
ax.set(xscale='log', yscale='log', ylim=(1e-7, 1e4), xlim=(0.001, 300))

Neukum production function

1 Answers1