Python Sed.readSED_flambda方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
    def loadwdSEDs(self, subset = None):

        """
        By default will load all seds in wd directory. The user can also define a subset of
        what's in the directory and load only those SEDs instead. Will skip over extraneous
        files in sed folder.

        @param [in] subset is the list of the subset of files wanted if one doesn't want all files
        in the kurucz directory.

        @param [out] sedListH is the set of model SED spectra objects for Hydrogen WDs to be passed onto
        the matching routines.

        @param [out] sedListHE is the set of model SED spectra objects for Helium WDs to be passed onto
        the matching routines.
        """

        if self.wdDir is None:
            try:
                self.wdDir = str(self.sEDDir + '/' +
                                 self.specMapDict['wd'] + '/')
            except:
                raise ValueError(str('self.wdDir is None. ' +
                                     'Add path to wddirectory.'))


        files = []

        if subset is None:
            for fileName in os.listdir(self.wdDir):
                files.append(fileName)
        else:
            for fileName in subset:
                files.append(fileName)

        numFiles = len(files)
        numOn = 0

        sedListH = []
        sedListHE = []

        for fileName in files:
            if numOn % 100 == 0:
                print('Loading %i of %i: WD SEDs' % (numOn, numFiles))

            try:
                spec = Sed()
                spec.readSED_flambda(str(self.wdDir + '/' + fileName))
                spec.name = fileName
                if fileName.split("_")[1] == 'He':
                    sedListHE.append(spec)
                else:
                    sedListH.append(spec)

            except:
                continue

            numOn += 1

        return sedListH, sedListHE

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
 def _calcColors(self, sedname='C.dat'):
     """
     Calculate the colors for a moving object with sed 'sedname'.
     """
     # Do we need to read in the LSST bandpasses?
     try:
         self.lsst
     except AttributeError:
         filterdir = os.getenv('LSST_THROUGHPUTS_BASELINE')
         filterlist = ('u', 'g', 'r', 'i', 'z', 'y')
         self.lsst ={}
         for f in filterlist:
             self.lsst[f] = Bandpass()
             self.lsst[f].readThroughput(os.path.join(filterdir, 'total_'+f+'.dat'))
         self.vband = Bandpass()
         self.vband.readThroughput('harris_V.dat')
         self.colors = {}
     # See if the sed's colors are in memory already.
     if sedname not in self.colors:
         moSed = Sed()
         moSed.readSED_flambda(sedname)
         vmag = moSed.calcMag(self.vband)
         self.colors[sedname] = {}
         for f in filterlist:
             self.colors[sedname][f] = moSed.calcMag(self.lsst[f]) - vmag
     return self.colors[sedname]

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
    def testNoSystematicUncertainty(self):
        """
        Test that systematic uncertainty is handled correctly when set to None.
        """
        m5 = [23.5, 24.3, 22.1, 20.0, 19.5, 21.7]
        photParams = PhotometricParameters(sigmaSys=0.0)

        obs_metadata = ObservationMetaData(
            unrefractedRA=23.0, unrefractedDec=45.0, m5=m5, bandpassName=self.filterNameList
        )

        magnitudes = []
        for bp in self.bpList:
            mag = self.starSED.calcMag(bp)
            magnitudes.append(mag)

        skySedList = []

        for bp, hardware, filterName in zip(self.bpList, self.hardwareList, self.filterNameList):
            skyDummy = Sed()
            skyDummy.readSED_flambda(os.path.join(lsst.utils.getPackageDir("throughputs"), "baseline", "darksky.dat"))
            normalizedSkyDummy = setM5(
                obs_metadata.m5[filterName],
                skyDummy,
                bp,
                hardware,
                seeing=LSSTdefaults().seeing(filterName),
                photParams=photParams,
            )

            skySedList.append(normalizedSkyDummy)

        sigmaList = snr.calcMagError_m5(numpy.array(magnitudes), numpy.array(self.bpList), numpy.array(m5), photParams)

        for i in range(len(self.bpList)):
            snrat = snr.calcSNR_sed(
                self.starSED,
                self.bpList[i],
                skySedList[i],
                self.hardwareList[i],
                seeing=LSSTdefaults().seeing(self.filterNameList[i]),
                photParams=PhotometricParameters(),
            )

            testSNR, gamma = snr.calcSNR_m5(
                numpy.array([magnitudes[i]]),
                [self.bpList[i]],
                numpy.array([m5[i]]),
                photParams=PhotometricParameters(sigmaSys=0.0),
            )

            self.assertAlmostEqual(
                snrat, testSNR[0], 10, msg="failed on calcSNR_m5 test %e != %e " % (snrat, testSNR[0])
            )

            control = snr.magErrorFromSNR(testSNR)

            msg = "%e is not %e; failed" % (sigmaList[i], control)

            self.assertAlmostEqual(sigmaList[i], control, 10, msg=msg)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [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 readSED_flambda [as 别名]
    def _getSedCopy(self, sedName):
        """
        Return a copy of the requested SED, either from the cached
        version or creating a new one and caching a copy for later
        reuse.
        """
        if sedName in self.uniqueSeds:
            # we have already read in this file; no need to do it again
            sed = Sed(wavelen=self.uniqueSeds[sedName].wavelen,
                      flambda=self.uniqueSeds[sedName].flambda,
                      fnu=self.uniqueSeds[sedName].fnu,
                      name=self.uniqueSeds[sedName].name)
        else:
            # load the SED of the object
            sed = Sed()
            sedFile = os.path.join(self.sedDir, sedName)
            sed.readSED_flambda(sedFile)

            flambdaCopy = copy.deepcopy(sed.flambda)

            #If the SED is zero inside of the bandpass, GalSim raises an error.
            #This sets a minimum flux value of 1.0e-30 so that the SED is never technically
            #zero inside of the bandpass.
            sed.flambda = numpy.array([ff if ff>1.0e-30 else 1.0e-30 for ff in flambdaCopy])
            sed.fnu = None

            #copy the unnormalized file to uniqueSeds so we don't have to read it in again
            sedCopy = Sed(wavelen=sed.wavelen, flambda=sed.flambda,
                          fnu=sed.fnu, name=sed.name)
            self.uniqueSeds[sedName] = sedCopy
        return sed

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [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 readSED_flambda [as 别名]
    def setUpClass(cls):
        lsstDefaults=LSSTdefaults()
        cls.dbName = 'uncertaintyTestDB.db'
        if os.path.exists(cls.dbName):
            os.unlink(cls.dbName)

        default_obs_metadata = makePhoSimTestDB(filename=cls.dbName, size=10, radius = 5.0)
        bandpass = ['u', 'g', 'r', 'i', 'z', 'y']
        m5 = lsstDefaults._m5.values()

        cls.obs_metadata = ObservationMetaData(
                                              pointingRA = default_obs_metadata.pointingRA,
                                              pointingDec = default_obs_metadata.pointingDec,
                                              rotSkyPos = default_obs_metadata.rotSkyPos,
                                              bandpassName = bandpass,
                                              m5 = m5
                                              )

        cls.obs_metadata.setBandpassM5andSeeing(bandpassName=bandpass, m5=m5)
        cls.driver = 'sqlite'
        cls.host = ''

        cls.skySeds = []
        cls.hardwareBandpasses = []
        cls.totalBandpasses = []
        cls.bandpasses = ['u', 'g', 'r', 'i', 'z', 'y']

        components = ['detector.dat', 'm1.dat', 'm2.dat', 'm3.dat',
                      'lens1.dat', 'lens2.dat', 'lens3.dat']

        for b in cls.bandpasses:
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughput(os.path.join(lsst.utils.getPackageDir('throughputs'),
                                                      'baseline', 'total_%s.dat' % b))
            cls.totalBandpasses.append(bandpassDummy)

        for b in cls.bandpasses:
            finalComponents = []
            for c in components:
                finalComponents.append(os.path.join(lsst.utils.getPackageDir('throughputs'), 'baseline', c))
            finalComponents.append(os.path.join(lsst.utils.getPackageDir('throughputs'), 'baseline', 'filter_%s.dat' %b))
            bandpassDummy = Bandpass()
            bandpassDummy.readThroughputList(finalComponents)
            cls.hardwareBandpasses.append(bandpassDummy)

        for i in range(len(cls.bandpasses)):
            sedDummy = Sed()
            sedDummy.readSED_flambda(os.path.join(lsst.utils.getPackageDir('throughputs'), 'baseline', 'darksky.dat'))
            normalizedSedDummy = setM5(cls.obs_metadata.m5[cls.bandpasses[i]], sedDummy,
                                       cls.totalBandpasses[i], cls.hardwareBandpasses[i],
                                       FWHMeff=lsstDefaults.FWHMeff(cls.bandpasses[i]),
                                       photParams=PhotometricParameters())

            cls.skySeds.append(normalizedSedDummy)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
    def loadBC03(self, subset = None):

        """
        This loads the Bruzual and Charlot SEDs that are currently in the SIMS_SED_LIBRARY.
        If the user wants to use different SEDs another loading method can be created and used in place
        of this.

        @param [in] subset is the list of the subset of files in the galDir that the user
        can specify if using all the SEDs in the directory is not desired.

        @param [out] sedList is the set of model SED spectra objects to be passed onto the matching routines.
        """

        if self.galDir is None:
            raise ValueError('self.galDir is None. Add path to galaxy directory.')

        files = []

        if subset is None:
            for fileName in os.listdir(self.galDir):
                files.append(fileName)
        else:
            for fileName in subset:
                files.append(fileName)

        numFiles = len(files)
        numOn = 0

        sedList = []

        for fileName in files:
            if numOn % 100 == 0:
                print('Loading %i of %i: BC Galaxy SEDs' % (numOn, numFiles))

            try:
                spec = Sed()
                spec.readSED_flambda(str(self.galDir + '/' + fileName))
                spec.name = fileName
                fileNameAsList = fileName.split('.')
                spec.type = fileNameAsList[0]
                spec.age = float(fileNameAsList[1])
                metallicity = fileNameAsList[2].split('Z')[0]
                #Final form is z/zSun
                spec.metallicity = float(metallicity) * (10 ** ((len(metallicity)-1)*-1))

            except:
                continue

            sedList.append(spec)

            numOn += 1

        return sedList

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
    def loadmltSEDs(self, subset = None):

        """
        By default will load all seds in mlt directory. The user can also define a subset of
        what's in the directory and load only those SEDs instead. Will skip over extraneous
        files in sed folder.

        @param [in] subset is the list of the subset of files wanted if one doesn't want all files
        in the mlt directory.

        @param [out] sedList is the set of model SED spectra objects to be passed onto the matching
        routines.
        """

        if self.mltDir is None:
            try:
                self.mltDir = str(self.sEDDir + '/' +
                                  self.specMapDict['mlt'] + '/')
            except:
                raise ValueError(str('self.mltDir is None. ' +
                                     'Add path to mlt directory.'))

        files = []

        if subset is None:
            for fileName in os.listdir(self.mltDir):
                files.append(fileName)
        else:
            for fileName in subset:
                files.append(fileName)

        numFiles = len(files)
        numOn = 0

        sedList = []

        for fileName in files:
            if numOn % 100 == 0:
                print('Loading %i of %i: MLT SEDs' % (numOn, numFiles))

            try:
                spec = Sed()
                spec.readSED_flambda(str(self.mltDir + '/' + fileName))
                spec.name = fileName

            except:
                continue

            sedList.append(spec)

            numOn += 1

        return sedList

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [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 readSED_flambda [as 别名]
def calcColors(sedname='C.dat'):
    # Calculate SSO colors.
    filterdir = os.getenv('LSST_THROUGHPUTS_BASELINE')
    filterlist = ('u', 'g', 'r', 'i', 'z', 'y')
    lsst ={}
    for f in filterlist:
        lsst[f] = Bandpass()
        lsst[f].readThroughput(os.path.join(filterdir, 'total_'+f+'.dat'))
    vband = Bandpass()
    vband.readThroughput('harris_V.dat')
    csed = Sed()
    csed.readSED_flambda(sedname)
    vmag = csed.calcMag(vband)
    dmags = {}
    for f in filterlist:
        dmags[f] = csed.calcMag(lsst[f]) - vmag
    return dmags

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
def setup_photUtils():
    throughputDir = os.getenv("LSST_THROUGHPUTS_DEFAULT")
    filterlist = ("u", "g", "r", "i", "z", "y")
    hardware = {}
    system = {}
    corecomponents = ["detector.dat", "lens1.dat", "lens2.dat", "lens3.dat", "m1.dat", "m2.dat", "m3.dat"]
    for f in filterlist:
        hardware[f] = Bandpass()
        system[f] = Bandpass()
        componentlist = copy.deepcopy(corecomponents)
        componentlist += ["filter_%s.dat" % f]
        hardware[f].readThroughputList(componentlist, rootDir=throughputDir)
        componentlist += ["atmos_10.dat"]
        system[f].readThroughputList(componentlist, rootDir=throughputDir)
    darksky = Sed()
    darksky.readSED_flambda(os.path.join(throughputDir, "darksky.dat"))
    return hardware, system, darksky

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [as 别名]
    def testNoSystematicUncertainty(self):
        """
        Test that systematic uncertainty is handled correctly when set to None.
        """
        m5_list = [23.5, 24.3, 22.1, 20.0, 19.5, 21.7]
        photParams= PhotometricParameters(sigmaSys=0.0)

        obs_metadata = ObservationMetaData(pointingRA=23.0, pointingDec=45.0,
                                           m5=m5_list, bandpassName=self.filterNameList)

        magnitude_list = []
        for bp in self.bpList:
            mag = self.starSED.calcMag(bp)
            magnitude_list.append(mag)

        skySedList = []

        for bp, hardware, filterName, mm, m5 in \
            zip(self.bpList, self.hardwareList, self.filterNameList, magnitude_list, m5_list):

            skyDummy = Sed()
            skyDummy.readSED_flambda(os.path.join(lsst.utils.getPackageDir('throughputs'),
                                     'baseline', 'darksky.dat'))

            normalizedSkyDummy = setM5(obs_metadata.m5[filterName], skyDummy,
                                       bp, hardware,
                                       FWHMeff=LSSTdefaults().FWHMeff(filterName),
                                       photParams=photParams)

            sigma, gamma = snr.calcMagError_m5(mm, bp, m5, photParams)


            snrat = snr.calcSNR_sed(self.starSED, bp, normalizedSkyDummy, hardware,
                              FWHMeff=LSSTdefaults().FWHMeff(filterName),
                              photParams=PhotometricParameters())

            testSNR, gamma = snr.calcSNR_m5(mm, bp, m5, photParams=PhotometricParameters(sigmaSys=0.0))

            self.assertAlmostEqual(snrat, testSNR, 10, msg = 'failed on calcSNR_m5 test %e != %e ' \
                                                               % (snrat, testSNR))

            control = snr.magErrorFromSNR(testSNR)

            msg = '%e is not %e; failed' % (sigma, control)

            self.assertAlmostEqual(sigma, control, 10, msg=msg)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import readSED_flambda [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 readSED_flambda [as 别名]
def read_quasar():
    # read quasar spectra and redshift
    homedir = os.getenv("HOME")
    quasardir = os.path.join(homedir, "seds/quasar")
    # read zero redshift quasar
    base = Sed()
    base.readSED_flambda(os.path.join(quasardir, "quasar.dat"))
    # redshift 
    #redshifts = [0, 0.1, 0.2, 0.3, 0.5, 0.8, 1.0, 1.3, 1.6, 1.9, 2.2, 2.5]
    #redshifts = numpy.array(redshifts)
    redshifts= numpy.arange(0, 2.8, 0.1)
    quasars = {}
    for z in redshifts:
        wavelen, flambda = base.redshiftSED(z, wavelen=base.wavelen, flambda=base.flambda)
        quasars[z] = Sed(wavelen=wavelen, flambda=flambda)
    print "# Generated %d quasars at redshifts between %f and %f" %(len(redshifts), redshifts.min(), redshifts.max())
    # resample onto the standard bandpass for Bandpass obj's and calculate fnu to speed later calculations
    for z in redshifts:
        quasars[z].synchronizeSED(wavelen_min=WMIN, wavelen_max=WMAX, wavelen_step=WSTEP)
    return quasars, redshifts