Python functions to calculate the integral of a Sérsic profile

sersicint.py

# Functions to calculate the integral of a Sérsic profile

# Sérsic profiles may be defined in various ways.
# Make sure you are using the correct quantities (e.g. not I0, mu_e, etc.)
# Ie : Intensity at the effective radius
# re : effective radius
# n : sersic index
# r : radius in same units as re

# Equations are taken from Graham and Driver (2005):
# https://ui.adsabs.harvard.edu/abs/2005PASA...22..118G/abstract

from numpy import pi, exp
from scipy.special import gamma, gammaincinv, gammainc

def b(n):
    # Normalisation constant
    return gammaincinv(2*n, 0.5)

def create_sersic_function(Ie, re, n):
    # Not required for integrals - provided for reference
    # This returns a "closure" function, which is fast to call repeatedly with different radii
    neg_bn = -b(n)
    reciprocal_n = 1.0/n
    f = neg_bn/re**reciprocal_n
    def sersic_wrapper(r):
        return Ie * exp(f * r ** reciprocal_n - neg_bn)
    return sersic_wrapper

def sersic_lum(Ie, re, n):
    # total luminosity (integrated to infinity)
    bn = b(n)
    g2n = gamma(2*n)
    return Ie * re**2 * 2*pi*n * exp(bn)/(bn**(2*n)) * g2n

def sersic_enc_lum(r, Ie, re, n):
    # luminosity enclosed within a radius r
    x = b(n) * (r/re)**(1.0/n)
    return sersic_lum(Ie, re, n) * gammainc(2*n, x)