Python Bandpass.setBandpass方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def buildVendorDetector(vendorDir, addLosses=True):
    """
    Builds a detector response from the files in vendorDir, by reading the *_QE.dat
      and *_Losses subdirectory for a single version of the detector.
    Returns a Bandpass object.
    If addLosses is True, the QE curve is multiplied by the losses in the *Losses.dat files.
    If addLosses is False, the QE curve does not have any losses included.
    """
    # Read the QE file.
    qefile = glob(os.path.join(vendorDir, '*_QE.dat'))
    if len(qefile) != 1:
        raise ValueError('Expected a single QE file in this directory, found: ', qefile)
    qefile = qefile[0]
    qe = Bandpass()
    qe.readThroughput(qefile)
    if addLosses:
        loss = _readLosses(vendorDir)
        wavelength, sb = qe.multiplyThroughputs(loss.wavelen, loss.sb)
        qe.setBandpass(wavelength, sb)
    # Verify that no values go significantly below zero.
    belowzero = np.where(qe.sb < 0)
    # If there are QE values significantly < 0, raise an exception.
    if qe.sb[belowzero] < belowZeroThreshhold:
        raise ValueError('Found values in QE response significantly below zero.')
    # If they are just small errors in interpolation, set to zero.
    qe.sb[belowzero] = 0
    return qe

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
    def testMags(self):
        """
        Test that the interpolated mags are similar to mags computed from interpolated spectra
        """

        throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
        filters = ['u', 'g', 'r', 'i', 'z', 'y']

        bps = []
        for filterName in filters:
            bp = np.loadtxt(os.path.join(throughPath, 'filter_%s.dat' % filterName),
                            dtype=zip(['wave', 'trans'], [float]*2))
            lsst_bp = Bandpass()
            lsst_bp.setBandpass(bp['wave'], bp['trans'])
            bps.append(lsst_bp)

        sm1 = self.sm_spec
        sm1.setRaDecMjd([36.], [-68.], 49353.18, degrees=True)
        mags1 = []
        for bp in bps:
            mags1.append(sm1.returnMags(bandpass=bp))
        mags1 = np.array(mags1)

        sm2 = self.sm_mags
        sm2.setRaDecMjd([36.], [-68.], 49353.18, degrees=True)
        mag2 = sm2.returnMags()
        for i, filtername in enumerate(filters):
            np.testing.assert_allclose(mags1[i, :], mag2[filtername], rtol=1e-4)

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def buildMirror(mirrorDir, addLosses=True):
    """
    Build a mirror throughput curve.
    Assumes there are *Losses.dat subdirectory with loss files
       and a m*_Ideal.dat file with the mirror throughput.
    Returns a bandpass object.
    If addLosses is True, the *_Ideal.dat file is multiplied by the *_Losses/*.dat files.
    """
    # Read the mirror reflectance curve.
    mirrorfile = glob(os.path.join(mirrorDir, 'm*Ideal.dat'))
    if len(mirrorfile) != 1:
        raise ValueError('Expected a single mirror file in directory %s, found: ' %mirrorDir, mirrorfile)
    mirrorfile = mirrorfile[0]
    mirror = Bandpass()
    mirror.readThroughput(mirrorfile)
    if addLosses:
        loss = _readLosses(mirrorDir)
        wavelen, sb = mirror.multiplyThroughputs(loss.wavelen, loss.sb)
        mirror.setBandpass(wavelen, sb)
    # Verify that no values go significantly below zero.
    belowzero = np.where(mirror.sb < 0)
    # If there are QE values significantly < 0, raise an exception.
    if mirror.sb[belowzero] < belowZeroThreshhold:
        raise ValueError('Found values in mirror response significantly below zero')
    # If they are just small errors in interpolation, set to zero.
    mirror.sb[belowzero] = 0
    return mirror

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
    def testMags(self):
        """
        Test that the interpolated mags are similar to mags computed from interpolated spectra
        """

        throughPath = os.path.join(getPackageDir("throughputs"), "baseline")
        filters = ["u", "g", "r", "i", "z", "y"]

        bps = []
        for filterName in filters:
            bp = np.loadtxt(
                os.path.join(throughPath, "filter_%s.dat" % filterName), dtype=zip(["wave", "trans"], [float] * 2)
            )
            lsst_bp = Bandpass()
            lsst_bp.setBandpass(bp["wave"], bp["trans"])
            bps.append(lsst_bp)

        sm1 = sb.SkyModel()
        sm1.setRaDecMjd([36.0], [-68.0], 49353.18, degrees=True)
        sm1.computeSpec()
        mags1 = []
        for bp in bps:
            mags1.append(sm1.computeMags(bandpass=bp))
        mags1 = np.array(mags1)

        sm2 = sb.SkyModel(mags=True)
        sm2.setRaDecMjd([36.0], [-68.0], 49353.18, degrees=True)
        sm2.computeSpec()
        mag2 = sm2.computeMags()
        np.testing.assert_allclose(mags1, mag2.T, rtol=1e-4)

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def lsst_filters():
    throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
    lsstKeys = ['u', 'g', 'r', 'i', 'z', 'y']
    bps = {}
    for key in lsstKeys:
        bp = np.loadtxt(os.path.join(throughPath, 'total_'+key+'.dat'),
                        dtype=zip(['wave', 'trans'], [float]*2))
        bpTemp = Bandpass()
        good = np.where(bp['trans'] > 0.)
        bpTemp.setBandpass(bp['wave'], bp['trans'], wavelen_min=bp['wave'][good].min(),
                           wavelen_max=bp['wave'][good].max())
        bps[key] = bpTemp
    return bps

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def addAerosol(atmosphere, X, tau0=0.05, alpha=1.0, wavelen0=440.0, plotAtmosphere=True):
    # Calculate the aerosol contribution -- sb with aerosols = sb*exp(-tau)
    tau = tau0 * np.power((wavelen0/atmosphere.wavelen), alpha)
    # Generate new atmosphere bandpass with aerosols.
    atmosphere_aerosol = Bandpass()
    atmosphere_aerosol.setBandpass(wavelen = atmosphere.wavelen,
                                   sb = atmosphere.sb * np.exp(-tau*X))

    if plotAtmosphere:
        # Plot for a check:
        atmodict = {'Original atmosphere':atmosphere,
                    'With aerosols': atmosphere_aerosol}
        bu.plotBandpasses(atmodict)

    return atmosphere_aerosol

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def stubb_fitlers(wave_min=350., wave_max=1050):
    """
    Define some narrow filters that overlap LSST u and y, and are in GAIA overlap.
    """
    throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
    bps = {}
    lsstKeys = ['u', 'y']
    bps = {}
    for key in lsstKeys:
        bp = np.loadtxt(os.path.join(throughPath, 'total_'+key+'.dat'),
                        dtype=zip(['wave', 'trans'], [float]*2))
        bpTemp = Bandpass()
        good = np.where((bp['trans'] > 0.) & (bp['wave'] > wave_min) & (bp['wave'] < wave_max))
        bpTemp.setBandpass(bp['wave'], bp['trans'], wavelen_min=bp['wave'][good].min(),
                           wavelen_max=bp['wave'][good].max())
        bps[key+'_truncated'] = bpTemp
    return bps

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def calcWDColors():
    """
    Calculate a few example WD colors. Values to go in stellarMags(). Here in case
    values need to be regenerated (different stars, bandpasses change, etc.)
    """

    try:
        from lsst.utils import getPackageDir
        import os
        from lsst.sims.photUtils import Bandpass, Sed
    except:
        'Need to setup sims_photUtils to generate WD magnitudes.'

    names = ['HeWD_25200_80', 'WD_11000_85', 'WD_3000_85']
    fns = ['bergeron_He_24000_80.dat_25200.gz',
           'bergeron_10500_85.dat_11000.gz', 'bergeron_2750_85.dat_3000.gz']
    wdDir = os.path.join(getPackageDir('sims_sed_library'), 'starSED/wDs/')
    files = [os.path.join(wdDir, filename) for filename in fns]

    # Read in the LSST bandpasses
    bpNames = ['u', 'g', 'r', 'i', 'z', 'y']
    bps = []
    throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
    for key in bpNames:
        bp = np.loadtxt(os.path.join(throughPath, 'filter_' + key + '.dat'),
                        dtype=list(zip(['wave', 'trans'], [float] * 2)))
        tempB = Bandpass()
        tempB.setBandpass(bp['wave'], bp['trans'])
        bps.append(tempB)

    # Read in the SEDs and compute mags
    mags = []
    for filename in files:
        star = Sed()
        star.readSED_flambda(filename)
        singleMags = [star.calcMag(band) for band in bps]
        mags.append([singleMags[i - 1] - singleMags[i] for i in range(1, 6)])

    for maglist, fn, name in zip(mags, fns, names):
        format = (name, fn) + tuple(maglist)
        print("['%s', '%s', %f, %f, %f, %f, %f]" % format)

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def buildLens(lensDir, addLosses=True):
    """
    Build the lens throughput curve from the files in lensDir.
    Returns a bandpass object.
    The coatings for the lens are in *_Coatings, the loss files are in *_Losses.
    The borosilicate glass throughput is in l*_Glass.dat; this file is smoothed using the
     savitzsky_golay function.
    The glass response is multiplied by the coatings and (if addLosses is True),
      also the loss curves.
    """
    lens = Bandpass()
    # Read the glass base file.
    glassfile = glob(os.path.join(lensDir, 'l*_Glass.dat'))
    if len(glassfile) != 1:
        raise ValueError('Expected a single glass file in this directory, found: ', glassfile)
    glassfile = glassfile[0]
    glass = Bandpass()
    glass.readThroughput(glassfile)
    # Smooth the glass response.
    smoothSb = savitzky_golay(glass.sb, 31, 3)
    lens = Bandpass()
    lens.setBandpass(glass.wavelen, smoothSb)
    # Read the broad band antireflective (BBAR) coatings files.
    bbars = _readCoatings(lensDir)
    # Multiply the bbars by the glass.
    wavelen, sb = lens.multiplyThroughputs(bbars.wavelen, bbars.sb)
    lens.setBandpass(wavelen, sb)
    # Add losses.
    if addLosses:
        loss = _readLosses(lensDir)
        wavelen, sb = lens.multiplyThroughputs(loss.wavelen, loss.sb)
        lens.setBandpass(wavelen, sb)
    # Verify that no values go significantly below zero.
    belowzero = np.where(lens.sb < 0)
    # If there are QE values significantly < 0, raise an exception.
    if lens.sb[belowzero] < belowZeroThreshhold:
        raise ValueError('Found values in lens throughput significantly below zero.')
    # If they are just small errors in interpolation, set to zero.
    lens.sb[belowzero] = 0
    return lens

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def packageLowerAtm():

    dataDir = getPackageDir('SIMS_SKYBRIGHTNESS_DATA')
    outDir = os.path.join(dataDir, 'ESO_Spectra/LowerAtm')

    # Read in all the spectra from ESO call and package into a single npz file

    files = glob.glob('LowerAtm/skytable*.fits')

    temp = pyfits.open(files[0])
    wave = temp[1].data['lam'].copy()*1e3

    airmasses = []
    nightTimes = []
    specs = []

    for i,filename in enumerate(files):
        fits = pyfits.open(filename)
        if np.max(fits[1].data['flux']) > 0:
            specs.append(fits[1].data['flux'].copy())
            header = fits[0].header['comment']
            for card in header:
                if 'SKYMODEL.TARGET.AIRMASS' in card:
                    airmasses.append(float(card.split('=')[-1]))
                elif 'SKYMODEL.TIME' in card:
                    nightTimes.append(float(card.split('=')[-1]))

    airmasses = np.array(airmasses)
    nigtTimes = np.array(nightTimes)

    nrec = airmasses.size
    nwave = wave.size

    dtype = [('airmass', 'float'),
             ('nightTimes', 'float'),
             ('spectra', 'float', (nwave)), ('mags', 'float', (6))]
    Spectra = np.zeros(nrec, dtype=dtype)
    Spectra['airmass'] = airmasses
    Spectra['nightTimes'] = nightTimes
    Spectra['spectra'] = specs

    hPlank = 6.626068e-27 # erg s
    cLight = 2.99792458e10 # cm/s

    # Convert spectra from ph/s/m2/micron/arcsec2 to erg/s/cm2/nm/arcsec2
    Spectra['spectra'] = Spectra['spectra']/(100.**2)*hPlank*cLight/(wave*1e-7)/1e3

    # Sort things since this might be helpful later
    Spectra.sort(order=['airmass','nightTimes'])

    # Load LSST filters
    throughPath = os.path.join(getPackageDir('throughputs'),'baseline')
    keys = ['u','g','r','i','z','y']
    nfilt = len(keys)
    filters = {}
    for filtername in keys:
        bp = np.loadtxt(os.path.join(throughPath, 'filter_'+filtername+'.dat'),
                        dtype=zip(['wave','trans'],[float]*2 ))
        tempB = Bandpass()
        tempB.setBandpass(bp['wave'],bp['trans'])
        filters[filtername] = tempB

    filterWave = np.array([filters[f].calcEffWavelen()[0] for f in keys ])

    for i,spectrum in enumerate(Spectra['spectra']):
        tempSed = Sed()
        tempSed.setSED(wave,flambda=spectrum)
        for j,filtName in enumerate(keys):
            try:
                Spectra['mags'][i][j] = tempSed.calcMag(filters[filtName])
            except:
                pass

    np.savez(os.path.join(outDir,'Spectra.npz'), wave = wave, spec=Spectra, filterWave=filterWave)

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
def packageZodiacal():

    dataDir = getPackageDir('SIMS_SKYBRIGHTNESS_DATA')
    outDir = os.path.join(dataDir, 'ESO_Spectra/Zodiacal')

    nside = 4

    # Read in all the spectra from ESO call and package into a single npz file

    files = glob.glob('Zodiacal/skytable*.fits')

    temp = pyfits.open(files[0])
    wave = temp[1].data['lam'].copy()*1e3

    airmasses = []
    eclLon = []
    eclLat = []
    specs = []

    for i,filename in enumerate(files):
        fits = pyfits.open(filename)
        if np.max(fits[1].data['flux']) > 0:
            specs.append(fits[1].data['flux'].copy())
            header = fits[0].header['comment']
            for card in header:
                if 'SKYMODEL.TARGET.AIRMASS' in card:
                    airmasses.append(float(card.split('=')[-1]))
                elif 'SKYMODEL.ECL.LON' in card:
                    eclLon.append(float(card.split('=')[-1]))
                elif 'SKYMODEL.ECL.LAT' in card:
                    eclLat.append(float(card.split('=')[-1]))


    airmasses = np.array(airmasses)
    eclLon = np.array(eclLon)
    eclLat = np.array(eclLat)

    wrapA = np.where(eclLon < 0.)
    eclLon[wrapA] = eclLon[wrapA]+360.

    uAM = np.unique(airmasses)
    nAM = uAM.size
    nwave = wave.size

    dtype = [('airmass', 'float'),
             ('hpid', 'int' ),
             ('spectra', 'float', (nwave)), ('mags', 'float', (6))]
    npix = hp.nside2npix(nside)
    Spectra = np.zeros(nAM*npix, dtype=dtype)
    for i,am in enumerate(uAM):
        Spectra['airmass'][i*npix:i*npix+npix] = am
        Spectra['hpid'][i*npix:i*npix+npix] = np.arange(npix)


    for am, lat, lon, spec in zip(airmasses,eclLat, eclLon, specs):
        hpid = hp.ang2pix(nside, np.radians(lat+90.), np.radians(lon) )
        good = np.where( (Spectra['airmass'] == am) & (Spectra['hpid'] == hpid))
        Spectra['spectra'][good] = spec.copy()

    hPlank = 6.626068e-27 # erg s
    cLight = 2.99792458e10 # cm/s

    # Convert spectra from ph/s/m2/micron/arcsec2 to erg/s/cm2/nm/arcsec2
    Spectra['spectra'] = Spectra['spectra']/(100.**2)*hPlank*cLight/(wave*1e-7)/1e3

    # Sort things since this might be helpful later
    Spectra.sort(order=['airmass', 'hpid'])

    # Load LSST filters
    throughPath = os.path.join(getPackageDir('throughputs'),'baseline')
    keys = ['u','g','r','i','z','y']
    nfilt = len(keys)
    filters = {}
    for filtername in keys:
        bp = np.loadtxt(os.path.join(throughPath, 'filter_'+filtername+'.dat'),
                        dtype=zip(['wave','trans'],[float]*2 ))
        tempB = Bandpass()
        tempB.setBandpass(bp['wave'],bp['trans'])
        filters[filtername] = tempB

    filterWave = np.array([filters[f].calcEffWavelen()[0] for f in keys ])

    for i,spectrum in enumerate(Spectra['spectra']):
        tempSed = Sed()
        tempSed.setSED(wave,flambda=spectrum)
        for j,filtName in enumerate(keys):
            try:
                Spectra['mags'][i][j] = tempSed.calcMag(filters[filtName])
            except:
                pass


    #span this over multiple files to store in github
    nbreak = 3
    nrec = np.size(Spectra)

    for i in np.arange(nbreak):
        np.savez(os.path.join(outDir,'zodiacalSpectra_'+str(i)+'.npz'), wave = wave,
                 spec=Spectra[i*nrec/nbreak:(i+1)*nrec/nbreak], filterWave=filterWave)

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

