Python Sed.addDust方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import addDust [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 addDust [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 addDust [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 addDust [as 别名]
    def test_mixed_stars(self):
        """
        Here we will test the (somewhat absurd) case of a catalog with two different bandpasses
        (lsst_ and cartoon_) in order to verify that gamma values are being cached correctly
        """

        lsst_u_band = Bandpass()
        lsst_u_band.readThroughput(os.path.join(getPackageDir('throughputs'), 'baseline', 'total_u.dat'))
        lsst_g_band = Bandpass()
        lsst_g_band.readThroughput(os.path.join(getPackageDir('throughputs'), 'baseline', 'total_g.dat'))

        obs = ObservationMetaData(bandpassName=['c_u', 'c_g', 'u', 'g'],
                                  m5=[25.0, 26.0, 15.0, 16.0])
        # make the difference in m5 between the two bandpass systems extreme
        # so that, in the unit test, we can be sure that the correct values
        # are being used for the correct getters

        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, column_outputs=['lsst_u', 'lsst_g',
                                                                 'sigma_lsst_u', 'sigma_lsst_g'])
        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)
            lsst_umag = spectrum.calcMag(lsst_u_band)
            self.assertAlmostEqual(lsst_umag, controlData['lsst_u'][ix], 3)
            lsst_gmag = spectrum.calcMag(lsst_g_band)
            self.assertAlmostEqual(lsst_gmag, controlData['lsst_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)

            lsst_umagError, gamma = calcMagError_m5(lsst_umag, lsst_u_band,
                                                    obs.m5['u'], PhotometricParameters())
            lsst_gmagError, gamma = calcMagError_m5(lsst_gmag, lsst_g_band,
                                                    obs.m5['g'], PhotometricParameters())

            self.assertAlmostEqual(lsst_umagError, controlData['sigma_lsst_u'][ix], 3)
            self.assertAlmostEqual(lsst_gmagError, controlData['sigma_lsst_g'][ix], 3)
            self.assertGreater(np.abs(lsst_umagError-umagError), 0.01)
            self.assertGreater(np.abs(lsst_gmagError-gmagError), 0.01)