Python Sed.calcADU方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcADU [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 calcADU [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 calcADU [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 calcADU [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 calcADU [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 calcADU [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 calcADU [as 别名]
# scaled properly, in other bandpasses).

mag_desired = 24.5
print "Now going to apply a scaling factor to the SED to set magnitude to %.4f" %(mag_desired)

# Calculate the scaling factor.
fluxnorm = star.calcFluxNorm(mag_desired, rband)
# Apply the scaling factor. 
star.multiplyFluxNorm(fluxnorm)

# Try the magnitude calculation again. 
mag = star.calcMag(rband)
print "After scaling, magnitude of SED is now %.4f (desired magnitude was %.4f)" %(mag, mag_desired)

# And let's calculate what the expected photon counts for LSST would be.
counts = star.calcADU(rband, expTime=30)
print "This would correspond to roughly %f counts in the LSST focal plane, in a 30s exposure." %(counts)

# For fun, let's see what else can happen.

ebv = 0.5
print ""
print "Let's try adding %.2f E(B-V) dust extinction to this star." %(ebv)
a, b = star.setupCCMab()
# You can use addCCMDust on the 'star' object itself, but I'm illustrating here how you could also
#  do otherwise - preserve the original 'star' object as is, and create a new Sed object that does
#  include the effects of dust ('dustystar').  Star's data will be unchanged by the dust. 
dustywavelen, dustyflambda = star.addCCMDust(a, b, ebv=ebv, wavelen=star.wavelen, flambda=star.flambda)
dustystar = Sed(wavelen=dustywavelen, flambda=dustyflambda)
magdust = dustystar.calcMag(rband)
print "With this dust, the magnitude of the star in this bandpass is now %.4f." %(magdust)

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

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import calcADU [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 calcADU [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 calcADU [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 calcADU [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')
#.........这里部分代码省略.........