Python Sed.multiplyFluxNorm方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [as 别名]
    def testStellarPhotometricUncertainties(self):
        """
        Test in the case of a catalog of stars
        """
        lsstDefaults = LSSTdefaults()
        starDB = testStarsDBObj(driver=self.driver, host=self.host, database=self.dbName)
        starCat = testStarCatalog(starDB, obs_metadata=self.obs_metadata)
        phot = PhotometryStars()

        ct = 0
        for line in starCat.iter_catalog():
            starSed = Sed()
            starSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                                 defaultSpecMap[line[14]]))
            imsimband = Bandpass()
            imsimband.imsimBandpass()
            fNorm = starSed.calcFluxNorm(line[15], imsimband)
            starSed.multiplyFluxNorm(fNorm)

            aV = numpy.float(line[16])
            a_int, b_int = starSed.setupCCMab()
            starSed.addCCMDust(a_int, b_int, A_v=aV)

            for i in range(len(self.bandpasses)):
                controlSigma = calcMagError_sed(starSed, self.totalBandpasses[i],
                                             self.skySeds[i],
                                             self.hardwareBandpasses[i],
                                             FWHMeff=lsstDefaults.FWHMeff(self.bandpasses[i]),
                                             photParams=PhotometricParameters())

                testSigma = line[8+i]
                self.assertAlmostEqual(controlSigma, testSigma, 4)
                ct += 1
        self.assertGreater(ct, 0)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [as 别名]
