python imf_cdf.py

Posted 2021-05-09

from scipy.special import hyp2f1
import numpy as np

def get_pdf(mass, alpha, beta, mstar=1):
    return mass**(-alpha) * (1 - np.exp(-(mass / mstar)**(alpha + beta)))


def get_cdf(mass, alpha, beta, mstar=1):
    f1 = lambda x, a, b, c: x**(-b) * (1 - np.exp(-(x / c)**(a + b)))

    lnx, dlnx = np.linspace(-5, 5, 1000, retstep=True)
    x, dx = np.exp(lnx), np.exp(lnx) * dlnx
    y = (f1(x, alpha, beta, mstar)).cumsum() * dx
    x = np.exp(lnx + 0.5 * dlnx)

    return np.interp(mass, x, y)


def get_pdf(mass, alpha, beta, mstar=1, flex=1):
    x, a, b, c, d = mass, alpha, beta, flex, mstar
    return x**a * (1 + (x / d)**c)**((b - a) / c)


def get_cdf(mass, alpha, beta, mstar=1, flex=1):
    x, a, b, c, d = mass, alpha, beta, flex, mstar
    return (x**(1 + a) * hyp2f1((1 + a) / c, (a - b) / c, (1 + a + c) / c, -(x / d)**c)) / (1 + a)