Python Sed.flambdaTofnu方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import flambdaTofnu [as 别名]
    def test_eq(self):
        """
        Test that __eq__ in Sed works correctly
        """
        sed_dir = os.path.join(getPackageDir('sims_photUtils'), 'tests',
                               'cartoonSedTestData', 'starSed', 'kurucz')
        list_of_seds = os.listdir(sed_dir)
        sedname1 = os.path.join(sed_dir, list_of_seds[0])
        sedname2 = os.path.join(sed_dir, list_of_seds[1])
        ss1 = Sed()
        ss1.readSED_flambda(sedname1)
        ss2 = Sed()
        ss2.readSED_flambda(sedname2)
        ss3 = Sed()
        ss3.readSED_flambda(sedname1)

        self.assertFalse(ss1 == ss2)
        self.assertTrue(ss1 != ss2)
        self.assertTrue(ss1 == ss3)
        self.assertFalse(ss1 != ss3)

        ss3.flambdaTofnu()

        self.assertFalse(ss1 == ss3)
        self.assertTrue(ss1 != ss3)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import flambdaTofnu [as 别名]
    def calcMagNorm(self, objectMags, sedObj, bandpassDict, mag_error = None,
                    redshift = None, filtRange = None):

        """
        This will find the magNorm value that gives the closest match to the magnitudes of the object
        using the matched SED. Uses scipy.optimize.leastsq to find the values of fluxNorm that minimizes
        the function: ((flux_obs - (fluxNorm*flux_model))/flux_error)**2.

        @param [in] objectMags are the magnitude values for the object with extinction matching that of
        the SED object. In the normal case using the selectSED routines above it will be dereddened mags.

        @param [in] sedObj is an Sed class instance that is set with the wavelength and flux of the
        matched SED

        @param [in] bandpassDict is a BandpassDict class instance with the Bandpasses set to those
        for the magnitudes given for the catalog object

        @param [in] mag_error are provided error values for magnitudes in objectMags. If none provided
        then this defaults to 1.0. This should be an array of the same length as objectMags.

        @param [in] redshift is the redshift of the object if the magnitude is observed

        @param [in] filtRange is a selected range of filters specified by their indices in the bandpassList
        to match up against. Used when missing data in some magnitude bands.

        @param [out] bestMagNorm is the magnitude normalization for the given magnitudes and SED
        """

        import scipy.optimize as opt

        sedTest = Sed()
        sedTest.setSED(sedObj.wavelen, flambda = sedObj.flambda)
        if redshift is not None:
            sedTest.redshiftSED(redshift)
        imSimBand = Bandpass()
        imSimBand.imsimBandpass()
        zp = -2.5*np.log10(3631)  #Note using default AB zeropoint
        flux_obs = np.power(10,(objectMags + zp)/(-2.5))
        sedTest.resampleSED(wavelen_match=bandpassDict.wavelenMatch)
        sedTest.flambdaTofnu()
        flux_model = sedTest.manyFluxCalc(bandpassDict.phiArray, bandpassDict.wavelenStep)
        if filtRange is not None:
            flux_obs = flux_obs[filtRange]
            flux_model = flux_model[filtRange]
        if mag_error is None:
            flux_error = np.ones(len(flux_obs))
        else:
            flux_error = np.abs(flux_obs*(np.log(10)/(-2.5))*mag_error)
        bestFluxNorm = opt.leastsq(lambda x: ((flux_obs - (x*flux_model))/flux_error), 1.0)[0][0]
        sedTest.multiplyFluxNorm(bestFluxNorm)
        bestMagNorm = sedTest.calcMag(imSimBand)
        return bestMagNorm

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import flambdaTofnu [as 别名]
def get_TotalMags(result, bandpasses=('u','g','r','i','z','y')):
    datadir = os.environ.get("SIMS_SED_LIBRARY_DIR")
    tpath = os.getenv('LSST_THROUGHPUTS_DEFAULT')
    bands = {"u":None, "g":None, "r":None, "i":None, "z":None, "y":None}
    for k in bands:
        bands[k] = Bandpass()
        bands[k].readThroughput(os.path.join(tpath, "total_%s.dat"%k))
    # Set up phi, the wavelength-normalized system response for each filter,
    # for each bandpass for manyMagCalc method.
    bplist = []
    for f in ['u','g','r','i','z','y']:
        bands[f].sbTophi()
        bplist.append(bands[f])
    ids = result['galid']
    diskfile = result['sedFilenameDisk']
    bulgefile = result['sedFilenameBulge']
    agnfile = result['sedFilenameAgn']

    diskmn = result['magNormDisk']
    bulgemn = result['magNormBulge']
    agnmn = result['magNormAgn']

    bulgeAv = result['internalAvBulge']
    diskAv = result['internalAvDisk']

    redshift = result['redshift']

    imsimband = Bandpass()
    imsimband.imsimBandpass()
    sedDict = {}
    retMags = dict([(k, []) for k in bands])
    a_int = None
    b_int = None
    tmpwavelen = None
    for id, df, dm, dav, bf, bm, bav, af, am, z in zip(ids, diskfile, diskmn, diskAv, 
            bulgefile, bulgemn, bulgeAv, agnfile, agnmn, redshift):
        tmpflux = None
        for comp in ((df, dm, dav, 'galaxySED', False), (bf, bm, bav, 'galaxySED', False), (af, am, None, 'agnSED', True)):
        #Zero out the AGN contribution
        #for comp in ((df, dm, dav, 'galaxySED', False), (bf, bm, bav, 'galaxySED', False), (af, 99.99, None, 'agnSED', True)):
            if not comp[0] == u'None':
                if sedDict.has_key(comp[0]):
                    sed = copy.deepcopy(sedDict[comp[0]])
                else:
                    sed = Sed()
                    print os.path.join(datadir,comp[3],comp[0])
                    sed.readSED_flambda(os.path.join(datadir,comp[3],comp[0]))
		    if comp[4]:
		        sed.resampleSED(wavelen_match=tmpwavelen)
                    sedDict[comp[0]] = sed
                if a_int is None:
                    phiarray, dlambda = sed.setupPhiArray(bplist)
                    a_int, b_int = sed.setupCCMab()
		    #Careful, this assumes that a disk or bulge sed is read
		    #before any agn sed
		    tmpwavelen = sed.wavelen
                fNorm = sed.calcFluxNorm(comp[1], imsimband)
                sed.multiplyFluxNorm(fNorm)
                #I guess this assumes rv=3.1??
                if comp[2]:
                    sed.addCCMDust(a_int, b_int, A_v=comp[2])
		wavelenArr=sed.wavelen
		if tmpflux is None:
		    tmpflux = sed.flambda
		else:
	            tmpflux += sed.flambda
	newgal = Sed(wavelen=wavelenArr, flambda=tmpflux)
        #a_mw, b_mw = sed.setupCCMab()
        #sed.addCCMDust(a_mw, b_mw, A_v=mwav)
        newgal.redshiftSED(z, dimming=True)
	newgal.resampleSED(wavelen_match=bplist[0].wavelen)
	newgal.flambdaTofnu()
        mags = newgal.manyMagCalc(phiarray, dlambda)
        for i,k in enumerate(['u','g','r','i','z','y']):
            retMags[k].append(mags[i])
    return retMags

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import flambdaTofnu [as 别名]
# Calculate dust extinction a_x and b_x vectors.
a_mw, b_mw = stars[starlist[0]].setupCCMab()
# Set up dictionary + arrays to hold calculated magnitude information.
mags1 = {}
for f in filterlist:
    mags1[f] = numpy.zeros(num_star, dtype="float")
# For each star (in num_star's), apply apply MW dust, fluxnorm & calculate mags.
for i in range(num_star):
    starname = starlist[star_name[i]]
    tmpstar = Sed(wavelen=stars[starname].wavelen, flambda=stars[starname].flambda)
    tmpstar.addCCMDust(a_mw, b_mw, ebv=ebv_mw[i])
    tmpstar.multiplyFluxNorm(fluxnorm[i])
    # Note that the stars have already been matched to the bandpass wavelength grid.
    # Just want to be sure that fnu is already calculated.
    tmpstar.flambdaTofnu()
    for f in filterlist:
        mags1[f][i] = tmpstar.calcMag(lsstbp[f])
dt, t = dtime(t)
print "Calculating dust/fluxnorm/%d magnitudes with some smart usage for %d stars took %f s" % (
    len(filterlist),
    num_star,
    dt,
)


# Test Sed.manyMagCalc :

# First: (re) calculate internal a/b on wavelength range required for dust extinction.
a_mw, b_mw = stars[starlist[0]].setupCCMab()
# Also: set up phi for each bandpass - ahead of time. And set up a list of bandpasses, then create phiarray