def setM5(m5target, skysed, totalBandpass, hardware,
          photParams,
          FWHMeff = None):
    """
    Take an SED representing the sky and normalize it so that
    m5 (the magnitude at which an object is detected in this
    bandpass at 5-sigma) is set to some specified value.

    The 5-sigma limiting magnitude (m5) for an observation is
    determined by a combination of the telescope and camera parameters
    (such as diameter of the mirrors and the readnoise) together with the
    sky background. This method (setM5) scales a provided sky background
    Sed so that an observation would have a target m5 value, for the
    provided hardware parameters. Using the resulting Sed in the
    'calcM5' method will return this target value for m5.

    @param [in] the desired value of m5

    @param [in] skysed is an instantiation of the Sed class representing
    sky emission

    @param [in] totalBandpass is an instantiation of the Bandpass class
    representing the total throughput of the telescope (instrumentation
    plus atmosphere)

    @param [in] hardware is an instantiation of the Bandpass class representing
    the throughput due solely to instrumentation.

    @param [in] photParams is an instantiation of the
    PhotometricParameters class that carries details about the
    photometric response of the telescope.

    @param [in] FWHMeff in arcseconds

    @param [out] returns an instantiation of the Sed class that is the skysed renormalized
    so that m5 has the desired value.

    Note that the returned SED will be renormalized such that calling the method
    self.calcADU(hardwareBandpass) on it will yield the number of counts per square
    arcsecond in a given bandpass.
    """

    #This is based on the LSST SNR document (v1.2, May 2010)
    #www.astro.washington.edu/users/ivezic/Astr511/LSST_SNRdoc.pdf

    if FWHMeff is None:
        FWHMeff = LSSTdefaults().FWHMeff('r')

    skyCountsTarget = calcSkyCountsPerPixelForM5(m5target, totalBandpass, FWHMeff=FWHMeff,
                                             photParams=photParams)

    skySedOut = Sed(wavelen=numpy.copy(skysed.wavelen),
                    flambda=numpy.copy(skysed.flambda))

    skyCounts = skySedOut.calcADU(hardware, photParams=photParams) \
                    * photParams.platescale * photParams.platescale
    skySedOut.multiplyFluxNorm(skyCountsTarget/skyCounts)

    return skySedOut

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [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 multiplyFluxNorm [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 multiplyFluxNorm [as 别名]
    def testVerboseSNR(self):
        """
        Make sure that calcSNR_sed has everything it needs to run in verbose mode
        """
        photParams = PhotometricParameters()

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

        snr.calcSNR_sed(spectrum, self.bpList[0], self.skySed,
                        self.hardwareList[0], photParams, FWHMeff=0.7, verbose=True)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [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 multiplyFluxNorm [as 别名]
def make_response_func(magnorm=16., filename='starSED/wDs/bergeron_14000_85.dat_14200.gz',
                       savefile='gaia_response.npz', noise=1, count_min=8.,
                       bluecut=700., redcut=650):
    """
    Declare some stars as "standards" and build a simple GAIA response function?

    Multiply GAIA observations by response function to get spectra in flambda units.
    """
    imsimBand = Bandpass()
    imsimBand.imsimBandpass()

    sed_dir = getPackageDir('sims_sed_library')
    filepath = os.path.join(sed_dir, filename)
    wd = Sed()
    wd.readSED_flambda(filepath)
    # Let's just use a flat spectrum
    wd.setFlatSED()
    fNorm = wd.calcFluxNorm(magnorm, imsimBand)
    wd.multiplyFluxNorm(fNorm)
    red_wd = copy.copy(wd)
    blue_wd = copy.copy(wd)
    gaia_obs = SED2GAIA(wd, noise=noise)
    red_wd.resampleSED(wavelen_match = gaia_obs['RP_wave'])
    blue_wd.resampleSED(wavelen_match = gaia_obs['BP_wave'])
    if noise == 1:
        red_response = red_wd.flambda / gaia_obs['noisySpec'][0]['RPNoisySpec']
        blue_response = blue_wd.flambda / gaia_obs['noisySpec'][0]['BPNoisySpec']
        too_low = np.where(gaia_obs['noisySpec'][0]['RPNoisySpec'] < count_min)
        red_response[too_low] = 0
        too_low = np.where(gaia_obs['noisySpec'][0]['BPNoisySpec'] < count_min)
        blue_response[too_low] = 0
    elif noise == 0:
        red_response = red_wd.flambda / gaia_obs['noiseFreeSpec']['RPNoiseFreeSpec']
        blue_response = blue_wd.flambda / gaia_obs['noiseFreeSpec']['BPNoiseFreeSpec']
        too_low = np.where(gaia_obs['noiseFreeSpec']['RPNoiseFreeSpec'] < count_min)
        red_response[too_low] = 0
        too_low = np.where(gaia_obs['noiseFreeSpec']['BPNoiseFreeSpec'] < count_min)
        blue_response[too_low] = 0

    blue_response[np.where(gaia_obs['BP_wave'] > bluecut)] = 0.
    red_response[np.where(gaia_obs['RP_wave'] < redcut)] = 0.

    # XXX check the mags of the original WD and the blue and red WD.

    np.savez(savefile, red_response=red_response, blue_response=blue_response,
             red_wavelen=gaia_obs['RP_wave'], blue_wavelen=gaia_obs['BP_wave'])

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [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 multiplyFluxNorm [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 multiplyFluxNorm [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 multiplyFluxNorm [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)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [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 = np.array(magList)

        # try for different normalizations of the skySED
        for fNorm in np.arange(1.0, 5.0, 1.0):
            self.skySed.multiplyFluxNorm(fNorm)

            for total, hardware, filterName, mm in \
                zip(self.bpList, self.hardwareList, self.filterNameList, magList):

                FWHMeff = defaults.FWHMeff(filterName)

                m5 = snr.calcM5(self.skySed, total, hardware, photParams, FWHMeff=FWHMeff)

                sigma_sed = snr.calcMagError_sed(spectrum, total, self.skySed,
                                                 hardware, photParams, FWHMeff=FWHMeff)

                sigma_m5, gamma = snr.calcMagError_m5(mm, total, m5, photParams)

                self.assertAlmostEqual(sigma_m5, sigma_sed, 3)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [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 multiplyFluxNorm [as 别名]
    def _quiescentMagnitudeGetter(self, bandpassDict, columnNameList, bandpassTag='lsst'):
        """
        Method that actually does the work calculating magnitudes for solar system objects.

        Because solar system objects have no dust extinction, this method works by loading
        each unique Sed once, normalizing it, calculating its magnitudes in the desired
        bandpasses, and then storing the normalizing magnitudes and the bandpass magnitudes
        in a dict.  Magnitudes for subsequent objects with identical Seds will be calculated
        by adding an offset to the magnitudes.  The offset is determined by comparing normalizing
        magnitues.

        @param [in] bandpassDict is an instantiation of BandpassDict representing the bandpasses
        to be integrated over

        @param [in] columnNameList is a list of the names of the columns being calculated
        by this getter

        @param [in] bandpassTag (optional) is a string indicating the name of the bandpass system
        (i.e. 'lsst', 'sdss', etc.).  This is in case the user wants to calculate the magnitudes
        in multiple systems simultaneously.  In that case, the dict will store magnitudes for each
        Sed in each magnitude system separately.

        @param [out] a numpy array of magnitudes corresponding to bandpassDict.
        """

        # figure out which of these columns we are actually calculating
        indices = [ii for ii, name in enumerate(columnNameList)
                   if name in self._actually_calculated_columns]

        if len(indices) == len(columnNameList):
            indices = None

        if not hasattr(self, '_ssmMagDict'):
            self._ssmMagDict = {}
            self._ssmMagNormDict = {}
            self._file_dir = getPackageDir('sims_sed_library')
            self._spec_map = defaultSpecMap
            self._normalizing_bandpass = Bandpass()
            self._normalizing_bandpass.imsimBandpass()

        sedNameList = self.column_by_name('sedFilename')
        magNormList = self.column_by_name('magNorm')

        if len(sedNameList)==0:
            # need to return something when InstanceCatalog goes through
            # it's "dry run" to determine what columns are required from
            # the database
            return np.zeros((len(bandpassDict.keys()),0))

        magListOut = []

        for sedName, magNorm in zip(sedNameList, magNormList):
            magTag = bandpassTag+'_'+sedName
            if sedName not in self._ssmMagNormDict or magTag not in self._ssmMagDict:
                dummySed = Sed()
                dummySed.readSED_flambda(os.path.join(self._file_dir, self._spec_map[sedName]))
                fnorm = dummySed.calcFluxNorm(magNorm, self._normalizing_bandpass)
                dummySed.multiplyFluxNorm(fnorm)
                magList = bandpassDict.magListForSed(dummySed, indices=indices)
                self._ssmMagDict[magTag] = magList
                self._ssmMagNormDict[sedName] = magNorm
            else:
                dmag = magNorm - self._ssmMagNormDict[sedName]
                magList = self._ssmMagDict[magTag] + dmag
            magListOut.append(magList)

        return np.array(magListOut).transpose()

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import multiplyFluxNorm [as 别名]
    def testObjectPlacement(self):
        """
        Test that GalSim places objects on the correct pixel by drawing
        images, reading them in, and then comparing the flux contained in
        circles of 2 fwhm radii about the object's expected positions with
        the actual expected flux of the objects.
        """
        scratchDir = os.path.join(getPackageDir('sims_GalSimInterface'), 'tests', 'scratchSpace')
        catName = os.path.join(scratchDir, 'placementCatalog.dat')
        imageRoot = os.path.join(scratchDir, 'placementImage')
        dbFileName = os.path.join(scratchDir, 'placementInputCatalog.dat')

        cameraDir = os.path.join(getPackageDir('sims_GalSimInterface'), 'tests', 'cameraData')
        camera = ReturnCamera(cameraDir)
        detector = camera[0]
        imageName = '%s_%s_u.fits' % (imageRoot, detector.getName())

        controlSed = Sed()
        controlSed.readSED_flambda(
                                   os.path.join(getPackageDir('sims_sed_library'),
                                               'flatSED','sed_flat.txt.gz')
                                   )

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

        controlBandpass = Bandpass()
        controlBandpass.imsimBandpass()

        ff = controlSed.calcFluxNorm(self.magNorm, uBandpass)
        controlSed.multiplyFluxNorm(ff)
        a_int, b_int = controlSed.setupCCMab()
        controlSed.addCCMDust(a_int, b_int, A_v=0.1, R_v=3.1)

        nSamples = 5
        numpy.random.seed(42)
        pointingRaList = numpy.random.random_sample(nSamples)*360.0
        pointingDecList = numpy.random.random_sample(nSamples)*180.0 - 90.0
        rotSkyPosList = numpy.random.random_sample(nSamples)*360.0
        fwhmList = numpy.random.random_sample(nSamples)*1.0 + 0.3

        actualCounts = None

        for pointingRA, pointingDec, rotSkyPos, fwhm in \
        zip(pointingRaList, pointingDecList, rotSkyPosList, fwhmList):


            obs = ObservationMetaData(unrefractedRA=pointingRA,
                                      unrefractedDec=pointingDec,
                                      boundType='circle',
                                      boundLength=4.0,
                                      mjd=49250.0,
                                      rotSkyPos=rotSkyPos)

            xDisplacementList = numpy.random.random_sample(nSamples)*60.0-30.0
            yDisplacementList = numpy.random.random_sample(nSamples)*60.0-30.0
            create_text_catalog(obs, dbFileName, xDisplacementList, yDisplacementList,
                                mag_norm=[self.magNorm]*len(xDisplacementList))
            db = placementFileDBObj(dbFileName, runtable='test')
            cat = placementCatalog(db, obs_metadata=obs)
            if actualCounts is None:
                actualCounts = controlSed.calcADU(uBandpass, cat.photParams)

            psf = SNRdocumentPSF(fwhm=fwhm)
            cat.setPSF(psf)
            cat.camera = camera

            cat.write_catalog(catName)
            cat.write_images(nameRoot=imageRoot)

            objRaList = []
            objDecList = []
            with open(catName, 'r') as inFile:
                for line in inFile:
                    if line[0] != '#':
                        words = line.split(';')
                        objRaList.append(numpy.radians(numpy.float(words[2])))
                        objDecList.append(numpy.radians(numpy.float(words[3])))

            objRaList = numpy.array(objRaList)
            objDecList = numpy.array(objDecList)

            self.check_placement(imageName, objRaList, objDecList,
                                [fwhm]*len(objRaList),
                                numpy.array([actualCounts]*len(objRaList)),
                                cat.photParams.gain, detector, camera, obs, epoch=2000.0)

            if os.path.exists(dbFileName):
                os.unlink(dbFileName)
            if os.path.exists(catName):
                os.unlink(catName)
            if os.path.exists(imageName):
                os.unlink(imageName)