Python Sed.calcFlux方法代码示例

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

        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)

        fluxList = testBpDict.fluxListForSedList(testSedList)
        self.assertEqual(fluxList.shape[0], nSed)
        self.assertEqual(fluxList.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):
                flux = dummySed.calcFlux(bpList[iy])
                self.assertAlmostEqual(flux/fluxList[ix][iy], 1.0, 2)

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

        fluxList = testBpDict.fluxListForSedList(testSedList)
        self.assertEqual(fluxList.shape[0], nSed)
        self.assertEqual(fluxList.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):
                flux = dummySed.calcFlux(bpList[iy])
                self.assertAlmostEqual(flux/fluxList[ix][iy], 1.0, 2)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcFlux [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 calcFlux [as 别名]
    def testFluxDictForSed(self):
        """
        Test that fluxDictForSed calculates the correct fluxes
        """

        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))

            fluxDict = testDict.fluxDictForSed(spectrum)
            for ix, (name, bp) in enumerate(zip(nameList, bpList)):
                fluxControl = spectrum.calcFlux(bp)
                self.assertAlmostEqual(fluxDict[name]/fluxControl, 1.0, 2)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcFlux [as 别名]
 def testSedMagErrors(self):
     """Test error handling at mag and adu calculation levels of sed."""
     sedwavelen = np.arange(self.wmin+50, self.wmax, 1)
     sedflambda = np.ones(len(sedwavelen))
     testsed = Sed(wavelen=sedwavelen, flambda=sedflambda)
     # Test handling in calcMag
     with warnings.catch_warnings(record=True) as w:
         mag = testsed.calcMag(self.testbandpass)
         self.assertEqual(len(w), 1)
         self.assertIn("non-overlap", str(w[-1].message))
     np.testing.assert_equal(mag, np.NaN)
     # Test handling in calcADU
     with warnings.catch_warnings(record=True) as w:
         adu = testsed.calcADU(self.testbandpass,
                               photParams=PhotometricParameters())
         self.assertEqual(len(w), 1)
         self.assertIn("non-overlap", str(w[-1].message))
     np.testing.assert_equal(adu, np.NaN)
     # Test handling in calcFlux
     with warnings.catch_warnings(record=True) as w:
         flux = testsed.calcFlux(self.testbandpass)
         self.assertEqual(len(w), 1)
         self.assertIn("non-overlap", str(w[-1].message))
     np.testing.assert_equal(flux, np.NaN)

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

