tdose_utilities.py 文件源码

def get_2DsersicIeff(value,reff,sersicindex,axisratio,boxiness=0.0,returnFtot=False):
    """
    Get the surface brightness value at the effective radius of a 2D sersic profile (given GALFIT Sersic parameters).
    Ieff is calculated using ewuations (4) and (5) in Peng et al. (2010), AJ 139:2097.
    This Ieff is what is referred to as 'amplitude' in astropy.modeling.models.Sersic2D
    used in tdose_utilities.gen_2Dsersic()

    --- INPUT ---
    value         If returnFtot=False "value" corresponds to Ftot of the profile (total flux for profile integrated
                  til r=infty) and Ieff will be returned.
                  If instead returnFtot=True "value" should provide Ieff so Ftot can be returned
    reff          Effective radius
    sersicindex   Sersic index of profile
    axisratio     Ratio between the minor and major axis (0<axisratio<1)
    boxiness      The boxiness of the profile
    returnFtot    If Ftot is not known, but Ieff is, set returnFtot=True to return Ftot instead (providing Ieff to "value")

    --- EXAMPLE OF USE ---
    Ieff        = 1.0
    reff        = 25.0
    sersicindex = 4.0
    axisratio   = 1.0
    Ftot_calc   = tu.get_2DsersicIeff(Ieff,reff,sersicindex,axisratio,returnFtot=True)
    Ieff_calc   = tu.get_2DsersicIeff(Ftot_calc,reff,sersicindex,axisratio)

    size = 1000
    x,y = np.meshgrid(np.arange(size), np.arange(size))
    mod = Sersic2D(amplitude = Ieff, r_eff = reff, n=sersicindex, x_0=size/2.0, y_0=size/2.0, ellip=1-axisratio, theta=-1)
    img = mod(x, y)
    hducube  = pyfits.PrimaryHDU(img)
    hdus = [hducube]
    hdulist = pyfits.HDUList(hdus)
    hdulist.writeto('/Volumes/DATABCKUP2/TDOSEextractions/models_cutouts/model_sersic_spherical.fits',clobber=True)

    """
    gam2n  = scipy.special.gamma(2.0*sersicindex)
    kappa  = scipy.special.gammaincinv(2.0*sersicindex,0.5)
    Rfct   = np.pi * (boxiness + 2.) / (4. * scipy.special.beta(1./(boxiness+2.),1.+1./(boxiness+2.)) )
    factor = 2.0 * np.pi * reff**2.0 * np.exp(kappa) * sersicindex * kappa**(-2*sersicindex) * gam2n * axisratio / Rfct

    if returnFtot:
        Ftot = value * factor
        return Ftot
    else:
        Ieff  = value / factor
        return Ieff
# = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = = =