Search code examples
pythonfunctionnumpymatplotlibastronomy

Defining and plotting a Schechter function: plot problems

I'm currently defining a function in python as:

def schechter_fit(logM, phi=5.96E-11, log_M0=11.03, alpha=-1.35, e=2.718281828):
    schechter = phi*(10**((alpha+1)*(logM-log_M0)))*(e**(pow(-10,logM-log_M0)))
    return schechter

schechter_range = numpy.linspace(10.0, 11.9, 10000)

And then plotting said function as:

import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid.axislines import SubplotZero

schechter_range = numpy.linspace(10, 12, 10000)

fig = plt.figure(1)
ax = SubplotZero(fig, 111)
fig.add_subplot(ax)

ax.plot(schechter_range, schechter_fit(schechter_range), 'k')

This is the graphical output I am receiving is just a blank plot with no curve plotted. There must be a problem with how I have defined the function, but I can't see the problem. The plot should look something like this:

enter image description here

I'm new to python functions so perhaps my equation isn't quite right. This is what I am looking to plot and the parameters I am starting with:

enter image description here

Solution

  • The function you describe returns a complex result over most of your input range. Here I added +0j to the input to allow for an imaginary result; if you don't do this you just get a bunch of nans (which mpl doesn't plot). Here are the plots:

    enter image description here

    import numpy
import matplotlib.pyplot as plt
from mpl_toolkits.axes_grid.axislines import SubplotZero

schechter_range = numpy.linspace(10, 12, 10000)

fig = plt.figure(1)
ax = SubplotZero(fig, 111)
fig.add_subplot(ax)

def schechter_fit(logM, phi=5.96E-11, log_M0=11.03, alpha=-1.35, e=2.718281828):
    schechter = phi*(10**((alpha+1)*(logM-log_M0)))*(e**(pow(-10,logM-log_M0)))
    return schechter

y = schechter_fit(schechter_range+0j)  # Note the +0j here to allow an imaginary result
ax.plot(schechter_range, y.real, 'b', label="Re Part")
ax.plot(schechter_range, y.imag, 'r', label="Im Part")
ax.legend()
plt.show()

    Now that you can see why the data is not plotting, and that complex numbers are being generated, and you know physically that you don't want that, it would be reasonable to figure out where these are coming from. Hopefully, it's obvious that these are originate from pow(-10,logM-log_M0), and from there it's clear that this is assuming the wrong operator precedence: the equation isn't pow(-10,logM-log_M0), but -pow(10,logM-log_M0). Making this corrections gives (after a log is taken, because I can see the log in the plot in the question): enter image description here

    I also extended the lower bound from 10 to 8, so the region of constant slope is clear and it better matches the graph shown in the question. This is still off by a factor on the y-axis, but I'm guessing that's a factor of (SFR/M*) that's not being applied correctly (it's difficult to know without seeing the context and the full y-axis).