#.........这里部分代码省略.........
        # milli-mag, our flaring model gives us the magnitudes
        # expected, given the light curves specified in
        # setUpClass()
        for mjd in (59580.0, 60000.0, 70000.0, 80000.0):

            obs = ObservationMetaData(mjd=mjd)

            quiet_cat = QuiescentCatalog(db, obs_metadata=obs)
            quiet_cat.write_catalog(quiet_cat_name)

            flare_cat = FlaringCatalogDummy(db, obs_metadata=obs)
            flare_cat.scratch_dir = self.scratch_dir
            flare_cat._mlt_lc_file = self.mlt_lc_name
            flare_cat.write_catalog(flare_cat_name)

            quiescent_data = np.genfromtxt(quiet_cat_name, dtype=dtype, delimiter=',')
            flaring_data = np.genfromtxt(flare_cat_name, dtype=dtype, delimiter=',')

            self.assertGreater(len(quiescent_data), 2)
            self.assertEqual(len(quiescent_data), len(flaring_data))
            self.assertIn(3, flaring_data['id'])

            for ix in range(len(flaring_data)):
                obj_id = flaring_data['id'][ix]
                self.assertEqual(obj_id, ix)


                msg = ('failed on object %d; mjd %.2f\n u_quiet %e u_flare %e\n g_quiet %e g_flare %e' %
                       (obj_id, mjd, quiescent_data['u'][obj_id], flaring_data['u'][obj_id],
                        quiescent_data['g'][obj_id], flaring_data['g'][obj_id]))

                self.assertEqual(quiescent_data['id'][obj_id], flaring_data['id'][obj_id], msg=msg)
                self.assertAlmostEqual(ss.magFromFlux(baseline_fluxes[obj_id][0]),
                                       quiescent_data['u'][obj_id], 3, msg=msg)
                self.assertAlmostEqual(ss.magFromFlux(baseline_fluxes[obj_id][1]),
                                       quiescent_data['g'][obj_id], 3, msg=msg)
                if obj_id != 3:

                    # the models below are as specified in the
                    # setUpClass() method
                    if obj_id == 0 or obj_id == 1:
                        amp = 1.0e42
                        dt = 3652.5
                        t_min = flare_cat._survey_start - t0_list[obj_id]

                        tt = mjd - t_min
                        while tt > dt:
                            tt -= dt

                        u_flux = amp*(1.0+np.power(np.sin(tt/100.0), 2))
                        g_flux = amp*(1.0+np.power(np.cos(tt/100.0), 2))
                    elif obj_id==2:
                        amp = 2.0e41
                        dt = 365.25
                        t_min = flare_cat._survey_start - t0_list[obj_id]

                        tt = mjd - t_min
                        while tt > dt:
                            tt -= dt
                        u_flux = amp*(1.0+np.power(np.sin(tt/50.0), 2))
                        g_flux = amp*(1.0+np.power(np.cos(tt/50.0), 2))

                    # calculate the multiplicative effect of dust on a 9000K
                    # black body
                    bb_sed = Sed(wavelen=bb_wavelen, flambda=bb_flambda)
                    u_bb_flux = bb_sed.calcFlux(bp_dict['u'])
                    g_bb_flux = bb_sed.calcFlux(bp_dict['g'])
                    a_x, b_x = bb_sed.setupCCM_ab()
                    bb_sed.addDust(a_x, b_x, A_v=av_list[obj_id])
                    u_bb_dusty_flux = bb_sed.calcFlux(bp_dict['u'])
                    g_bb_dusty_flux = bb_sed.calcFlux(bp_dict['g'])

                    dust_u = u_bb_dusty_flux/u_bb_flux
                    dust_g = g_bb_dusty_flux/g_bb_flux

                    area = 4.0*np.pi*np.power(distance_list[obj_id], 2)
                    tot_u_flux = baseline_fluxes[obj_id][0] + u_flux*dust_u/area
                    tot_g_flux = baseline_fluxes[obj_id][1] + g_flux*dust_g/area

                    self.assertAlmostEqual(ss.magFromFlux(tot_u_flux), flaring_data['u'][obj_id],
                                           3, msg=msg)
                    self.assertAlmostEqual(ss.magFromFlux(tot_g_flux), flaring_data['g'][obj_id],
                                           3, msg=msg)

                    self.assertGreater(np.abs(flaring_data['g'][obj_id]-quiescent_data['g'][obj_id]),
                                       0.001, msg=msg)
                    self.assertGreater(np.abs(flaring_data['u'][obj_id]-quiescent_data['u'][obj_id]),
                                       0.001, msg=msg)
                else:
                    self.assertAlmostEqual(flaring_data['g'][obj_id],
                                           quiescent_data['g'][obj_id]+3*(mjd-59580.0)/10000.0,
                                           3, msg=msg)
                    self.assertAlmostEqual(flaring_data['u'][obj_id],
                                           quiescent_data['u'][obj_id]+2*(mjd-59580.0)/10000.0,
                                           3, msg=msg)

        if os.path.exists(quiet_cat_name):
            os.unlink(quiet_cat_name)
        if os.path.exists(flare_cat_name):
            os.unlink(flare_cat_name)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcFlux [as 别名]
    def testIndicesOnFlux(self):
        """
        Test that, when you pass a list of indices into the calcFluxList
        methods, you get the correct fluxes out.
        """

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

        # first try it with a single Sed
        wavelen = numpy.arange(10.0,2000.0,1.0)
        flux = (wavelen*2.0-5.0)*1.0e-6
        spectrum = Sed(wavelen=wavelen, flambda=flux)
        indices = [1,2,5]

        fluxList = testBpDict.fluxListForSed(spectrum, indices=indices)
        ctNaN = 0
        for ix, (name, bp, fluxTest) in enumerate(zip(nameList, bpList, fluxList)):
            if ix in indices:
                fluxControl = spectrum.calcFlux(bp)
                self.assertAlmostEqual(fluxTest/fluxControl, 1.0, 2)
            else:
                ctNaN += 1
                self.assertTrue(numpy.isnan(fluxTest))

        self.assertEqual(ctNaN, 4)

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

        # now try a SedList without a wavelenMatch
        testSedList = SedList(sedNameList, magNormList,
                                    internalAvList=internalAvList,
                                    redshiftList=redshiftList,
                                    galacticAvList=galacticAvList)

        fluxList = testBpDict.fluxListForSedList(testSedList, indices=indices)
        self.assertEqual(fluxList.shape[0], nSed)
        self.assertEqual(fluxList.shape[1], nBandpasses)

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

            ctNaN = 0
            for iy, bp in enumerate(testBpDict):
                if iy in indices:
                    flux = dummySed.calcFlux(testBpDict[bp])
                    self.assertAlmostEqual(flux/fluxList[ix][iy], 1.0, 2)
                else:
                    ctNaN += 1
                    self.assertTrue(numpy.isnan(fluxList[ix][iy]))

            self.assertEqual(ctNaN, 4)

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

        fluxList = testBpDict.fluxListForSedList(testSedList, indices=indices)
        self.assertEqual(fluxList.shape[0], nSed)
        self.assertEqual(fluxList.shape[1], nBandpasses)

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

            ctNaN = 0
            for iy, bp in enumerate(testBpDict):
                if iy in indices:
                    flux = dummySed.calcFlux(testBpDict[bp])
                    self.assertAlmostEqual(flux/fluxList[ix][iy], 1.0, 2)
                else:
                    ctNaN +=  1
                    self.assertTrue(numpy.isnan(fluxList[ix][iy]))

            self.assertEqual(ctNaN, 4)