Luminosity of Sérsic function in annulus

sersic_lum.py

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

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

# Total luminosity of a 2D Sérsic profile
def sersic_total_lum(Ie, re, n):
    bn = b(n)
    g2n = gamma(2*n)
    return Ie * re**2 * 2*pi*n * exp(bn)/(bn**(2*n)) * g2n

# Luminosity of a 2D Sérsic enclosed within a radius, r, 
def sersic_enclosed_lum(r, Ie, re, n):
    x = b(n) * (r/re)**(1.0/n)
    return sersic_total_lum(Ie, re, n) * gammainc(2*n, x)

# Luminosity of a 2D Sérsic enclosed between two radii, r1 and r2 
def sersic_annular_lum(r1, r2, Ie, re, n):
    lum_r2 = sersic_enclosed_lum(r2, Ie, re, n)
    lum_r1 = sersic_enclosed_lum(r1, Ie, re, n)
    return abs(lum_r2 - lum_r1)

# Determine the Ie required for a Sérsic profile with a given total luminosity
def sersic_Ie_from_total_lum(lum, re, n):
    return lum / sersic_total_lum(1, re, n)

# Determine the Ie required for a Sérsic profile with a given annular luminosity
def sersic_Ie_from_annular_lum(lum, re, n):
    return lum / sersic_annular_lum(1, re, n)