sersic.py 文件源码

def bg_lum(pp, qq, reff=1.0, lum=1.0):
    # Tested 2011-08-31
    """Exact deprojection of Sersic profile using Meijer G function as
    described by Baes + Gentile arxiv:1009.4713.

    pp and qq are the numerator and denominator of the Sersic index
      (both integers) so that n=pp/qq, 
    reff is the projected half-light radius
    lum is the total luminosity.  

    This returns a function that takes a radius and returns a
    luminosity density.

    >>> lum = luminosity(5,3)
    >>> lum(1.1)
    >>> lum([1.1, 2.2, 3.3])

    """
    if not (pp==int(pp) and qq==int(qq)): raise RuntimeError
    pp, qq = int(pp), int(qq)
    i0, bb = bg_constants(pp, qq)

    # Solution gets slow for larger p,q, so make sure that fraction is reduced
    the_gcf = euclid_gcf(pp,qq)    
    if the_gcf != 1: return bg_lum(pp/the_gcf, qq/the_gcf)        

    # a and b vectors are specified: [[num1, num2, num3],[denom1,denom2,denom3]]    
    avect = [[], [xx/(1.0*qq) for xx in range(1,qq)]]
    bvect = [[xx/(2.0*pp) for xx in range(1,2*pp)] +
             [xx/(2.0*qq) for xx in range(1,2*qq,2)], []]
    factor = 2*i0*np.sqrt(pp*qq)/(reff*(2*np.pi)**pp)

    def luminosity(rr):
        if np.iterable(rr): return np.array([luminosity(r) for r in rr])
        ss = rr/reff
        zz = (bb/(2*pp))**(2*pp) * ss**(2*qq)
        return lum*((factor/ss)*mpmath.meijerg(avect, bvect, zz))

    return luminosity