Python Sed.setFlatSED方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
    def __init__(self, m5Col='fiveSigmaDepth', units='mag', maps=['DustMap'],
                 lsstFilter='r', wavelen_min=None , wavelen_max=None , wavelen_step=1., **kwargs ):
        """
        Args:
            m5Col (str): Column name that ('fiveSigmaDepth')
            units (str): units of the metric ('mag')
            maps (list): List of maps to use with the metric (['DustMap'])
            lsstFilter (str): Which LSST filter to calculate m5 for
            wavelen_min (float): Minimum wavength of your filter (None)
            wavelen_max (float): (None)
            wavelen_step (float): (1.)
            **kwargs:
        """

        waveMins={'u':330.,'g':403.,'r':552.,'i':691.,'z':818.,'y':950.}
        waveMaxes={'u':403.,'g':552.,'r':691.,'i':818.,'z':922.,'y':1070.}

        if lsstFilter is not None:
            wavelen_min = waveMins[lsstFilter]
            wavelen_max = waveMaxes[lsstFilter]

        self.m5Col = m5Col
        super(ExgalM5, self).__init__(col=[self.m5Col],
                                      maps=maps, units=units, **kwargs)

        testsed = Sed()
        testsed.setFlatSED(wavelen_min = wavelen_min,
                           wavelen_max = wavelen_max, wavelen_step = 1)
        self.a,self.b = testsed.setupCCMab()
        self.R_v = 3.1
        self.Coaddm5Metric = Coaddm5Metric(m5Col=m5Col)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [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 setFlatSED [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 setFlatSED [as 别名]
    def __init__(self, m5Col='fiveSigmaDepth', units='mag', maps=['DustMap'],
                 lsstFilter='r', wavelen_min=None , wavelen_max=None , wavelen_step=1., **kwargs ):
        """

        """

        waveMins={'u':330.,'g':403.,'r':552.,'i':691.,'z':818.,'y':950.}
        waveMaxes={'u':403.,'g':552.,'r':691.,'i':818.,'z':922.,'y':1070.}

        if lsstFilter is not None:
            wavelen_min = waveMins[lsstFilter]
            wavelen_max = waveMaxes[lsstFilter]

        self.m5Col = m5Col
        super(ExgalM5, self).__init__(col=[self.m5Col],
                                      maps=maps, units=units, **kwargs)

        testsed = Sed()
        testsed.setFlatSED(wavelen_min = wavelen_min,
                           wavelen_max = wavelen_max, wavelen_step = 1)
        self.a,self.b = testsed.setupCCMab()
        self.R_v = 3.1
        self.Coaddm5Metric = Coaddm5Metric(m5Col=m5Col)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
    def testApplication(self):
        """
        Test that PhotometricParameters get properly propagated into
        Sed methods.  We will test this using Sed.calcADU, since the ADU
        scale linearly with the appropriate parameter.
        """

        testSed = Sed()
        testSed.setFlatSED()

        testBandpass = Bandpass()
        testBandpass.readThroughput(os.path.join(lsst.utils.getPackageDir('throughputs'),
                                                 'baseline', 'total_g.dat'))

        control = testSed.calcADU(testBandpass,
                                  photParams=PhotometricParameters())

        testCase = PhotometricParameters(exptime=30.0)

        test = testSed.calcADU(testBandpass, photParams=testCase)

        self.assertGreater(control, 0.0)
        self.assertEqual(control, 0.5*test)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [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 setFlatSED [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 setFlatSED [as 别名]
def calcM5(hardware, system, atmos, title='m5'):
    """
    Calculate m5 values for all filters in hardware and system.
    Prints all values that go into "table 2" of the overview paper.
    Returns dictionary of m5 values.
    """
    # photParams stores default values for the exposure time, nexp, size of the primary,
    #  readnoise, gain, platescale, etc.
    # See https://github.com/lsst/sims_photUtils/blob/master/python/lsst/sims/photUtils/PhotometricParameters.py
    photParams = PhotometricParameters(gain=1)
    photParams_infinity = PhotometricParameters(readnoise=0, darkcurrent=0,
                                                othernoise=0, gain=1)
    # lsstDefaults stores default values for the FWHMeff.
    # See https://github.com/lsst/sims_photUtils/blob/master/python/lsst/sims/photUtils/LSSTdefaults.py
    lsstDefaults = LSSTdefaults()
    darksky = Sed()
    darksky.readSED_flambda(os.path.join('../siteProperties', 'darksky.dat'))
    flatSed = Sed()
    flatSed.setFlatSED()
    m5 = {}
    Tb = {}
    Sb = {}
    kAtm = {}
    Cm = {}
    dCm_infinity = {}
    sourceCounts = {}
    skyCounts = {}
    skyMag = {}
    gamma = {}
    for f in system:
        m5[f] = SignalToNoise.calcM5(darksky, system[f], hardware[f], photParams,
                                     FWHMeff=lsstDefaults.FWHMeff(f))
        fNorm = flatSed.calcFluxNorm(m5[f], system[f])
        flatSed.multiplyFluxNorm(fNorm)
        sourceCounts[f] = flatSed.calcADU(system[f], photParams=photParams)
        # Calculate the Skycounts expected in this bandpass.
        skyCounts[f] = (darksky.calcADU(hardware[f], photParams=photParams)
                        * photParams.platescale**2)
        # Calculate the sky surface brightness.
        skyMag[f] = darksky.calcMag(hardware[f])
        # Calculate the gamma value.
        gamma[f] = SignalToNoise.calcGamma(system[f], m5[f], photParams)
        # Calculate the "Throughput Integral" (this is the hardware + atmosphere)
        dwavelen = np.mean(np.diff(system[f].wavelen))
        Tb[f] = np.sum(system[f].sb / system[f].wavelen) * dwavelen
        # Calculate the "Sigma" 'system integral' (this is the hardware only)
        Sb[f] = np.sum(hardware[f].sb / hardware[f].wavelen) * dwavelen
        # Calculate km - atmospheric extinction in a particular bandpass
        kAtm[f] = -2.5*np.log10(Tb[f] / Sb[f])
        # Calculate the Cm and Cm_Infinity values.
        # m5 = Cm + 0.5*(msky - 21) + 2.5log10(0.7/FWHMeff) + 1.25log10(t/30) - km(X-1.0)
        # Exptime should be 30 seconds and X=1.0
        exptime = photParams.exptime * photParams.nexp
        if exptime != 30.0:
            print "Whoa, exposure time was not as expected - got %s not 30 seconds. Please edit Cm calculation." %(exptime)
        # Assumes atmosphere used in system throughput is X=1.0
        X = 1.0
        Cm[f] = (m5[f] - 0.5*(skyMag[f] - 21) - 2.5*np.log10(0.7/lsstDefaults.FWHMeff(f)))
        # Calculate Cm_Infinity by setting readout noise to zero.
        m5inf = SignalToNoise.calcM5(darksky, system[f], hardware[f],  photParams_infinity,
                                     FWHMeff=lsstDefaults.FWHMeff(f))
        Cm_infinity = (m5inf - 0.5*(skyMag[f] - 21)
                       - 2.5*np.log10(0.7/lsstDefaults.FWHMeff(f)))
        dCm_infinity[f] = Cm_infinity - Cm[f]
    print title
    print 'Filter FWHMeff FWHMgeom SkyMag SkyCounts Tb Sb kAtm Gamma Cm dCm_infinity m5 SourceCounts'
    for f in ('u', 'g' ,'r', 'i', 'z', 'y'):
        print '%s %.2f %.2f %.2f %.1f %.3f %.3f %.4f %.6f %.2f %.2f %.2f %.2f'\
           %(f, lsstDefaults.FWHMeff(f),
             SignalToNoise.FWHMeff2FWHMgeom(lsstDefaults.FWHMeff(f)),
             skyMag[f], skyCounts[f], Tb[f], Sb[f], kAtm[f],
             gamma[f], Cm[f], dCm_infinity[f], m5[f], sourceCounts[f])

    # Show what these look like individually (add sky & m5 limits on throughput curves)
    plt.figure()
    for f in filterlist:
        plt.plot(system[f].wavelen, system[f].sb, color=filtercolors[f], linewidth=2, label=f)
    plt.plot(atmosphere.wavelen, atmosphere.sb, 'k:', label='X=1.0')
    plt.legend(loc='center right', fontsize='smaller')
    plt.xlim(300, 1100)
    plt.ylim(0, 1)
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Throughput')
    plt.title('System Throughputs')
    plt.grid(True)

    plt.figure()
    ax = plt.gca()
    # Add dark sky
    ax2 = ax.twinx()
    plt.sca(ax2)
    skyab = -2.5*np.log10(darksky.fnu) - darksky.zp
    ax2.plot(darksky.wavelen, skyab,
             'k-', linewidth=0.8, label='Dark sky mags')
    ax2.set_ylabel('AB mags')
    ax2.set_ylim(24, 14)
    plt.sca(ax)
    # end of dark sky
    handles = []
    for f in filterlist:
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
def calcM5s(hardware, system, atmos, title='m5'):
    photParams = PhotometricParameters()
    lsstDefaults = LSSTdefaults()
    darksky = Sed()
    darksky.readSED_flambda(os.path.join(os.getenv('SYSENG_THROUGHPUTS_DIR'), 'siteProperties', 'darksky.dat'))
    flatSed = Sed()
    flatSed.setFlatSED()
    m5 = {}
    sourceCounts = {}
    skyCounts = {}
    skyMag = {}
    gamma = {}
    for f in system:
        m5[f] = SignalToNoise.calcM5(darksky, system[f], hardware[f], photParams, seeing=lsstDefaults.seeing(f))
        fNorm = flatSed.calcFluxNorm(m5[f], system[f])
        flatSed.multiplyFluxNorm(fNorm)
        sourceCounts[f] = flatSed.calcADU(system[f], photParams=photParams)
        # Calculate the Skycounts expected in this bandpass.
        skyCounts[f] = darksky.calcADU(hardware[f], photParams=photParams) * photParams.platescale**2
        # Calculate the sky surface brightness.
        skyMag[f] = darksky.calcMag(hardware[f])
        # Calculate the gamma value.
        gamma[f] = SignalToNoise.calcGamma(system[f], m5[f], photParams)
    print title
    print 'Filter m5 SourceCounts SkyCounts SkyMag Gamma'
    for f in ('u', 'g' ,'r', 'i', 'z', 'y'):
        print '%s %.2f %.1f %.2f %.2f %.6f' %(f, m5[f], sourceCounts[f], skyCounts[f], skyMag[f], gamma[f])

    # Show what these look like individually (add sky & m5 limits on throughput curves)
    plt.figure()
    ax = plt.gca()
    # Add dark sky
    ax2 = ax.twinx()
    plt.sca(ax2)
    skyab = -2.5*np.log10(darksky.fnu) - darksky.zp
    ax2.plot(darksky.wavelen, skyab,
             'k-', linewidth=0.8, label='Dark sky mags')
    ax2.set_ylabel('AB mags')
    ax2.set_ylim(24, 10)
    plt.sca(ax)
    # end of dark sky
    handles = []
    for f in filterlist:
        plt.plot(system[f].wavelen, system[f].sb, color=filtercolors[f], linewidth=2)
        myline = mlines.Line2D([], [], color=filtercolors[f], linestyle='-', linewidth=2,
                               label = '%s: m5 %.1f (sky %.1f)' %(f, m5[f], skyMag[f]))
        handles.append(myline)
    plt.plot(atmos.wavelen, atmos.sb, 'k:', label='Atmosphere, X=1.2')
    # Add legend for dark sky.
    myline = mlines.Line2D([], [], color='k', linestyle='-', label='Dark sky AB mags')
    handles.append(myline)
    # end of dark sky legend line
    plt.legend(loc=(0.01, 0.69), handles=handles, fancybox=True, numpoints=1, fontsize='small')
    plt.ylim(0, 1)
    plt.xlim(300, 1100)
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Fractional Throughput Response')
    if title == 'Vendor combo':
        title = ''
    plt.title('System total response curves %s' %(title))
    if title == '':
        plt.savefig('throughputs.pdf', format='pdf', dpi=600)
    return m5

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
def calcM5(hardware, system, atmos, title="m5"):
    effarea = np.pi * (6.423 / 2.0 * 100.0) ** 2
    photParams = PhotometricParameters(effarea=effarea)
    lsstDefaults = LSSTdefaults()
    darksky = Sed()
    darksky.readSED_flambda(os.path.join("../siteProperties", "darksky.dat"))
    flatSed = Sed()
    flatSed.setFlatSED()
    m5 = {}
    sourceCounts = {}
    skyCounts = {}
    skyMag = {}
    gamma = {}
    for f in system:
        m5[f] = SignalToNoise.calcM5(darksky, system[f], hardware[f], photParams, FWHMeff=lsstDefaults.FWHMeff(f))
        fNorm = flatSed.calcFluxNorm(m5[f], system[f])
        flatSed.multiplyFluxNorm(fNorm)
        sourceCounts[f] = flatSed.calcADU(system[f], photParams=photParams)
        # Calculate the Skycounts expected in this bandpass.
        skyCounts[f] = darksky.calcADU(hardware[f], photParams=photParams) * photParams.platescale ** 2
        # Calculate the sky surface brightness.
        skyMag[f] = darksky.calcMag(hardware[f])
        # Calculate the gamma value.
        gamma[f] = SignalToNoise.calcGamma(system[f], m5[f], photParams)
    print title
    print "Filter m5 SourceCounts SkyCounts SkyMag Gamma"
    for f in ("u", "g", "r", "i", "z", "y"):
        print "%s %.2f %.1f %.2f %.2f %.6f" % (f, m5[f], sourceCounts[f], skyCounts[f], skyMag[f], gamma[f])

    # Show what these look like individually (add sky & m5 limits on throughput curves)
    plt.figure()
    ax = plt.gca()
    # Add dark sky
    ax2 = ax.twinx()
    plt.sca(ax2)
    skyab = -2.5 * np.log10(darksky.fnu) - darksky.zp
    ax2.plot(darksky.wavelen, skyab, "k-", linewidth=0.8, label="Dark sky mags")
    ax2.set_ylabel("AB mags")
    ax2.set_ylim(24, 14)
    plt.sca(ax)
    # end of dark sky
    handles = []
    for f in filterlist:
        plt.plot(system[f].wavelen, system[f].sb, color=filtercolors[f], linewidth=2)
        myline = mlines.Line2D(
            [],
            [],
            color=filtercolors[f],
            linestyle="-",
            linewidth=2,
            label="%s: m5 %.1f (sky %.1f)" % (f, m5[f], skyMag[f]),
        )
        handles.append(myline)
    plt.plot(atmos.wavelen, atmos.sb, "k:", label="Atmosphere, X=1.0 with aerosols")
    # Add legend for dark sky.
    myline = mlines.Line2D([], [], color="k", linestyle="-", label="Dark sky AB mags")
    handles.append(myline)
    # end of dark sky legend line
    plt.legend(loc=(0.01, 0.69), handles=handles, fancybox=True, numpoints=1, fontsize="small")
    plt.ylim(0, 1)
    plt.xlim(300, 1100)
    plt.xlabel("Wavelength (nm)")
    plt.ylabel("Fractional Throughput Response")
    if title == "Vendor combo":
        title = ""
    plt.title("System total response curves %s" % (title))
    plt.savefig("../plots/system+sky" + title + ".png", format="png", dpi=600)
    return m5

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
def calc_adu(mag, bandpass):
    sed = Sed()
    sed.setFlatSED()
    fluxNorm = sed.calcFluxNorm(mag, bandpass)
    sed.multiplyFluxNorm(fluxNorm)
    return sed.calcADU(bandpass, fake_phot_params())

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
def calcM5(hardware, system, atmos, title='m5', return_t2_values=False):
    """
    Calculate m5 values for all filters in hardware and system.
    Prints all values that go into "table 2" of the overview paper.
    Returns dictionary of m5 values.
    """
    # photParams stores default values for the exposure time, nexp, size of the primary,
    #  readnoise, gain, platescale, etc.
    # See https://github.com/lsst/sims_photUtils/blob/master/python/lsst/sims/photUtils/PhotometricParameters.py
    effarea = np.pi * (6.423/2.*100.)**2
    photParams_zp = PhotometricParameters(exptime=1, nexp=1, gain=1, effarea=effarea,
                                          readnoise=8.8, othernoise=0, darkcurrent=0.2)
    photParams = PhotometricParameters(gain=1.0, effarea=effarea, readnoise=8.8, othernoise=0, darkcurrent=0.2)
    photParams_infinity = PhotometricParameters(gain=1.0, readnoise=0, darkcurrent=0,
                                                othernoise=0, effarea=effarea)
    # lsstDefaults stores default values for the FWHMeff.
    # See https://github.com/lsst/sims_photUtils/blob/master/python/lsst/sims/photUtils/LSSTdefaults.py
    lsstDefaults = LSSTdefaults()
    darksky = Sed()
    darksky.readSED_flambda(os.path.join(getPackageDir('syseng_throughputs'),
                                         'siteProperties', 'darksky.dat'))
    flatSed = Sed()
    flatSed.setFlatSED()
    m5 = {}
    Tb = {}
    Sb = {}
    kAtm = {}
    Cm = {}
    dCm_infinity = {}
    sourceCounts = {}
    skyCounts = {}
    skyMag = {}
    gamma = {}
    zpT = {}
    FWHMgeom = {}
    FWHMeff = {}
    for f in system:
        zpT[f] = system[f].calcZP_t(photParams_zp)
        m5[f] = SignalToNoise.calcM5(darksky, system[f], hardware[f], photParams, FWHMeff=lsstDefaults.FWHMeff(f))
        fNorm = flatSed.calcFluxNorm(m5[f], system[f])
        flatSed.multiplyFluxNorm(fNorm)
        sourceCounts[f] = flatSed.calcADU(system[f], photParams=photParams)
        # Calculate the Skycounts expected in this bandpass.
        skyCounts[f] = (darksky.calcADU(hardware[f], photParams=photParams)
                        * photParams.platescale**2)
        # Calculate the sky surface brightness.
        skyMag[f] = darksky.calcMag(hardware[f])
        # Calculate the gamma value.
        gamma[f] = SignalToNoise.calcGamma(system[f], m5[f], photParams)
        # Calculate the "Throughput Integral" (this is the hardware + atmosphere)
        dwavelen = np.mean(np.diff(system[f].wavelen))
        Tb[f] = np.sum(system[f].sb / system[f].wavelen) * dwavelen
        # Calculate the "Sigma" 'system integral' (this is the hardware only)
        Sb[f] = np.sum(hardware[f].sb / hardware[f].wavelen) * dwavelen
        # Calculate km - atmospheric extinction in a particular bandpass
        kAtm[f] = -2.5*np.log10(Tb[f] / Sb[f])
        # Calculate the Cm and Cm_Infinity values.
        # m5 = Cm + 0.5*(msky - 21) + 2.5log10(0.7/FWHMeff) + 1.25log10(t/30) - km(X-1.0)
        # Assumes atmosphere used in system throughput is X=1.0
        X = 1.0
        Cm[f] = (m5[f] - 0.5*(skyMag[f] - 21) - 2.5*np.log10(0.7/lsstDefaults.FWHMeff(f))
                 - 1.25*np.log10((photParams.exptime*photParams.nexp)/30.0) + kAtm[f]*(X-1.0))
        # Calculate Cm_Infinity by setting readout noise to zero.
        m5inf = SignalToNoise.calcM5(darksky, system[f], hardware[f],  photParams_infinity,
                                     FWHMeff=lsstDefaults.FWHMeff(f))
        Cm_infinity = (m5inf - 0.5*(skyMag[f] - 21) - 2.5*np.log10(0.7/lsstDefaults.FWHMeff(f))
                       - 1.25*np.log10((photParams.exptime*photParams.nexp)/30.0) + kAtm[f]*(X-1.0))
        dCm_infinity[f] = Cm_infinity - Cm[f]
    print 'Filter FWHMeff FWHMgeom SkyMag SkyCounts Zp_t Tb Sb kAtm Gamma Cm dCm_infinity m5 SourceCounts'
    for f in ('u', 'g' ,'r', 'i', 'z', 'y'):
        FWHMeff[f] = lsstDefaults.FWHMeff(f)
        FWHMgeom[f] = SignalToNoise.FWHMeff2FWHMgeom(lsstDefaults.FWHMeff(f))
        print '%s %.2f %.2f %.2f %.1f %.2f %.3f %.3f %.4f %.6f %.2f %.2f %.2f %.2f'\
           % (f, FWHMeff[f], FWHMgeom[f],
              skyMag[f], skyCounts[f], zpT[f], Tb[f], Sb[f], kAtm[f],
              gamma[f], Cm[f], dCm_infinity[f], m5[f], sourceCounts[f])
    if return_t2_values:
        return {'FHWMeff': FWHMeff, 'FWHMgeom': FWHMgeom, 'skyMag': skyMag, 'skycounts': skyCounts,
                'zpT': zpT, 'Tb': Tb, 'Sb': Sb, 'kAtm': kAtm,
                'gamma': gamma, 'Cm': Cm, 'dCm_infinity': dCm_infinity,
                'm5': m5, 'sourceCounts': sourceCounts}

    for f in filterlist:
        m5_cm = Cm[f] + 0.5*(skyMag[f] - 21.0) + 2.5*np.log10(0.7/lsstDefaults.FWHMeff(f))
        if m5_cm - m5[f] > 0.001:
            raise ValueError('Cm calculation for %s band is incorrect! m5_cm != m5_snr' %f)

    # Show what these look like individually (add sky & m5 limits on throughput curves)
    plt.figure()
    for f in filterlist:
        plt.plot(system[f].wavelen, system[f].sb, color=filtercolors[f], linewidth=2, label=f)
    plt.plot(atmosphere.wavelen, atmosphere.sb, 'k:', label='X=1.0')
    plt.legend(loc='center right', fontsize='smaller')
    plt.xlim(300, 1100)
    plt.ylim(0, 1)
    plt.xlabel('Wavelength (nm)')
    plt.ylabel('Throughput')
    plt.title('System Throughputs')
    plt.grid(True)
    plt.savefig('../plots/throughputs.png', format='png')
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import setFlatSED [as 别名]
                          wavelen_max=wavelen_max,
                          wavelen_step = wavelen_step)
    for f in filterlist:
        mags1[f][i] = tmpgal.calcMag(lsstbp[f])
dt, t = dtime(t)
print "Calculating dust/redshift/dust/fluxnorm/%d magnitudes for %d galaxies took %f s" \
      %(len(filterlist), num_gal, dt)

# For next test: want to also do all the same steps, but in an optimized form. This means
# doing some things that Sed does 'behind the scenes' explicitly, but also means the code may be a little
# harder to read at first.
# First: calculate internal a/b on wavelength range required for internal dust extinction.
a_int, b_int = gals[gallist[0]].setupCCMab()  # this is a/b on native galaxy sed range. 
# Next: calculate milky way a/b on wavelength range required for calculating magnitudes - i.e. 300 to 1200 nm.
tmpgal = Sed()
tmpgal.setFlatSED(wavelen_min=wavelen_min, wavelen_max=wavelen_max, wavelen_step = wavelen_step)
a_mw, b_mw = tmpgal.setupCCMab()  # so this is a/b on native MW bandpass range. 
# Also: set up phi for each bandpass - ahead of time. And set up a list of bandpasses, then create phiarray 
# and dlambda to set up for manyMagCalc method.
bplist = []
for f in filterlist:
    lsstbp[f].sbTophi()
    bplist.append(lsstbp[f])
phiarray, dlambda = tmpgal.setupPhiArray(bplist)
# Set up dictionary + arrays to hold calculated magnitude information. 
mags2 = {}
for f in filterlist:
    mags2[f] = numpy.zeros(num_gal, dtype='float')
# For each galaxy (in num_gal's), apply internal dust, redshift, resample to 300-1200 nm, apply MW dust on
#   shorter (and standardized) wavelength range, fluxnorm, and then calculate mags using manyMagCalc. 
for i in range(num_gal):