#.........这里部分代码省略.........
                moonSpec.append(fits[1].data['flux'].copy())
                header = fits[0].header['comment']
                for card in header:
                    if 'SKYMODEL.MOON.SUN.SEP' in card:
                        moonSunSep.append(float(card.split('=')[-1]))
                    elif 'SKYMODEL.TARGET.AIRMASS' in card:
                        #moonAM.append( 1./np.cos(np.radians(90.-float(card.split('=')[-1]))) )
                        moonAM.append( float(card.split('=')[-1]) )
                    elif 'SKYMODEL.MOON.TARGET.SEP' in card:
                        moonTargetSep.append(float(card.split('=')[-1]))
                    elif 'SKYMODEL.MOON.ALT' in card:
                        moonAlt.append(float(card.split('=')[-1]))
        except:
            print filename, ' Failed'

    import healpy as hp
    from lsst.sims.utils import haversine

    nside = 4
    lat, az = hp.pix2ang(nside, np.arange(hp.nside2npix(nside)))
    alt = np.pi/2.-lat
    airmass = 1./np.cos(np.pi/2.-alt)

    # Only need low airmass and then 1/2 to sky
    good = np.where( (az >= 0) & (az <= np.pi) & (airmass <=2.6) & (airmass >= 1.) )
    airmass = airmass[good]
    alt=alt[good]
    az = az[good]

    moonAM = np.array(moonAM)
    moonAlt = np.array(moonAlt)
    moonSunSep = np.array(moonSunSep)
    moonTargetSep = np.array(moonTargetSep)
    moonAzDiff = moonTargetSep*0
    targetAlt = np.pi/2.-np.arccos(1./moonAM)
    # Compute the azimuth difference given the moon-target-seperation
    # Let's just do a stupid loop:
    for i in np.arange(targetAlt.size):
        possibleDistances = haversine(0., np.radians(moonAlt[i]),  az, az*0+targetAlt[i])
        diff = np.abs(possibleDistances - np.radians(moonTargetSep[i]))
        good = np.where(diff == diff.min())
        moonAzDiff[i] = az[good][0]
        # ok, now I have an alt and az, I can convert that back to a healpix id.

        hpid.append(hp.ang2pix(nside, np.pi/2.-targetAlt[i], moonAzDiff[i]))

    nrec = moonAM.size
    nwave = moonWave.size

    dtype = [('hpid', 'int'),
             ('moonAltitude', 'float'),
             ('moonSunSep', 'float'),
             ('spectra', 'float', (nwave)), ('mags', 'float', (6))]
    moonSpectra = np.zeros(nrec, dtype=dtype)
    moonSpectra['hpid'] = hpid
    moonSpectra['moonAltitude'] = moonAlt
    moonSpectra['moonSunSep'] = moonSunSep
    moonSpectra['spectra'] = moonSpec

    hPlank = 6.626068e-27  # erg s
    cLight = 2.99792458e10 # cm/s

    # Convert spectra from ph/s/m2/micron/arcsec2 to erg/s/cm2/nm/arcsec2
    moonSpectra['spectra'] = moonSpectra['spectra']/(100.**2)*hPlank*cLight/(moonWave*1e-7)/1e3

    # Sort things since this might be helpful later
    moonSpectra.sort(order=['moonSunSep','moonAltitude', 'hpid'])

    # Crop off the incomplete ones

    good =np.where((moonSpectra['moonAltitude'] >= 0) & (moonSpectra['moonAltitude'] < 89.) )
    moonSpectra = moonSpectra[good]

    # Load LSST filters
    throughPath = os.path.join(getPackageDir('throughputs'),'baseline')
    keys = ['u','g','r','i','z','y']
    nfilt = len(keys)
    filters = {}
    for filtername in keys:
        bp = np.loadtxt(os.path.join(throughPath, 'filter_'+filtername+'.dat'),
                        dtype=zip(['wave','trans'],[float]*2 ))
        tempB = Bandpass()
        tempB.setBandpass(bp['wave'],bp['trans'])
        filters[filtername] = tempB

    filterWave = np.array([filters[f].calcEffWavelen()[0] for f in keys ])

    for i,spectrum in enumerate(moonSpectra['spectra']):
        tempSed = Sed()
        tempSed.setSED(moonWave,flambda=spectrum)
        for j,filtName in enumerate(keys):
            try:
                moonSpectra['mags'][i][j] = tempSed.calcMag(filters[filtName])
            except:
                pass

    nbreak=5
    nrec = np.size(moonSpectra)
    for i in np.arange(nbreak):
        np.savez(os.path.join(outDir,'moonSpectra_'+str(i)+'.npz'), wave = moonWave, spec=moonSpectra[i*nrec/nbreak:(i+1)*nrec/nbreak], filterWave=filterWave)

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
    def __init__(self, mags=False, darkSkyMags=None, fitResults=None):
        """
        Read the Solar spectrum into a handy object and compute mags in different filters
        mags:  If true, only return the LSST filter magnitudes, otherwise return the full spectrum

        darkSkyMags = dict of the zenith dark sky values to be assumed. The twilight fits are
        done relative to the dark sky level.
        fitResults = dict of twilight parameters based on twilightFunc. Keys should be filter names.
        """

        if darkSkyMags is None:
            darkSkyMags = {'u': 22.8, 'g': 22.3, 'r': 21.2,
                           'i': 20.3, 'z': 19.3, 'y': 18.0,
                           'B': 22.35, 'G': 21.71, 'R': 21.3}

        self.mags = mags

        dataDir = getPackageDir('sims_skybrightness_data')

        solarSaved = np.load(os.path.join(dataDir, 'solarSpec/solarSpec.npz'))
        self.solarSpec = Sed(wavelen=solarSaved['wave'], flambda=solarSaved['spec'])
        solarSaved.close()

        canonFilters = {}
        fnames = ['blue_canon.csv', 'green_canon.csv', 'red_canon.csv']

        # Filter names, from bluest to reddest.
        self.filterNames = ['B', 'G', 'R']

        for fname, filterName in zip(fnames, self.filterNames):
            bpdata = np.genfromtxt(os.path.join(dataDir, 'Canon/', fname), delimiter=', ',
                                   dtype=list(zip(['wave', 'through'], [float]*2)))
            bpTemp = Bandpass()
            bpTemp.setBandpass(bpdata['wave'], bpdata['through'])
            canonFilters[filterName] = bpTemp

        # Tack on the LSST filters
        throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
        lsstKeys = ['u', 'g', 'r', 'i', 'z', 'y']
        for key in lsstKeys:
            bp = np.loadtxt(os.path.join(throughPath, 'filter_'+key+'.dat'),
                            dtype=list(zip(['wave', 'trans'], [float]*2)))
            tempB = Bandpass()
            tempB.setBandpass(bp['wave'], bp['trans'])
            canonFilters[key] = tempB
            self.filterNames.append(key)

        # MAGIC NUMBERS from fitting the all-sky camera:
        # Code to generate values in sims_skybrightness/examples/fitTwiSlopesSimul.py
        # Which in turn uses twilight maps from sims_skybrightness/examples/buildTwilMaps.py
        # values are of the form:
        # 0: ratio of f^z_12 to f_dark^z
        # 1: slope of curve wrt sun alt
        # 2: airmass term (10^(arg[2]*(X-1)))
        # 3: azimuth term.
        # 4: zenith dark sky flux (erg/s/cm^2)

        # For z and y, just assuming the shape parameter fits are similar to the other bands.
        # Looks like the diode is not sensitive enough to detect faint sky.
        # Using the Patat et al 2006 I-band values for z and modeified a little for y as a temp fix.
        if fitResults is None:
            self.fitResults = {'B': [7.56765633e+00, 2.29798055e+01, 2.86879956e-01,
                                     3.01162143e-01, 2.58462036e-04],
                               'G': [2.38561156e+00, 2.29310648e+01, 2.97733083e-01,
                                     3.16403197e-01, 7.29660095e-04],
                               'R': [1.75498017e+00, 2.22011802e+01, 2.98619033e-01,
                                     3.28880254e-01, 3.24411056e-04],
                               'z': [2.29, 24.08, 0.3, 0.3, -666],
                               'y': [2.0, 24.08, 0.3, 0.3, -666]}

            # XXX-completely arbitrary fudge factor to make things brighter in the blue
            # Just copy the blue and say it's brighter.
            self.fitResults['u'] = [16., 2.29622121e+01, 2.85862729e-01,
                                    2.99902574e-01, 2.32325117e-04]
        else:
            self.fitResults = fitResults

        # Take out any filters that don't have fit results
        self.filterNames = [key for key in self.filterNames if key in self.fitResults]

        self.effWave = []
        self.solarMag = []
        for filterName in self.filterNames:
            self.effWave.append(canonFilters[filterName].calcEffWavelen()[0])
            self.solarMag.append(self.solarSpec.calcMag(canonFilters[filterName]))

        ord = np.argsort(self.effWave)
        self.filterNames = np.array(self.filterNames)[ord]
        self.effWave = np.array(self.effWave)[ord]
        self.solarMag = np.array(self.solarMag)[ord]

        # update the fit results to be zeropointed properly
        for key in self.fitResults:
            f0 = 10.**(-0.4*(darkSkyMags[key]-np.log10(3631.)))
            self.fitResults[key][-1] = f0

        self.solarWave = self.solarSpec.wavelen
        self.solarFlux = self.solarSpec.flambda
        # This one isn't as bad as the model grids, maybe we could get away with computing the magnitudes
        # in the __call__ each time.
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
# Crop off the incomplete ones

