Python Sed.calcMag方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
 def _calcColors(self, sedname='C.dat'):
     """
     Calculate the colors for a moving object with sed 'sedname'.
     """
     # Do we need to read in the LSST bandpasses?
     try:
         self.lsst
     except AttributeError:
         filterdir = os.getenv('LSST_THROUGHPUTS_BASELINE')
         filterlist = ('u', 'g', 'r', 'i', 'z', 'y')
         self.lsst ={}
         for f in filterlist:
             self.lsst[f] = Bandpass()
             self.lsst[f].readThroughput(os.path.join(filterdir, 'total_'+f+'.dat'))
         self.vband = Bandpass()
         self.vband.readThroughput('harris_V.dat')
         self.colors = {}
     # See if the sed's colors are in memory already.
     if sedname not in self.colors:
         moSed = Sed()
         moSed.readSED_flambda(sedname)
         vmag = moSed.calcMag(self.vband)
         self.colors[sedname] = {}
         for f in filterlist:
             self.colors[sedname][f] = moSed.calcMag(self.lsst[f]) - vmag
     return self.colors[sedname]

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def testMagListForSedList(self):
        """
        Test that magListForSedList calculates the correct magnitude
        """

        nBandpasses = 7
        bpNameList, bpList = self.getListOfBandpasses(nBandpasses)
        testBpDict = BandpassDict(bpList, bpNameList)

        nSed = 20
        sedNameList = self.getListOfSedNames(nSed)
        magNormList = self.rng.random_sample(nSed)*5.0 + 15.0
        internalAvList = self.rng.random_sample(nSed)*0.3 + 0.1
        redshiftList = self.rng.random_sample(nSed)*5.0
        galacticAvList = self.rng.random_sample(nSed)*0.3 + 0.1

        # first, test on an SedList without a wavelenMatch
        testSedList = SedList(sedNameList, magNormList,
                              fileDir=self.sedDir,
                              internalAvList=internalAvList,
                              redshiftList=redshiftList,
                              galacticAvList=galacticAvList)

        magList = testBpDict.magListForSedList(testSedList)
        self.assertEqual(magList.shape[0], nSed)
        self.assertEqual(magList.shape[1], nBandpasses)

        for ix, sedObj in enumerate(testSedList):
            dummySed = Sed(wavelen=copy.deepcopy(sedObj.wavelen),
                           flambda=copy.deepcopy(sedObj.flambda))

            for iy, bp in enumerate(testBpDict):
                mag = dummySed.calcMag(bpList[iy])
                self.assertAlmostEqual(mag, magList[ix][iy], 2)

        # now use wavelenMatch
        testSedList = SedList(sedNameList, magNormList,
                              fileDir=self.sedDir,
                              internalAvList=internalAvList,
                              redshiftList=redshiftList,
                              galacticAvList=galacticAvList,
                              wavelenMatch=testBpDict.wavelenMatch)

        magList = testBpDict.magListForSedList(testSedList)
        self.assertEqual(magList.shape[0], nSed)
        self.assertEqual(magList.shape[1], nBandpasses)

        for ix, sedObj in enumerate(testSedList):
            dummySed = Sed(wavelen=copy.deepcopy(sedObj.wavelen),
                           flambda=copy.deepcopy(sedObj.flambda))

            for iy, bp in enumerate(testBpDict):
                mag = dummySed.calcMag(bpList[iy])
                self.assertAlmostEqual(mag, magList[ix][iy], 2)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def returnMags(self, bandpass=None):
        """
        Convert the computed spectra to magnitudes using the supplied bandpasses,
        or, if self.mags=True, just return the mags in the LSST filters

        If mags=True when initialized, return mags returns an structured array with
        dtype names u,g,r,i,z,y.
        """
        if self.mags:
            if bandpass:
                warnings.warn('Ignoring set bandpasses and returning LSST ugrizy.')
            mags = -2.5*np.log10(self.spec)+np.log10(3631.)
            # Mask out high airmass
            mags[self.mask] *= np.nan
            # Convert to a structured array
            mags = np.core.records.fromarrays(mags.transpose(),
                                              names='u,g,r,i,z,y',
                                              formats='float,'*6)
        else:
            mags = np.zeros(self.npts, dtype=float)-666
            tempSed = Sed()
            isThrough = np.where(bandpass.sb > 0)
            minWave = bandpass.wavelen[isThrough].min()
            maxWave = bandpass.wavelen[isThrough].max()
            inBand = np.where((self.wave >= minWave) & (self.wave <= maxWave))
            for i, ra in enumerate(self.ra):
                # Check that there is flux in the band, otherwise calcMag fails
                if np.max(self.spec[i, inBand]) > 0:
                    tempSed.setSED(self.wave, flambda=self.spec[i, :])
                    mags[i] = tempSed.calcMag(bandpass)

            # Mask out high airmass
            mags[self.mask] *= np.nan
        return mags

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
def calcColors(sedname='C.dat'):
    # Calculate SSO colors.
    filterdir = os.getenv('LSST_THROUGHPUTS_BASELINE')
    filterlist = ('u', 'g', 'r', 'i', 'z', 'y')
    lsst ={}
    for f in filterlist:
        lsst[f] = Bandpass()
        lsst[f].readThroughput(os.path.join(filterdir, 'total_'+f+'.dat'))
    vband = Bandpass()
    vband.readThroughput('harris_V.dat')
    csed = Sed()
    csed.readSED_flambda(sedname)
    vmag = csed.calcMag(vband)
    dmags = {}
    for f in filterlist:
        dmags[f] = csed.calcMag(lsst[f]) - vmag
    return dmags

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
 def get_mags(self,source, phase ):
     
     sed = Sed()
     print "==========================================="
     print phase
     print "==========================================="
     if phase > -20 and phase < 50 :
         sourceflux = source.flux(phase=phase, wave=self.rband.wavelen*10.)
         sed.setSED(wavelen=self.rband.wavelen, flambda=sourceflux/10.)
     else:
         sed.setSED(wavelen=self.rband.wavelen, flambda=flambda)
     #sed.redshiftSED(redshift=_z[i], dimming=True)
     return [sed.calcMag(bandpass=self.uband),
             sed.calcMag(bandpass=self.gband),
             sed.calcMag(bandpass=self.rband),
             sed.calcMag(bandpass=self.iband),
             sed.calcMag(bandpass=self.zband)]

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def testSignalToNoise(self):
        """
        Test that calcSNR_m5 and calcSNR_sed give similar results
        """
        defaults = LSSTdefaults()
        photParams = PhotometricParameters()

        m5 = []
        for i in range(len(self.hardwareList)):
            m5.append(
                snr.calcM5(
                    self.skySed,
                    self.bpList[i],
                    self.hardwareList[i],
                    photParams,
                    seeing=defaults.seeing(self.filterNameList[i]),
                )
            )

        sedDir = lsst.utils.getPackageDir("sims_sed_library")
        sedDir = os.path.join(sedDir, "starSED", "kurucz")
        fileNameList = os.listdir(sedDir)

        numpy.random.seed(42)
        offset = numpy.random.random_sample(len(fileNameList)) * 2.0

        for ix, name in enumerate(fileNameList):
            if ix > 100:
                break
            spectrum = Sed()
            spectrum.readSED_flambda(os.path.join(sedDir, name))
            ff = spectrum.calcFluxNorm(m5[2] - offset[ix], self.bpList[2])
            spectrum.multiplyFluxNorm(ff)
            magList = []
            controlList = []
            magList = []
            for i in range(len(self.bpList)):
                controlList.append(
                    snr.calcSNR_sed(
                        spectrum,
                        self.bpList[i],
                        self.skySed,
                        self.hardwareList[i],
                        photParams,
                        defaults.seeing(self.filterNameList[i]),
                    )
                )

                magList.append(spectrum.calcMag(self.bpList[i]))

            testList, gammaList = snr.calcSNR_m5(
                numpy.array(magList), numpy.array(self.bpList), numpy.array(m5), photParams
            )

            for tt, cc in zip(controlList, testList):
                msg = "%e != %e " % (tt, cc)
                self.assertTrue(numpy.abs(tt / cc - 1.0) < 0.001, msg=msg)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [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 calcMag [as 别名]
    def test_stars(self):
        obs = ObservationMetaData(bandpassName=['c_u', 'c_g'], m5=[25.0, 26.0])

        db_dtype = np.dtype([('id', np.int),
                             ('raJ2000', np.float),
                             ('decJ2000', np.float),
                             ('sedFilename', str, 100),
                             ('magNorm', np.float),
                             ('galacticAv', np.float)])

        inputDir = os.path.join(getPackageDir('sims_catUtils'), 'tests', 'testData')
        inputFile = os.path.join(inputDir, 'IndicesTestCatalogStars.txt')
        db = fileDBObject(inputFile, dtype=db_dtype, runtable='test', idColKey='id')
        cat = CartoonStars(db, obs_metadata=obs)
        with lsst.utils.tests.getTempFilePath('.txt') as catName:
            cat.write_catalog(catName)
            dtype = np.dtype([(name, np.float) for name in cat.column_outputs])
            controlData = np.genfromtxt(catName, dtype=dtype, delimiter=',')

        db_columns = db.query_columns(['id', 'raJ2000', 'decJ2000', 'sedFilename', 'magNorm', 'galacticAv'])

        sedDir = os.path.join(getPackageDir('sims_sed_library'), 'starSED', 'kurucz')

        for ix, line in enumerate(next(db_columns)):
            spectrum = Sed()
            spectrum.readSED_flambda(os.path.join(sedDir, line[3]))
            fnorm = spectrum.calcFluxNorm(line[4], self.normband)
            spectrum.multiplyFluxNorm(fnorm)
            a_x, b_x = spectrum.setupCCM_ab()
            spectrum.addDust(a_x, b_x, A_v=line[5])
            umag = spectrum.calcMag(self.uband)
            self.assertAlmostEqual(umag, controlData['cartoon_u'][ix], 3)
            gmag = spectrum.calcMag(self.gband)
            self.assertAlmostEqual(gmag, controlData['cartoon_g'][ix], 3)
            umagError, gamma = calcMagError_m5(umag, self.uband, obs.m5['c_u'], PhotometricParameters())
            gmagError, gamma = calcMagError_m5(gmag, self.gband, obs.m5['c_g'], PhotometricParameters())
            self.assertAlmostEqual(umagError, controlData['sigma_cartoon_u'][ix], 3)
            self.assertAlmostEqual(gmagError, controlData['sigma_cartoon_g'][ix], 3)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def test_mags_vs_flux(self):
        """
        Verify that the relationship between Sed.calcMag() and Sed.calcFlux()
        is as expected
        """
        wavelen = np.arange(100.0, 1500.0, 1.0)
        flambda = np.exp(-0.5*np.power((wavelen-500.0)/100.0,2))
        sb = (wavelen-100.0)/1400.0

        ss = Sed(wavelen=wavelen, flambda=flambda)
        bp = Bandpass(wavelen=wavelen, sb=sb)

        mag = ss.calcMag(bp)
        flux = ss.calcFlux(bp)

        self.assertAlmostEqual(ss.magFromFlux(flux)/mag, 1.0, 10)
        self.assertAlmostEqual(ss.fluxFromMag(mag)/flux, 1.0, 10)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def testMagDictForSed(self):
        """
        Test that magDictForSed calculates the correct magnitude
        """

        wavelen = numpy.arange(10.0,2000.0,1.0)
        flux = (wavelen*2.0-5.0)*1.0e-6
        spectrum = Sed(wavelen=wavelen, flambda=flux)

        for nBp in range(3, 10, 1):

            nameList, bpList = self.getListOfBandpasses(nBp)
            testDict = BandpassDict(bpList, nameList)
            self.assertFalse(len(testDict.values()[0].wavelen)==len(spectrum.wavelen))

            magDict = testDict.magDictForSed(spectrum)
            for ix, (name, bp) in enumerate(zip(nameList, bpList)):
                magControl = spectrum.calcMag(bp)
                self.assertAlmostEqual(magDict[name], magControl, 5)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def computeMags(self, bandpass=None):
        """After the spectra have been computed, optionally convert to mags"""
        if self.mags:
            mags = -2.5*np.log10(self.spec)+np.log10(3631.)
        else:
            mags = np.zeros(self.npts, dtype=float)-666
            tempSed = Sed()
            isThrough = np.where(bandpass.sb > 0)
            minWave = bandpass.wavelen[isThrough].min()
            maxWave = bandpass.wavelen[isThrough].max()
            inBand = np.where( (self.wave >= minWave) & (self.wave <= maxWave))
            for i, ra in enumerate(self.ra):
                if np.max(self.spec[i,inBand]) > 0:
                    tempSed.setSED(self.wave, flambda=self.spec[i,:])
                    # Need to try/except because the spectra might be zero in the filter
                    # XXX-upgrade this to check if it's zero
                    mags[i] = tempSed.calcMag(bandpass)

        return mags

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def returnMags(self, bandpasses=None):
        """
        Convert the computed spectra to a magnitude using the supplied bandpass,
        or, if self.mags=True, return the mags in the LSST filters

        If mags=True when initialized, return mags returns an structured array with
        dtype names u,g,r,i,z,y.

        bandpasses: optional dictionary with bandpass name keys and bandpass object values.

        """
        if self.azs is None:
            raise ValueError('No coordinates set. Use setRaDecMjd, setRaDecAltAzMjd, or setParams methods before calling returnMags.')

        if self.mags:
            if bandpasses:
                warnings.warn('Ignoring set bandpasses and returning LSST ugrizy.')
            mags = -2.5*np.log10(self.spec)+np.log10(3631.)
            # Mask out high airmass
            mags[self.mask] *= np.nan
            mags = mags.swapaxes(0, 1)
            magsBack = {}
            for i, f in enumerate(self.filterNames):
                magsBack[f] = mags[i]
        else:
            magsBack = {}
            for key in bandpasses:
                mags = np.zeros(self.npts, dtype=float)-666
                tempSed = Sed()
                isThrough = np.where(bandpasses[key].sb > 0)
                minWave = bandpasses[key].wavelen[isThrough].min()
                maxWave = bandpasses[key].wavelen[isThrough].max()
                inBand = np.where((self.wave >= minWave) & (self.wave <= maxWave))
                for i, ra in enumerate(self.ra):
                    # Check that there is flux in the band, otherwise calcMag fails
                    if np.max(self.spec[i, inBand]) > 0:
                        tempSed.setSED(self.wave, flambda=self.spec[i, :])
                        mags[i] = tempSed.calcMag(bandpasses[key])
                # Mask out high airmass
                mags[self.mask] *= np.nan
                magsBack[key] = mags
        return magsBack

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
def calcWDColors():
    """
    Calculate a few example WD colors. Values to go in stellarMags(). Here in case
    values need to be regenerated (different stars, bandpasses change, etc.)
    """

    try:
        from lsst.utils import getPackageDir
        import os
        from lsst.sims.photUtils import Bandpass, Sed
    except:
        'Need to setup sims_photUtils to generate WD magnitudes.'

    names = ['HeWD_25200_80', 'WD_11000_85', 'WD_3000_85']
    fns = ['bergeron_He_24000_80.dat_25200.gz',
           'bergeron_10500_85.dat_11000.gz', 'bergeron_2750_85.dat_3000.gz']
    wdDir = os.path.join(getPackageDir('sims_sed_library'), 'starSED/wDs/')
    files = [os.path.join(wdDir, filename) for filename in fns]

    # Read in the LSST bandpasses
    bpNames = ['u', 'g', 'r', 'i', 'z', 'y']
    bps = []
    throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
    for key in bpNames:
        bp = np.loadtxt(os.path.join(throughPath, 'filter_' + key + '.dat'),
                        dtype=list(zip(['wave', 'trans'], [float] * 2)))
        tempB = Bandpass()
        tempB.setBandpass(bp['wave'], bp['trans'])
        bps.append(tempB)

    # Read in the SEDs and compute mags
    mags = []
    for filename in files:
        star = Sed()
        star.readSED_flambda(filename)
        singleMags = [star.calcMag(band) for band in bps]
        mags.append([singleMags[i - 1] - singleMags[i] for i in range(1, 6)])

    for maglist, fn, name in zip(mags, fns, names):
        format = (name, fn) + tuple(maglist)
        print("['%s', '%s', %f, %f, %f, %f, %f]" % format)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def testSignalToNoise(self):
        """
        Test that calcSNR_m5 and calcSNR_sed give similar results
        """
        defaults = LSSTdefaults()
        photParams = PhotometricParameters()

        m5 = []
        for i in range(len(self.hardwareList)):
            m5.append(snr.calcM5(self.skySed, self.bpList[i],
                      self.hardwareList[i],
                      photParams, FWHMeff=defaults.FWHMeff(self.filterNameList[i])))

        sedDir = os.path.join(lsst.utils.getPackageDir('sims_photUtils'),
                              'tests/cartoonSedTestData/starSed/')
        sedDir = os.path.join(sedDir, 'kurucz')
        fileNameList = os.listdir(sedDir)

        rng = np.random.RandomState(42)
        offset = rng.random_sample(len(fileNameList))*2.0

        for ix, name in enumerate(fileNameList):
            if ix > 100:
                break
            spectrum = Sed()
            spectrum.readSED_flambda(os.path.join(sedDir, name))
            ff = spectrum.calcFluxNorm(m5[2]-offset[ix], self.bpList[2])
            spectrum.multiplyFluxNorm(ff)
            for i in range(len(self.bpList)):
                control_snr = snr.calcSNR_sed(spectrum, self.bpList[i],
                                              self.skySed,
                                              self.hardwareList[i],
                                              photParams, defaults.FWHMeff(self.filterNameList[i]))

                mag = spectrum.calcMag(self.bpList[i])

                test_snr, gamma = snr.calcSNR_m5(mag, self.bpList[i], m5[i], photParams)
                self.assertLess((test_snr-control_snr)/control_snr, 0.001)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcMag [as 别名]
    def testMagError(self):
        """
        Make sure that calcMagError_sed and calcMagError_m5
        agree to within 0.001
        """
        defaults = LSSTdefaults()
        photParams = PhotometricParameters()

        # create a cartoon spectrum to test on
        spectrum = Sed()
        spectrum.setFlatSED()
        spectrum.multiplyFluxNorm(1.0e-9)

        # find the magnitudes of that spectrum in our bandpasses
        magList = []
        for total in self.bpList:
            magList.append(spectrum.calcMag(total))
        magList = numpy.array(magList)

        # try for different normalizations of the skySED
        for fNorm in numpy.arange(1.0, 5.0, 1.0):
            self.skySed.multiplyFluxNorm(fNorm)
            m5List = []
            magSed = []
            for total, hardware, filterName in zip(self.bpList, self.hardwareList, self.filterNameList):

                seeing = defaults.seeing(filterName)

                m5List.append(snr.calcM5(self.skySed, total, hardware, photParams, seeing=seeing))

                magSed.append(snr.calcMagError_sed(spectrum, total, self.skySed, hardware, photParams, seeing=seeing))

            magSed = numpy.array(magSed)

            magM5 = snr.calcMagError_m5(magList, self.bpList, numpy.array(m5List), photParams)

            numpy.testing.assert_array_almost_equal(magM5, magSed, decimal=3)