Python Sed.redshiftSED方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
 def testRedshiftName(self):
     testsed = Sed(self.testsed.wavelen, self.testsed.flambda, name=self.testsed.name)
     redshift = .2
     testsed.redshiftSED(redshift=redshift)
     newname = testsed.name + '_Z' + '%.2f' % (redshift)
     testsed.name = newname
     self.assertEqual(testsed.name, newname)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [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 redshiftSED [as 别名]
def calcADUwrapper(sedName=None, magNorm=None, redshift=None, internalAv=None, internalRv=None,
                   galacticAv=None, galacticRv=None, bandpass=None):
    """
    Read in an SED and calculat the number of ADU produced by that SED in a specified bandpass

    Parameters
    ----------
    sedName is a string specifying the file name of the SED

    magNorm is the normalizing magnitude of the SED in the imsimBandpass

    redshift is the redshift of the SED

    internalAv is the Av due to internal dust of the source (if a galaxy)

    internalRv is the Rv due to internal dust of the source (if a galaxy)

    galacticAv is the Av due to Milky Way dust between observer and source

    galacticRv is the Rv due to Milky Way dust between observer and source

    bandpass is an intantiation of Bandpass representing the band in which the ADUs are measured

    Returns
    -------
    A float representing the number of ADUs measured in the bandpass
    """

    imsimband = Bandpass()
    imsimband.imsimBandpass()
    sed = Sed()
    sed.readSED_flambda(sedName)
    fNorm = sed.calcFluxNorm(magNorm, imsimband)
    sed.multiplyFluxNorm(fNorm)
    if internalAv is not None and internalRv is not None:
        if internalAv != 0.0 and internalRv != 0.0:
            a_int, b_int = sed.setupCCM_ab()
            sed.addDust(a_int, b_int, A_v=internalAv, R_v=internalRv)

    if redshift is not None and redshift != 0.0:
        sed.redshiftSED(redshift, dimming=True)

    a_int, b_int = sed.setupCCM_ab()
    sed.addDust(a_int, b_int, A_v=galacticAv, R_v=galacticRv)

    adu = sed.calcADU(bandpass, photParams=PhotometricParameters())

    return adu

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
def calcADUwrapper(sedName=None, magNorm=None, redshift=None, internalAv=None, internalRv=None,
                   galacticAv=None, galacticRv=None, bandpass=None):

    imsimband = Bandpass()
    imsimband.imsimBandpass()
    sed = Sed()
    sed.readSED_flambda(sedName)
    fNorm = sed.calcFluxNorm(magNorm, imsimband)
    sed.multiplyFluxNorm(fNorm)
    if internalAv is not None and internalRv is not None:
        if internalAv != 0.0 and internalRv != 0.0:
            a_int, b_int = sed.setupCCMab()
            sed.addCCMDust(a_int, b_int, A_v=internalAv, R_v=internalRv)
    
    if redshift is not None and redshift != 0.0:
        sed.redshiftSED(redshift, dimming=True)
    
    a_int, b_int = sed.setupCCMab()
    sed.addCCMDust(a_int, b_int, A_v=galacticAv, R_v=galacticRv)
    
    adu = sed.calcADU(bandpass, photParams=PhotometricParameters())
    
    return adu

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
def read_quasar():
    # read quasar spectra and redshift
    homedir = os.getenv("HOME")
    quasardir = os.path.join(homedir, "seds/quasar")
    # read zero redshift quasar
    base = Sed()
    base.readSED_flambda(os.path.join(quasardir, "quasar.dat"))
    # redshift 
    #redshifts = [0, 0.1, 0.2, 0.3, 0.5, 0.8, 1.0, 1.3, 1.6, 1.9, 2.2, 2.5]
    #redshifts = numpy.array(redshifts)
    redshifts= numpy.arange(0, 2.8, 0.1)
    quasars = {}
    for z in redshifts:
        wavelen, flambda = base.redshiftSED(z, wavelen=base.wavelen, flambda=base.flambda)
        quasars[z] = Sed(wavelen=wavelen, flambda=flambda)
    print "# Generated %d quasars at redshifts between %f and %f" %(len(redshifts), redshifts.min(), redshifts.max())
    # resample onto the standard bandpass for Bandpass obj's and calculate fnu to speed later calculations
    for z in redshifts:
        quasars[z].synchronizeSED(wavelen_min=WMIN, wavelen_max=WMAX, wavelen_step=WSTEP)
    return quasars, redshifts

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]