good = np.where((moonSpectra['moonAltitude'] >= 0) & (moonSpectra['moonAltitude'] < 89.))
moonSpectra = moonSpectra[good]


# Load LSST filters
throughPath = os.getenv('LSST_THROUGHPUTS_BASELINE')
keys = ['u', 'g', 'r', 'i', 'z', 'y']
nfilt = len(keys)
filters = {}
for filtername in keys:
    bp = np.loadtxt(os.path.join(throughPath, 'filter_'+filtername+'.dat'),
                    dtype=list(zip(['wave', 'trans'], [float]*2)))
    tempB = Bandpass()
    tempB.setBandpass(bp['wave'], bp['trans'])
    filters[filtername] = tempB

filterWave = np.array([filters[f].calcEffWavelen()[0] for f in keys])

for i, spectrum in enumerate(moonSpectra['spectra']):
    tempSed = Sed()
    tempSed.setSED(moonWave, flambda=spectrum)
    for j, filtName in enumerate(keys):
        try:
            moonSpectra['mags'][i][j] = tempSed.calcMag(filters[filtName])
        except:
            pass


nbreak = 5

# 需要导入模块: from lsst.sims.photUtils import Bandpass [as 别名]
# 或者: from lsst.sims.photUtils.Bandpass import setBandpass [as 别名]
wd = Sed()
wd.readSED_flambda(filepath)
magnorm = 16
fNorm = wd.calcFluxNorm(magnorm, imsimBand)
wd.multiplyFluxNorm(fNorm)

throughPath = os.path.join(getPackageDir('throughputs'), 'baseline')
lsstKeys = ['u', 'g', 'r', 'i', 'z', 'y']
# lsstKeys = ['r']
bps = {}
for key in lsstKeys:
    bp = np.loadtxt(os.path.join(throughPath, 'total_'+key+'.dat'),
                    dtype=zip(['wave', 'trans'], [float]*2))
    bpTemp = Bandpass()
    good = np.where(bp['trans'] > 0.)
    bpTemp.setBandpass(bp['wave'], bp['trans'], wavelen_min=bp['wave'][good].min(),
                       wavelen_max=bp['wave'][good].max())
    bps[key] = bpTemp

# Generate a GAIA observation
gaia_obs = SED2GAIA(wd)

# Read it in as an SED object
final_sed = ulysses2SED(data=gaia_obs, noisy=False, response=response)
print 're-binned mag, orignal mag, diff'
for key in lsstKeys:
    print key, final_sed.calcMag(bps[key]), wd.calcMag(bps[key]), \
        final_sed.calcMag(bps[key])-wd.calcMag(bps[key])