#.........这里部分代码省略.........
        @param [out] magNormMatches are the magnitude normalizations for the given magnitudes and
        matched SED.

        @param [out] matchErrors contains the Mean Squared Error between the colors of each object and 
        the colors of the matched SED.
        """

        #Set up photometry to calculate model Mags
        if bandpassDict is None:
            galPhot = BandpassDict.loadTotalBandpassesFromFiles(['u','g','r','i','z'],
                                            bandpassDir = os.path.join(lsst.utils.getPackageDir('throughputs'),'sdss'),
                                            bandpassRoot = 'sdss_')
        else:
            galPhot = bandpassDict

        #Calculate ebv from ra, dec coordinates if needed
        if reddening == True:
            #Check that catRA and catDec are included
            if catRA is None or catDec is None:
                raise RuntimeError("Reddening is True, but catRA and catDec are not included.")
            calcEBV = ebv()
            raDec = np.array((catRA,catDec))
            #If only matching one object need to reshape for calculateEbv
            if len(raDec.shape) == 1:
                raDec = raDec.reshape((2,1))
            ebvVals = calcEBV.calculateEbv(equatorialCoordinates = raDec)
            objMags = self.deReddenMags(ebvVals, catMags, extCoeffs)
        else:
            objMags = catMags

        minRedshift = np.round(np.min(catRedshifts), dzAcc)
        maxRedshift = np.round(np.max(catRedshifts), dzAcc)
        dz = np.power(10., (-1*dzAcc))

        redshiftRange = np.round(np.arange(minRedshift - dz, maxRedshift + (2*dz), dz), dzAcc)
        numRedshifted = 0
        sedMatches = [None] * len(catRedshifts)
        magNormMatches = [None] * len(catRedshifts)
        matchErrors = [None] * len(catRedshifts)
        redshiftIndex = np.argsort(catRedshifts)

        numOn = 0
        notMatched = 0
        lastRedshift = -100
        print('Starting Matching. Arranged by redshift value.')
        for redshift in redshiftRange:

            if numRedshifted % 10 == 0:
                print('%i out of %i redshifts gone through' % (numRedshifted, len(redshiftRange)))
            numRedshifted += 1

            colorSet = []
            for galSpec in sedList:
                sedColors = []
                fileSED = Sed()
                fileSED.setSED(wavelen = galSpec.wavelen, flambda = galSpec.flambda)
                fileSED.redshiftSED(redshift)
                sedColors = self.calcBasicColors([fileSED], galPhot, makeCopy = True)
                colorSet.append(sedColors)
            colorSet = np.transpose(colorSet)
            for currentIndex in redshiftIndex[numOn:]:
                matchMags = objMags[currentIndex]
                if lastRedshift < np.round(catRedshifts[currentIndex],dzAcc) <= redshift:
                    colorRange = np.arange(0, len(galPhot)-1)
                    matchColors = []
                    for colorNum in colorRange:
                        matchColors.append(matchMags[colorNum] - matchMags[colorNum+1])
                    #This is done to handle objects with incomplete magnitude data
                    filtNums = np.arange(0, len(galPhot))
                    if np.isnan(np.amin(matchColors))==True:
                        colorRange = np.where(np.isnan(matchColors)==False)[0]
                        filtNums = np.unique([colorRange, colorRange+1]) #Pick right filters in calcMagNorm
                    if len(colorRange) == 0:
                        print('Could not match object #%i. No magnitudes for two adjacent bandpasses.' \
                              % (currentIndex))
                        notMatched += 1
                        #Don't need to assign 'None' here in result array, b/c 'None' is default value
                    else:
                        distanceArray = [np.zeros(len(sedList))]
                        for colorNum in colorRange:
                            distanceArray += np.power((colorSet[colorNum] - matchColors[colorNum]),2)
                        matchedSEDNum = np.nanargmin(distanceArray)
                        sedMatches[currentIndex] = sedList[matchedSEDNum].name
                        magNormVal = self.calcMagNorm(np.array(matchMags), sedList[matchedSEDNum], 
                                                      galPhot, mag_error = mag_error,
                                                      redshift = catRedshifts[currentIndex],
                                                      filtRange = filtNums)
                        magNormMatches[currentIndex] = magNormVal
                        matchErrors[currentIndex] = (distanceArray[0,matchedSEDNum]/len(colorRange))
                    numOn += 1
                else:
                    break
            lastRedshift = redshift

        print('Done Matching. Matched %i of %i catalog objects to SEDs' % (len(catMags)-notMatched,
                                                                           len(catMags)))
        if notMatched > 0:
            print('%i objects did not get matched.' % (notMatched))

        return sedMatches, magNormMatches, matchErrors

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
counts = star.calcADU(rband, expTime=30)
print "This would correspond to roughly %f counts in the LSST focal plane, in a 30s exposure." %(counts)

# For fun, let's see what else can happen.

ebv = 0.5
print ""
print "Let's try adding %.2f E(B-V) dust extinction to this star." %(ebv)
a, b = star.setupCCMab()
# You can use addCCMDust on the 'star' object itself, but I'm illustrating here how you could also
#  do otherwise - preserve the original 'star' object as is, and create a new Sed object that does
#  include the effects of dust ('dustystar').  Star's data will be unchanged by the dust. 
dustywavelen, dustyflambda = star.addCCMDust(a, b, ebv=ebv, wavelen=star.wavelen, flambda=star.flambda)
dustystar = Sed(wavelen=dustywavelen, flambda=dustyflambda)
magdust = dustystar.calcMag(rband)
print "With this dust, the magnitude of the star in this bandpass is now %.4f." %(magdust)

redshift = 0.2
print "What if this star was at a redshift of %f?" %(redshift)
# Here (unlike above with the dust), I'm applying the redshift to the 'star' object itself.
#  Star's data will be changed by the redshifting. 
star.redshiftSED(redshift, dimming=True)
magredshift = star.calcMag(rband)
print ""
print "Redshifted to %.2f, and adding cosmological dimming, the magnitude is now %.4f" \
      %(redshift, magredshift)
print " (this was a pretty hot star, so redshifting brings more flux into this particular bandpass.)"
print ""

# There is more functionality in the class. Please see the code. 

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
# create the a,b arrays for all the gals (because we resampled the gals onto the
#  same wavelength range we can just calculate a/b once, and this is slow)

a_gal, b_gal = gals[galaxykeys[0]].setupCCMab()

# pretend we want to read mags into an array .. you could just as easily put it into a
# dictionary or list, with small variations in the code
mags = n.empty(len(galaxykeys), dtype='float')

for i in range(len(galaxykeys)):
    # make a copy of the original SED if you want to 'reuse' the SED for multiple magnitude
    # calculations with various fluxnorms, redshifts and dusts
    tmpgal = Sed(wavelen=gals[galaxykeys[i]].wavelen, flambda=gals[galaxykeys[i]].flambda)
    # add the dust internal to the distant galaxy
    tmpgal.addCCMDust(a_gal, b_gal, ebv=ebv_gal[i])
    # redshift the galaxy
    tmpgal.redshiftSED(redshifts[i], dimming=False)
    # add the dust from our milky way - have to recalculate a/b because now wavelenghts
    # for each galaxy are *different*
    a_mw, b_mw = tmpgal.setupCCMab()
    tmpgal.addCCMDust(a_mw, b_mw, ebv=ebv_mw[i])
    tmpgal.multiplyFluxNorm(fluxnorm[i])
    mags[i] = tmpgal.calcMag(rband)


# show results
print "#sedname      fluxnorm     redshift  ebv_gal   ebv_mw  magnitude "
for i in range(len(galaxykeys)):
    print "%s %.5g  %.3f %.5f %.5f %.5f" %(galaxykeys[i], fluxnorm[i], redshifts[i], ebv_gal[i], ebv_mw[i], mags[i])
    

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [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 redshiftSED [as 别名]
    def testAddingNonesToList(self):
        """
        Test what happens if you add SEDs to an SedList that have None for
        one or more of the physical parameters (i.e. galacticAv, internalAv, or redshift)
        """
        imsimBand = Bandpass()
        imsimBand.imsimBandpass()
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                 redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                 wavelenMatch=wavelen_match)


        sedNameList_1 = self.getListOfSedNames(nSed)
        magNormList_1 = list(numpy.random.random_sample(nSed)*5.0 + 15.0)
        internalAvList_1 = list(numpy.random.random_sample(nSed)*0.3 + 0.1)
        redshiftList_1 = list(numpy.random.random_sample(nSed)*5.0)
        galacticAvList_1 = list(numpy.random.random_sample(nSed)*0.3 + 0.1)

        internalAvList_1[0] = None
        redshiftList_1[1] = None
        galacticAvList_1[2] = None

        internalAvList_1[3] = None
        redshiftList_1[3] = None

        internalAvList_1[4] = None
        galacticAvList_1[4] = None

        redshiftList_1[5] = None
        galacticAvList_1[5] = None

        internalAvList_1[6] = None
        redshiftList_1[6] = None
        galacticAvList_1[6] = None

        testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                  internalAvList=internalAvList_1,
                                  galacticAvList=galacticAvList_1,
                                  redshiftList=redshiftList_1)

        self.assertEqual(len(testList), 2*nSed)
        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

        for ix in range(len(sedNameList_0)):
            self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
            self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
            self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)


        for ix in range(len(sedNameList_1)):
            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix+nSed], 10)
            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix+nSed], 10)
            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix+nSed], 10)

        for ix, (name, norm, iav, gav, zz) in \
        enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                      galacticAvList_0, redshiftList_0)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)


        for ix, (name, norm, iav, gav, zz) in \
        enumerate(zip(sedNameList_1, magNormList_1, internalAvList_1, \
                      galacticAvList_1, redshiftList_1)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            if iav is not None:
                a_coeff, b_coeff = sedControl.setupCCMab()
                sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
    def testGalaxyPhotometricUncertainties(self):
        """
        Test in the case of a catalog of galaxies
        """
        lsstDefaults = LSSTdefaults()
        phot = PhotometryGalaxies()
        galDB = testGalaxyTileDBObj(driver=self.driver, host=self.host, database=self.dbName)
        galCat = testGalaxyCatalog(galDB, obs_metadata=self.obs_metadata)
        imsimband = Bandpass()
        imsimband.imsimBandpass()
        ct = 0
        for line in galCat.iter_catalog():
            bulgeSedName = line[50]
            diskSedName = line[51]
            agnSedName = line[52]
            magNormBulge = line[53]
            magNormDisk = line[54]
            magNormAgn = line[55]
            avBulge = line[56]
            avDisk = line[57]
            redshift = line[58]

            bulgeSed = Sed()
            bulgeSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                     defaultSpecMap[bulgeSedName]))
            fNorm=bulgeSed.calcFluxNorm(magNormBulge, imsimband)
            bulgeSed.multiplyFluxNorm(fNorm)

            diskSed = Sed()
            diskSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                    defaultSpecMap[diskSedName]))
            fNorm = diskSed.calcFluxNorm(magNormDisk, imsimband)
            diskSed.multiplyFluxNorm(fNorm)

            agnSed = Sed()
            agnSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                   defaultSpecMap[agnSedName]))
            fNorm = agnSed.calcFluxNorm(magNormAgn, imsimband)
            agnSed.multiplyFluxNorm(fNorm)

            a_int, b_int = bulgeSed.setupCCMab()
            bulgeSed.addCCMDust(a_int, b_int, A_v=avBulge)

            a_int, b_int = diskSed.setupCCMab()
            diskSed.addCCMDust(a_int, b_int, A_v=avDisk)

            bulgeSed.redshiftSED(redshift, dimming=True)
            diskSed.redshiftSED(redshift, dimming=True)
            agnSed.redshiftSED(redshift, dimming=True)

            bulgeSed.resampleSED(wavelen_match=self.totalBandpasses[0].wavelen)
            diskSed.resampleSED(wavelen_match=bulgeSed.wavelen)
            agnSed.resampleSED(wavelen_match=bulgeSed.wavelen)

            numpy.testing.assert_almost_equal(bulgeSed.wavelen, diskSed.wavelen)
            numpy.testing.assert_almost_equal(bulgeSed.wavelen, agnSed.wavelen)

            fl = bulgeSed.flambda + diskSed.flambda + agnSed.flambda

            totalSed = Sed(wavelen=bulgeSed.wavelen, flambda=fl)

            sedList = [totalSed, bulgeSed, diskSed, agnSed]

            for i, spectrum in enumerate(sedList):
                if i==0:
                    msgroot = 'failed on total'
                elif i==1:
                    msgroot = 'failed on bulge'
                elif i==2:
                    msgroot = 'failed on disk'
                elif i==3:
                    msgroot = 'failed on agn'

                for j, b in enumerate(self.bandpasses):
                    controlSigma = calcMagError_sed(spectrum, self.totalBandpasses[j],
                                             self.skySeds[j],
                                             self.hardwareBandpasses[j],
                                             FWHMeff=lsstDefaults.FWHMeff(b),
                                             photParams=PhotometricParameters())

                    testSigma = line[26+(i*6)+j]
                    msg = '%e neq %e; ' % (testSigma, controlSigma) + msgroot
                    self.assertAlmostEqual(testSigma, controlSigma, 10, msg=msg)
                    ct += 1

        self.assertGreater(ct, 0)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
    def testAddingToList(self):
        """
        Test that we can add Seds to an already instantiated SedList
        """
        imsimBand = Bandpass()
        imsimBand.imsimBandpass()
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                 redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                 wavelenMatch=wavelen_match)


        # experiment with adding different combinations of physical parameter lists
        # as None and not None
        for addIav in [True, False]:
            for addRedshift in [True, False]:
                for addGav in [True, False]:

                    testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                             redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                             wavelenMatch=wavelen_match)

                    sedNameList_1 = self.getListOfSedNames(nSed)
                    magNormList_1 = numpy.random.random_sample(nSed)*5.0 + 15.0

                    if addIav:
                        internalAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1
                    else:
                        internalAvList_1 = None

                    if addRedshift:
                        redshiftList_1 = numpy.random.random_sample(nSed)*5.0
                    else:
                        redshiftList_1 = None

                    if addGav:
                        galacticAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1
                    else:
                        galacticAvList_1 = None


                    testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                              internalAvList=internalAvList_1,
                                              galacticAvList=galacticAvList_1,
                                              redshiftList=redshiftList_1)

                    self.assertEqual(len(testList), 2*nSed)
                    numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

                    for ix in range(len(sedNameList_0)):
                        self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
                        self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
                        self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)


                    for ix in range(len(sedNameList_1)):
                        if addIav:
                            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix+nSed], 10)
                        else:
                            self.assertTrue(testList.internalAvList[ix+nSed] is None)

                        if addGav:
                            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix+nSed], 10)
                        else:
                            self.assertTrue(testList.galacticAvList[ix+nSed] is None)

                        if addRedshift:
                            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix+nSed], 10)
                        else:
                            self.assertTrue(testList.redshiftList[ix+nSed] is None)

                    for ix, (name, norm, iav, gav, zz) in \
                      enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                                  galacticAvList_0, redshiftList_0)):

                        sedControl = Sed()
                        sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

                        fnorm = sedControl.calcFluxNorm(norm, imsimBand)
                        sedControl.multiplyFluxNorm(fnorm)

                        a_coeff, b_coeff = sedControl.setupCCMab()
                        sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

                        sedControl.redshiftSED(zz, dimming=True)
                        sedControl.resampleSED(wavelen_match=wavelen_match)

                        a_coeff, b_coeff = sedControl.setupCCMab()
                        sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

                        sedTest = testList[ix]

                        numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
                        numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
    def _calculateGalSimSeds(self):
        """
        Apply any physical corrections to the objects' SEDS (redshift them, apply dust, etc.).
        Return a list of Sed objects containing the SEDS
        """

        sedList = []
        actualSEDnames = self.column_by_name('sedFilepath')
        redshift = self.column_by_name('redshift')
        internalAv = self.column_by_name('internalAv')
        internalRv = self.column_by_name('internalRv')
        galacticAv = self.column_by_name('galacticAv')
        galacticRv = self.column_by_name('galacticRv')
        magNorm = self.column_by_name('magNorm')

        #for setting magNorm
        imsimband = Bandpass()
        imsimband.imsimBandpass()

        outputNames=[]

        for (sedName, zz, iAv, iRv, gAv, gRv, norm) in \
            zip(actualSEDnames, redshift, internalAv, internalRv, galacticAv, galacticRv, magNorm):

            if is_null(sedName):
                sedList.append(None)
            else:
                if sedName in self.uniqueSeds:
                    #we have already read in this file; no need to do it again
                    sed = Sed(wavelen=self.uniqueSeds[sedName].wavelen,
                              flambda=self.uniqueSeds[sedName].flambda,
                              fnu=self.uniqueSeds[sedName].fnu,
                              name=self.uniqueSeds[sedName].name)
                else:
                    #load the SED of the object
                    sed = Sed()
                    sedFile = os.path.join(self.sedDir, sedName)
                    sed.readSED_flambda(sedFile)

                    flambdaCopy = copy.deepcopy(sed.flambda)

                    #If the SED is zero inside of the bandpass, GalSim raises an error.
                    #This sets a minimum flux value of 1.0e-30 so that the SED is never technically
                    #zero inside of the bandpass.
                    sed.flambda = numpy.array([ff if ff>1.0e-30 else 1.0e-30 for ff in flambdaCopy])
                    sed.fnu = None

                    #copy the unnormalized file to uniqueSeds so we don't have to read it in again
                    sedCopy = Sed(wavelen=sed.wavelen, flambda=sed.flambda,
                                  fnu=sed.fnu, name=sed.name)
                    self.uniqueSeds[sedName] = sedCopy

                #normalize the SED
                #Consulting the file sed.py in GalSim/galsim/ it appears that GalSim expects
                #its SEDs to ultimately be in units of ergs/nm so that, when called, they can
                #be converted to photons/nm (see the function __call__() and the assignment of
                #self._rest_photons in the __init__() of galsim's sed.py file).  Thus, we need
                #to read in our SEDs, normalize them, and then multiply by the exposure time
                #and the effective area to get from ergs/s/cm^2/nm to ergs/nm.
                #
                #The gain parameter should convert between photons and ADU (so: it is the
                #traditional definition of "gain" -- electrons per ADU -- multiplied by the
                #quantum efficiency of the detector).  Because we fold the quantum efficiency
                #of the detector into our total_[u,g,r,i,z,y].dat bandpass files
                #(see the readme in the THROUGHPUTS_DIR/baseline/), we only need to multiply
                #by the electrons per ADU gain.
                #
                #We will take these parameters from an instantiation of the PhotometricParameters
                #class (which can be reassigned by defining a daughter class of this class)
                #
                fNorm = sed.calcFluxNorm(norm, imsimband)
                sed.multiplyFluxNorm(fNorm)

                #apply dust extinction (internal)
                if iAv != 0.0 and iRv != 0.0:
                    a_int, b_int = sed.setupCCMab()
                    sed.addCCMDust(a_int, b_int, A_v=iAv, R_v=iRv)

                #22 June 2015
                #apply redshift; there is no need to apply the distance modulus from
                #sims/photUtils/CosmologyWrapper; magNorm takes that into account
                #however, magNorm does not take into account cosmological dimming
                if zz != 0.0:
                    sed.redshiftSED(zz, dimming=True)

                #apply dust extinction (galactic)
                a_int, b_int = sed.setupCCMab()
                sed.addCCMDust(a_int, b_int, A_v=gAv, R_v=gRv)
                sedList.append(sed)

        return sedList

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
    def testAlternateNormalizingBandpass(self):
        """
        A reiteration of testAddingToList, but testing with a non-imsimBandpass
        normalizing bandpass
        """
        normalizingBand = Bandpass()
        normalizingBand.readThroughput(os.path.join(getPackageDir('throughputs'),'baseline','total_r.dat'))
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0,
                           normalizingBandpass=normalizingBand,
                           internalAvList=internalAvList_0,
                           redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                           wavelenMatch=wavelen_match)


        sedNameList_1 = self.getListOfSedNames(nSed)
        magNormList_1 = numpy.random.random_sample(nSed)*5.0 + 15.0

        internalAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1

        redshiftList_1 = numpy.random.random_sample(nSed)*5.0

        galacticAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1


        testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                  internalAvList=internalAvList_1,
                                  galacticAvList=galacticAvList_1,
                                  redshiftList=redshiftList_1)

        self.assertEqual(len(testList), 2*nSed)
        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

        for ix in range(len(sedNameList_0)):
            self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
            self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
            self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)


        for ix in range(len(sedNameList_1)):
            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix+nSed], 10)
            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix+nSed], 10)
            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix+nSed], 10)

        for ix, (name, norm, iav, gav, zz) in \
          enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                     galacticAvList_0, redshiftList_0)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, normalizingBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)

        for ix, (name, norm, iav, gav, zz) in \
            enumerate(zip(sedNameList_1, magNormList_1, internalAvList_1, \
                          galacticAvList_1, redshiftList_1)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, normalizingBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)

            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix+nSed]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import redshiftSED [as 别名]
def run_lsst_baker(redshift_grid):
    """
    Accepts a redshift grid (a numpy array).  Reads in all of the SEDs from
    ../data/sed/ and the calculates their colors in all of the lsst bands,
    given those redshifts

    Returns a structured numpy array containing
    sedname
    redshift
    ug
    gr
    ri
    iz
    zy
    time (the time spent calculating the magnitudes)
    """
    bandpass_dir = os.path.join('..', 'data', 'bandpass')
    bandpass_list = []
    for bp_name in ('total_u.dat', 'total_g.dat', 'total_r.dat', 'total_i.dat',
                    'total_z.dat', 'total_y.dat'):

        bp = Bandpass()
        bp.readThroughput(os.path.join(bandpass_dir, bp_name))
        bandpass_list.append(bp)

    sed_dir = os.path.join('..', 'data', 'sed')
    list_of_sed_names = sorted([nn for nn in os.listdir(sed_dir)
                                if 'spec' in nn and 'ang' not in nn])

    n_rows = len(list_of_sed_names)*len(redshift_grid)
    dtype = np.dtype([('sedname', str, 300), ('redshift', np.float),
                      ('ug', np.float), ('gr', np.float), ('ri', np.float),
                      ('iz', np.float), ('zy', np.float), ('time', np.float)])

    dummy = ('aaaa', 1.0, 1, 1, 1, 1, 1, 1.0)
    output_array = np.array([dummy]*n_rows, dtype=dtype)

    i_row = 0
    for redshift in redshift_grid:
        for sed_name in list_of_sed_names:
            ss = Sed()
#            t_start = time.clock()
            t_start = time.time()
            ss.readSED_flambda(os.path.join(sed_dir, sed_name))
            ss.redshiftSED(redshift)
            local_list = []
            for bp in bandpass_list:
                mm = ss.calcMag(bp)
                local_list.append(mm)
            #print local_list[0], local_list[1], local_list[0]-local_list[1]
#            time_spent = time.clock()-t_start
            time_spent = time.time()-t_start
            output_array[i_row][0] = sed_name
            output_array[i_row][1] = redshift
            output_array[i_row][2] = local_list[0]-local_list[1]
            output_array[i_row][3] = local_list[1]-local_list[2]
            output_array[i_row][4] = local_list[2]-local_list[3]
            output_array[i_row][5] = local_list[3]-local_list[4]
            output_array[i_row][6] = local_list[4]-local_list[5]
            output_array[i_row][7] = time_spent
            i_row += 1

    return output_array