Python Sed.resampleSED方法代码示例

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
 def testSedBandpassMatch(self):
     """Test errors when bandpass and sed do not completely overlap in wavelength range."""
     # Test case where they do match (no error message)
     sedwavelen = np.arange(self.wmin, self.wmax+.5, 1)
     sedflambda = np.ones(len(sedwavelen))
     testsed = Sed(wavelen=sedwavelen, flambda=sedflambda)
     print('')
     # Test that no warning is made.
     with warnings.catch_warnings(record=True) as wa:
         w, f = testsed.resampleSED(wavelen_match=self.testbandpass.wavelen,
                                    wavelen=testsed.wavelen, flux=testsed.flambda)
         self.assertEqual(len(wa), 0)
     np.testing.assert_equal(w, testsed.wavelen)
     np.testing.assert_equal(f, testsed.flambda)
     # Test that warning is given for non-overlap at either top or bottom end of wavelength range.
     sedwavelen = np.arange(self.wmin, self.wmax - 50, 1)
     sedflambda = np.ones(len(sedwavelen))
     testsed = Sed(wavelen=sedwavelen, flambda=sedflambda)
     with warnings.catch_warnings(record=True) as wa:
         testsed.resampleSED(wavelen_match=self.testbandpass.wavelen)
         self.assertEqual(len(wa), 1)
         self.assertIn('non-overlap', str(wa[-1].message))
     np.testing.assert_equal(testsed.flambda[-1:], np.NaN)
     sedwavelen = np.arange(self.wmin+50, self.wmax, 1)
     sedflambda = np.ones(len(sedwavelen))
     testsed = Sed(wavelen=sedwavelen, flambda=sedflambda)
     with warnings.catch_warnings(record=True) as wa:
         testsed.resampleSED(wavelen_match=self.testbandpass.wavelen)
         self.assertEqual(len(wa), 1)
         self.assertIn('non-overlap', str(wa[-1].message))
     np.testing.assert_equal(testsed.flambda[0], np.NaN)
     np.testing.assert_equal(testsed.flambda[49], np.NaN)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def magListForSed(self, sedobj, indices=None):
        """
        Return a list of magnitudes for a single Sed object.

        @param [in] sedobj is an Sed object.  Its wavelength grid can be arbitrary.  If necessary,
        a copy will be created and resampled onto the wavelength grid of the Bandpasses before
        magnitudes are calculated.  The original Sed will be unchanged.

        @param [in] indices is an optional list of indices indicating which bandpasses to actually
        calculate magnitudes for.  Other magnitudes will be listed as numpy.NaN (i.e. this method will
        return as many magnitudes as were loaded with the loadBandpassesFromFiles methods; it will
        just return numpy.NaN for magnitudes you did not actually ask for)

        @param [out] magList is a list of magnitudes in the bandpasses stored in this BandpassDict
        """

        if sedobj.wavelen is not None:

            # If the Sed's wavelength grid agrees with self._wavelen_match to one part in
            # 10^6, just use the Sed as-is.  Otherwise, copy it and resample it onto
            # self._wavelen_match
            if sedobj._needResample(wavelen_match=self._wavelen_match):
                dummySed = Sed(wavelen=sedobj.wavelen, flambda=sedobj.flambda)
                dummySed.resampleSED(force=True, wavelen_match=self._bandpassDict.values()[0].wavelen)
            else:
                dummySed = sedobj

            return numpy.array(self._magListForSed(dummySed, indices=indices))

        else:
            return numpy.array([numpy.NaN]*len(self._bandpassDict))

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def calcMagNorm(self, objectMags, sedObj, bandpassDict, mag_error = None,
                    redshift = None, filtRange = None):

        """
        This will find the magNorm value that gives the closest match to the magnitudes of the object
        using the matched SED. Uses scipy.optimize.leastsq to find the values of fluxNorm that minimizes
        the function: ((flux_obs - (fluxNorm*flux_model))/flux_error)**2.

        @param [in] objectMags are the magnitude values for the object with extinction matching that of
        the SED object. In the normal case using the selectSED routines above it will be dereddened mags.

        @param [in] sedObj is an Sed class instance that is set with the wavelength and flux of the
        matched SED

        @param [in] bandpassDict is a BandpassDict class instance with the Bandpasses set to those
        for the magnitudes given for the catalog object

        @param [in] mag_error are provided error values for magnitudes in objectMags. If none provided
        then this defaults to 1.0. This should be an array of the same length as objectMags.

        @param [in] redshift is the redshift of the object if the magnitude is observed

        @param [in] filtRange is a selected range of filters specified by their indices in the bandpassList
        to match up against. Used when missing data in some magnitude bands.

        @param [out] bestMagNorm is the magnitude normalization for the given magnitudes and SED
        """

        import scipy.optimize as opt

        sedTest = Sed()
        sedTest.setSED(sedObj.wavelen, flambda = sedObj.flambda)
        if redshift is not None:
            sedTest.redshiftSED(redshift)
        imSimBand = Bandpass()
        imSimBand.imsimBandpass()
        zp = -2.5*np.log10(3631)  #Note using default AB zeropoint
        flux_obs = np.power(10,(objectMags + zp)/(-2.5))
        sedTest.resampleSED(wavelen_match=bandpassDict.wavelenMatch)
        sedTest.flambdaTofnu()
        flux_model = sedTest.manyFluxCalc(bandpassDict.phiArray, bandpassDict.wavelenStep)
        if filtRange is not None:
            flux_obs = flux_obs[filtRange]
            flux_model = flux_model[filtRange]
        if mag_error is None:
            flux_error = np.ones(len(flux_obs))
        else:
            flux_error = np.abs(flux_obs*(np.log(10)/(-2.5))*mag_error)
        bestFluxNorm = opt.leastsq(lambda x: ((flux_obs - (x*flux_model))/flux_error), 1.0)[0][0]
        sedTest.multiplyFluxNorm(bestFluxNorm)
        bestMagNorm = sedTest.calcMag(imSimBand)
        return bestMagNorm

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def fluxListForSed(self, sedobj, indices=None):
        """
        Return a list of Fluxes for a single Sed object.

        @param [in] sedobj is an Sed object.   Its wavelength grid can be arbitrary. If necessary,
        a copy will be created and resampled onto the wavelength grid of the Bandpasses before
        fluxes are calculated.  The original Sed will be unchanged.

        @param [in] indices is an optional list of indices indicating which bandpasses to actually
        calculate fluxes for.  Other fluxes will be listed as numpy.NaN (i.e. this method will
        return as many fluxes as were loaded with the loadBandpassesFromFiles methods; it will
        just return numpy.NaN for fluxes you did not actually ask for)

        @param [out] fluxList is a list of fluxes in the bandpasses stored in this BandpassDict

        Note on units: Fluxes calculated this way will be the flux density integrated over the
        weighted response curve of the bandpass.  See equaiton 2.1 of the LSST Science Book

        http://www.lsst.org/scientists/scibook
        """

        if sedobj.wavelen is not None:

            # If the Sed's wavelength grid agrees with self._wavelen_match to one part in
            # 10^6, just use the Sed as-is.  Otherwise, copy it and resample it onto
            # self._wavelen_match
            if sedobj._needResample(wavelen_match=self._wavelen_match):
                dummySed = Sed(wavelen=sedobj.wavelen, flambda=sedobj.flambda)
                dummySed.resampleSED(force=True, wavelen_match=self._bandpassDict.values()[0].wavelen)
            else:
                dummySed = sedobj

            return numpy.array(self._fluxListForSed(dummySed, indices=indices))

        else:
            return numpy.array([numpy.NaN]*len(self._bandpassDict))

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def applyIGM(self, redshift, sedobj):

        """
        Apply IGM extinction to already redshifted sed with redshift
        between zMin and zMax defined by range of lookup tables

        @param [in] redshift is the redshift of the incoming SED object

        @param [in] sedobj is the SED object to which IGM extinction will be applied. This object
        will be modified as a result of this.
        """

        if self.IGMisInitialized == False:
            self.initializeIGM()

        #First make sure redshift is in range of lookup tables.
        if (redshift < self.zMin) or (redshift > self.zMax):
            warnings.warn(str("IGM Lookup tables only applicable for " + str(self.zMin) + " < z < " + str(self.zMax) + ". No action taken"))
            return

        #Now read in closest two lookup tables for given redshift
        lowerSed = Sed()
        upperSed = Sed()
        for lower, upper in zip(self.zRange[:-1], self.zRange[1:]):
            if lower <= redshift <= upper:
                lowerSed.setSED(self.meanLookups['%.1f' % lower][:,0],
                                flambda = self.meanLookups['%.1f' % lower][:,1])
                upperSed.setSED(self.meanLookups['%.1f' % upper][:,0],
                                flambda = self.meanLookups['%.1f' % lower][:,1])
                break

        #Redshift lookup tables to redshift of source, i.e. if source redshift is 1.78 shift lookup
        #table for 1.7 and lookup table for 1.8 to up and down to 1.78, respectively
        zLowerShift = ((1.0 + redshift)/(1.0 + lower)) - 1.0
        zUpperShift = ((1.0 + redshift)/(1.0 + upper)) - 1.0
        lowerSed.redshiftSED(zLowerShift)
        upperSed.redshiftSED(zUpperShift)

        #Resample lower and upper transmission data onto same wavelength grid.
        minWavelen = 300. #All lookup tables are usable above 300nm
        maxWavelen = np.amin([lowerSed.wavelen[-1], upperSed.wavelen[-1]]) - 0.01
        lowerSed.resampleSED(wavelen_min = minWavelen, wavelen_max = maxWavelen, wavelen_step = 0.01)
        upperSed.resampleSED(wavelen_match = lowerSed.wavelen)

        #Now insert this into a transmission array of 1.0 beyond the limits of current application
        #So that we can get an sed back that extends to the longest wavelengths of the incoming sed
        finalWavelen = np.arange(300., sedobj.wavelen[-1]+0.01, 0.01)
        finalFlambdaExtended = np.ones(len(finalWavelen))

        #Weighted Average of Transmission from each lookup table to get final transmission
        #table at desired redshift
        dzGrid = self.zDelta #Step in redshift between transmission lookup table files
        finalSed = Sed()
        finalFlambda = (lowerSed.flambda*(1.0 - ((redshift - lower)/dzGrid)) +
                        upperSed.flambda*(1.0 - ((upper - redshift)/dzGrid)))
        finalFlambdaExtended[0:len(finalFlambda)] = finalFlambda
        finalSed.setSED(wavelen = finalWavelen, flambda = finalFlambdaExtended)

        #Resample incoming sed to new grid so that we don't get warnings from multiplySED
        #about matching wavelength grids
        sedobj.resampleSED(wavelen_match=finalSed.wavelen)

        #Now multiply transmission curve by input SED to get final result and make it the new flambda
        #data in the original sed which also is now on a new grid starting at 300 nm
        test = sedobj.multiplySED(finalSed)
        sedobj.flambda = test.flambda

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
def get_TotalMags(result, bandpasses=('u','g','r','i','z','y')):
    datadir = os.environ.get("SIMS_SED_LIBRARY_DIR")
    tpath = os.getenv('LSST_THROUGHPUTS_DEFAULT')
    bands = {"u":None, "g":None, "r":None, "i":None, "z":None, "y":None}
    for k in bands:
        bands[k] = Bandpass()
        bands[k].readThroughput(os.path.join(tpath, "total_%s.dat"%k))
    # Set up phi, the wavelength-normalized system response for each filter,
    # for each bandpass for manyMagCalc method.
    bplist = []
    for f in ['u','g','r','i','z','y']:
        bands[f].sbTophi()
        bplist.append(bands[f])
    ids = result['galid']
    diskfile = result['sedFilenameDisk']
    bulgefile = result['sedFilenameBulge']
    agnfile = result['sedFilenameAgn']

    diskmn = result['magNormDisk']
    bulgemn = result['magNormBulge']
    agnmn = result['magNormAgn']

    bulgeAv = result['internalAvBulge']
    diskAv = result['internalAvDisk']

    redshift = result['redshift']

    imsimband = Bandpass()
    imsimband.imsimBandpass()
    sedDict = {}
    retMags = dict([(k, []) for k in bands])
    a_int = None
    b_int = None
    tmpwavelen = None
    for id, df, dm, dav, bf, bm, bav, af, am, z in zip(ids, diskfile, diskmn, diskAv, 
            bulgefile, bulgemn, bulgeAv, agnfile, agnmn, redshift):
        tmpflux = None
        for comp in ((df, dm, dav, 'galaxySED', False), (bf, bm, bav, 'galaxySED', False), (af, am, None, 'agnSED', True)):
        #Zero out the AGN contribution
        #for comp in ((df, dm, dav, 'galaxySED', False), (bf, bm, bav, 'galaxySED', False), (af, 99.99, None, 'agnSED', True)):
            if not comp[0] == u'None':
                if sedDict.has_key(comp[0]):
                    sed = copy.deepcopy(sedDict[comp[0]])
                else:
                    sed = Sed()
                    print os.path.join(datadir,comp[3],comp[0])
                    sed.readSED_flambda(os.path.join(datadir,comp[3],comp[0]))
		    if comp[4]:
		        sed.resampleSED(wavelen_match=tmpwavelen)
                    sedDict[comp[0]] = sed
                if a_int is None:
                    phiarray, dlambda = sed.setupPhiArray(bplist)
                    a_int, b_int = sed.setupCCMab()
		    #Careful, this assumes that a disk or bulge sed is read
		    #before any agn sed
		    tmpwavelen = sed.wavelen
                fNorm = sed.calcFluxNorm(comp[1], imsimband)
                sed.multiplyFluxNorm(fNorm)
                #I guess this assumes rv=3.1??
                if comp[2]:
                    sed.addCCMDust(a_int, b_int, A_v=comp[2])
		wavelenArr=sed.wavelen
		if tmpflux is None:
		    tmpflux = sed.flambda
		else:
	            tmpflux += sed.flambda
	newgal = Sed(wavelen=wavelenArr, flambda=tmpflux)
        #a_mw, b_mw = sed.setupCCMab()
        #sed.addCCMDust(a_mw, b_mw, A_v=mwav)
        newgal.redshiftSED(z, dimming=True)
	newgal.resampleSED(wavelen_match=bplist[0].wavelen)
	newgal.flambdaTofnu()
        mags = newgal.manyMagCalc(phiarray, dlambda)
        for i,k in enumerate(['u','g','r','i','z','y']):
            retMags[k].append(mags[i])
    return retMags

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
import numpy as np
import lsst.sims.photUtils.Sed as Sed
import os

dataDir = os.getenv('SIMS_SKYBRIGHTNESS_DATA_DIR')

data = np.genfromtxt(os.path.join(dataDir, 'solarSpec/solarSpec.dat'), dtype=zip(['microns','Irr'],[float]*2))
#data['Irr'] = data['Irr']*1 #convert W/m2/micron to erg/s/cm2/nm (HA, it's the same!)

sun = Sed()
sun.setSED(data['microns']*1e3, flambda=data['Irr'])

# Match the wavelenth spacing and range to the ESO spectra
airglowSpec = np.load(os.path.join(dataDir, 'ESO_Spectra/Airglow/airglowSpectra.npz'))
sun.resampleSED(wavelen_match=airglowSpec['wave'])

np.savez(os.path.join(dataDir,'solarSpec/solarSpec.npz'), wave=sun.wavelen, spec=sun.flambda)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def testGalaxyPhotometricUncertainties(self):
        """
        Test in the case of a catalog of galaxies
        """
        lsstDefaults = LSSTdefaults()
        phot = PhotometryGalaxies()
        galDB = testGalaxyTileDBObj(driver=self.driver, host=self.host, database=self.dbName)
        galCat = testGalaxyCatalog(galDB, obs_metadata=self.obs_metadata)
        imsimband = Bandpass()
        imsimband.imsimBandpass()
        ct = 0
        for line in galCat.iter_catalog():
            bulgeSedName = line[50]
            diskSedName = line[51]
            agnSedName = line[52]
            magNormBulge = line[53]
            magNormDisk = line[54]
            magNormAgn = line[55]
            avBulge = line[56]
            avDisk = line[57]
            redshift = line[58]

            bulgeSed = Sed()
            bulgeSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                     defaultSpecMap[bulgeSedName]))
            fNorm=bulgeSed.calcFluxNorm(magNormBulge, imsimband)
            bulgeSed.multiplyFluxNorm(fNorm)

            diskSed = Sed()
            diskSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                    defaultSpecMap[diskSedName]))
            fNorm = diskSed.calcFluxNorm(magNormDisk, imsimband)
            diskSed.multiplyFluxNorm(fNorm)

            agnSed = Sed()
            agnSed.readSED_flambda(os.path.join(lsst.utils.getPackageDir('sims_sed_library'),
                                   defaultSpecMap[agnSedName]))
            fNorm = agnSed.calcFluxNorm(magNormAgn, imsimband)
            agnSed.multiplyFluxNorm(fNorm)

            a_int, b_int = bulgeSed.setupCCMab()
            bulgeSed.addCCMDust(a_int, b_int, A_v=avBulge)

            a_int, b_int = diskSed.setupCCMab()
            diskSed.addCCMDust(a_int, b_int, A_v=avDisk)

            bulgeSed.redshiftSED(redshift, dimming=True)
            diskSed.redshiftSED(redshift, dimming=True)
            agnSed.redshiftSED(redshift, dimming=True)

            bulgeSed.resampleSED(wavelen_match=self.totalBandpasses[0].wavelen)
            diskSed.resampleSED(wavelen_match=bulgeSed.wavelen)
            agnSed.resampleSED(wavelen_match=bulgeSed.wavelen)

            numpy.testing.assert_almost_equal(bulgeSed.wavelen, diskSed.wavelen)
            numpy.testing.assert_almost_equal(bulgeSed.wavelen, agnSed.wavelen)

            fl = bulgeSed.flambda + diskSed.flambda + agnSed.flambda

            totalSed = Sed(wavelen=bulgeSed.wavelen, flambda=fl)

            sedList = [totalSed, bulgeSed, diskSed, agnSed]

            for i, spectrum in enumerate(sedList):
                if i==0:
                    msgroot = 'failed on total'
                elif i==1:
                    msgroot = 'failed on bulge'
                elif i==2:
                    msgroot = 'failed on disk'
                elif i==3:
                    msgroot = 'failed on agn'

                for j, b in enumerate(self.bandpasses):
                    controlSigma = calcMagError_sed(spectrum, self.totalBandpasses[j],
                                             self.skySeds[j],
                                             self.hardwareBandpasses[j],
                                             FWHMeff=lsstDefaults.FWHMeff(b),
                                             photParams=PhotometricParameters())

                    testSigma = line[26+(i*6)+j]
                    msg = '%e neq %e; ' % (testSigma, controlSigma) + msgroot
                    self.assertAlmostEqual(testSigma, controlSigma, 10, msg=msg)
                    ct += 1

        self.assertGreater(ct, 0)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def testFlush(self):
        """
        Test that the flush method of SedList behaves properly
        """
        imsimBand = Bandpass()
        imsimBand.imsimBandpass()
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                 redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                 wavelenMatch=wavelen_match)

        self.assertEqual(len(testList), nSed)
        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

        for ix in range(len(sedNameList_0)):
            self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
            self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
            self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)

        for ix, (name, norm, iav, gav, zz) in \
        enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                      galacticAvList_0, redshiftList_0)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)



        testList.flush()

        sedNameList_1 = self.getListOfSedNames(nSed/2)
        magNormList_1 = numpy.random.random_sample(nSed/2)*5.0 + 15.0
        internalAvList_1 = numpy.random.random_sample(nSed/2)*0.3 + 0.1
        redshiftList_1 = numpy.random.random_sample(nSed/2)*5.0
        galacticAvList_1 = numpy.random.random_sample(nSed/2)*0.3 + 0.1

        testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                  internalAvList=internalAvList_1,
                                  galacticAvList=galacticAvList_1,
                                  redshiftList=redshiftList_1)

        self.assertEqual(len(testList), nSed/2)
        self.assertEqual(len(testList.redshiftList), nSed/2)
        self.assertEqual(len(testList.internalAvList), nSed/2)
        self.assertEqual(len(testList.galacticAvList), nSed/2)
        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)


        for ix in range(len(sedNameList_1)):
            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix], 10)
            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix], 10)
            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix], 10)

        for ix, (name, norm, iav, gav, zz) in \
        enumerate(zip(sedNameList_1, magNormList_1, internalAvList_1, \
                      galacticAvList_1, redshiftList_1)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
#.........这里部分代码省略.........

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

#.........这里部分代码省略.........
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))
            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=av)

            sedControl.redshiftSED(zz, dimming=False)

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)


        ################ now add galacticAv
        sedNameList = self.getListOfSedNames(nSed)
        magNormList = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList = numpy.random.random_sample(nSed)*5.0
        galacticAvList = numpy.random.random_sample(nSed)*0.3 + 0.1
        testList = SedList(sedNameList, magNormList, internalAvList=internalAvList,
                                 redshiftList=redshiftList, galacticAvList=galacticAvList)
        self.assertTrue(testList.wavelenMatch is None)
        self.assertTrue(testList.cosmologicalDimming is True)
        for avControl, avTest in zip(internalAvList, testList.internalAvList):
            self.assertAlmostEqual(avControl, avTest, 10)

        for zControl, zTest in zip(redshiftList, testList.redshiftList):
            self.assertAlmostEqual(zControl, zTest, 10)

        for avControl, avTest in zip(galacticAvList, testList.galacticAvList):
            self.assertAlmostEqual(avControl, avTest, 10)

        for name, norm, av, zz, gav, sedTest in \
        zip(sedNameList, magNormList, internalAvList, \
            redshiftList, galacticAvList, testList):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))
            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=av)

            sedControl.redshiftSED(zz, dimming=True)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)


        ################ now use a wavelen_match
        sedNameList = self.getListOfSedNames(nSed)
        magNormList = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList = numpy.random.random_sample(nSed)*5.0
        galacticAvList = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList, magNormList, internalAvList=internalAvList,
                                 redshiftList=redshiftList, galacticAvList=galacticAvList,
                                 wavelenMatch=wavelen_match)

        self.assertTrue(testList.cosmologicalDimming is True)
        for avControl, avTest in zip(internalAvList, testList.internalAvList):
            self.assertAlmostEqual(avControl, avTest, 10)

        for zControl, zTest in zip(redshiftList, testList.redshiftList):
            self.assertAlmostEqual(zControl, zTest, 10)

        for avControl, avTest in zip(galacticAvList, testList.galacticAvList):
            self.assertAlmostEqual(avControl, avTest, 10)

        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

        for name, norm, av, zz, gav, sedTest in \
        zip(sedNameList, magNormList, internalAvList, \
            redshiftList, galacticAvList, testList):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=av)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def testAlternateNormalizingBandpass(self):
        """
        A reiteration of testAddingToList, but testing with a non-imsimBandpass
        normalizing bandpass
        """
        normalizingBand = Bandpass()
        normalizingBand.readThroughput(os.path.join(getPackageDir('throughputs'),'baseline','total_r.dat'))
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0,
                           normalizingBandpass=normalizingBand,
                           internalAvList=internalAvList_0,
                           redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                           wavelenMatch=wavelen_match)


        sedNameList_1 = self.getListOfSedNames(nSed)
        magNormList_1 = numpy.random.random_sample(nSed)*5.0 + 15.0

        internalAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1

        redshiftList_1 = numpy.random.random_sample(nSed)*5.0

        galacticAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1


        testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                  internalAvList=internalAvList_1,
                                  galacticAvList=galacticAvList_1,
                                  redshiftList=redshiftList_1)

        self.assertEqual(len(testList), 2*nSed)
        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

        for ix in range(len(sedNameList_0)):
            self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
            self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
            self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)


        for ix in range(len(sedNameList_1)):
            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix+nSed], 10)
            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix+nSed], 10)
            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix+nSed], 10)

        for ix, (name, norm, iav, gav, zz) in \
          enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                     galacticAvList_0, redshiftList_0)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, normalizingBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)

        for ix, (name, norm, iav, gav, zz) in \
            enumerate(zip(sedNameList_1, magNormList_1, internalAvList_1, \
                          galacticAvList_1, redshiftList_1)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, normalizingBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)

            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix+nSed]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def testAddingNonesToList(self):
        """
        Test what happens if you add SEDs to an SedList that have None for
        one or more of the physical parameters (i.e. galacticAv, internalAv, or redshift)
        """
        imsimBand = Bandpass()
        imsimBand.imsimBandpass()
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                 redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                 wavelenMatch=wavelen_match)


        sedNameList_1 = self.getListOfSedNames(nSed)
        magNormList_1 = list(numpy.random.random_sample(nSed)*5.0 + 15.0)
        internalAvList_1 = list(numpy.random.random_sample(nSed)*0.3 + 0.1)
        redshiftList_1 = list(numpy.random.random_sample(nSed)*5.0)
        galacticAvList_1 = list(numpy.random.random_sample(nSed)*0.3 + 0.1)

        internalAvList_1[0] = None
        redshiftList_1[1] = None
        galacticAvList_1[2] = None

        internalAvList_1[3] = None
        redshiftList_1[3] = None

        internalAvList_1[4] = None
        galacticAvList_1[4] = None

        redshiftList_1[5] = None
        galacticAvList_1[5] = None

        internalAvList_1[6] = None
        redshiftList_1[6] = None
        galacticAvList_1[6] = None

        testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                  internalAvList=internalAvList_1,
                                  galacticAvList=galacticAvList_1,
                                  redshiftList=redshiftList_1)

        self.assertEqual(len(testList), 2*nSed)
        numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

        for ix in range(len(sedNameList_0)):
            self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
            self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
            self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)


        for ix in range(len(sedNameList_1)):
            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix+nSed], 10)
            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix+nSed], 10)
            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix+nSed], 10)

        for ix, (name, norm, iav, gav, zz) in \
        enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                      galacticAvList_0, redshiftList_0)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

            sedControl.redshiftSED(zz, dimming=True)
            sedControl.resampleSED(wavelen_match=wavelen_match)

            a_coeff, b_coeff = sedControl.setupCCMab()
            sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

            sedTest = testList[ix]

            numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
            numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
            numpy.testing.assert_array_equal(sedControl.fnu, sedTest.fnu)


        for ix, (name, norm, iav, gav, zz) in \
        enumerate(zip(sedNameList_1, magNormList_1, internalAvList_1, \
                      galacticAvList_1, redshiftList_1)):

            sedControl = Sed()
            sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

            fnorm = sedControl.calcFluxNorm(norm, imsimBand)
            sedControl.multiplyFluxNorm(fnorm)

            if iav is not None:
                a_coeff, b_coeff = sedControl.setupCCMab()
                sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
    def testAddingToList(self):
        """
        Test that we can add Seds to an already instantiated SedList
        """
        imsimBand = Bandpass()
        imsimBand.imsimBandpass()
        nSed = 10
        sedNameList_0 = self.getListOfSedNames(nSed)
        magNormList_0 = numpy.random.random_sample(nSed)*5.0 + 15.0
        internalAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        redshiftList_0 = numpy.random.random_sample(nSed)*5.0
        galacticAvList_0 = numpy.random.random_sample(nSed)*0.3 + 0.1
        wavelen_match = numpy.arange(300.0, 1500.0, 10.0)
        testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                 redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                 wavelenMatch=wavelen_match)


        # experiment with adding different combinations of physical parameter lists
        # as None and not None
        for addIav in [True, False]:
            for addRedshift in [True, False]:
                for addGav in [True, False]:

                    testList = SedList(sedNameList_0, magNormList_0, internalAvList=internalAvList_0, \
                                             redshiftList=redshiftList_0, galacticAvList=galacticAvList_0,
                                             wavelenMatch=wavelen_match)

                    sedNameList_1 = self.getListOfSedNames(nSed)
                    magNormList_1 = numpy.random.random_sample(nSed)*5.0 + 15.0

                    if addIav:
                        internalAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1
                    else:
                        internalAvList_1 = None

                    if addRedshift:
                        redshiftList_1 = numpy.random.random_sample(nSed)*5.0
                    else:
                        redshiftList_1 = None

                    if addGav:
                        galacticAvList_1 = numpy.random.random_sample(nSed)*0.3 + 0.1
                    else:
                        galacticAvList_1 = None


                    testList.loadSedsFromList(sedNameList_1, magNormList_1,
                                              internalAvList=internalAvList_1,
                                              galacticAvList=galacticAvList_1,
                                              redshiftList=redshiftList_1)

                    self.assertEqual(len(testList), 2*nSed)
                    numpy.testing.assert_array_equal(wavelen_match, testList.wavelenMatch)

                    for ix in range(len(sedNameList_0)):
                        self.assertAlmostEqual(internalAvList_0[ix], testList.internalAvList[ix], 10)
                        self.assertAlmostEqual(galacticAvList_0[ix], testList.galacticAvList[ix], 10)
                        self.assertAlmostEqual(redshiftList_0[ix], testList.redshiftList[ix], 10)


                    for ix in range(len(sedNameList_1)):
                        if addIav:
                            self.assertAlmostEqual(internalAvList_1[ix], testList.internalAvList[ix+nSed], 10)
                        else:
                            self.assertTrue(testList.internalAvList[ix+nSed] is None)

                        if addGav:
                            self.assertAlmostEqual(galacticAvList_1[ix], testList.galacticAvList[ix+nSed], 10)
                        else:
                            self.assertTrue(testList.galacticAvList[ix+nSed] is None)

                        if addRedshift:
                            self.assertAlmostEqual(redshiftList_1[ix], testList.redshiftList[ix+nSed], 10)
                        else:
                            self.assertTrue(testList.redshiftList[ix+nSed] is None)

                    for ix, (name, norm, iav, gav, zz) in \
                      enumerate(zip(sedNameList_0, magNormList_0, internalAvList_0, \
                                  galacticAvList_0, redshiftList_0)):

                        sedControl = Sed()
                        sedControl.readSED_flambda(os.path.join(self.sedDir, name+'.gz'))

                        fnorm = sedControl.calcFluxNorm(norm, imsimBand)
                        sedControl.multiplyFluxNorm(fnorm)

                        a_coeff, b_coeff = sedControl.setupCCMab()
                        sedControl.addCCMDust(a_coeff, b_coeff, A_v=iav)

                        sedControl.redshiftSED(zz, dimming=True)
                        sedControl.resampleSED(wavelen_match=wavelen_match)

                        a_coeff, b_coeff = sedControl.setupCCMab()
                        sedControl.addCCMDust(a_coeff, b_coeff, A_v=gav)

                        sedTest = testList[ix]

                        numpy.testing.assert_array_equal(sedControl.wavelen, sedTest.wavelen)
                        numpy.testing.assert_array_equal(sedControl.flambda, sedTest.flambda)
#.........这里部分代码省略.........

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
def ulysses2SED(data=None, workdir='output', wavefile='Ulysses_GaiaBPRP_meanSpecWavelength.txt',
                specfile='Ulysses_GaiaBPRP_noiseFreeSpectra.txt',
                noiseRoot='Ulysses_GaiaBPRP_noisyPhotSpec', noisy=True,
                response=None, wavelen_step=1.0, switch=675.):
    """
    Read in some ulysses output and return a single SED object.
    """
    if data is None:
        data = read_ulysses(workdir=workdir, wavefile=wavefile, specfile=specfile,
                            noiseRoot=noiseRoot)
    if response is None:
        response = gaia_response()

    if noisy:
        datakey = 'noisySpec'
        key2 = 'NoisySpec'
        red_spec = response.apply(data[datakey][0]['RP'+key2], blue=False)
        blue_spec = response.apply(data[datakey][0]['BP'+key2], blue=True)
    else:
        datakey = 'noiseFreeSpec'
        key2 = 'NoiseFreeSpec'
        red_spec = response.apply(data[datakey]['RP'+key2], blue=False)
        blue_spec = response.apply(data[datakey]['BP'+key2], blue=True)

    red_sed = Sed(wavelen=data['RP_wave'], flambda=red_spec)  # * 1e3)
    blue_sed = Sed(wavelen=data['BP_wave'], flambda=blue_spec)  # * 1e3)

    wavelen_min = data['BP_wave'].min()
    wavelen_max = data['RP_wave'].max()

    # Rebin the red and blue to a common wavelength array
    red_sed.resampleSED(wavelen_min=wavelen_min, wavelen_max=wavelen_max, wavelen_step=wavelen_step)
    blue_sed.resampleSED(wavelen_min=wavelen_min, wavelen_max=wavelen_max, wavelen_step=wavelen_step)

    # Clean up nan's from resampling
    red_sed.flambda[np.isnan(red_sed.flambda)] = 0.
    blue_sed.flambda[np.isnan(blue_sed.flambda)] = 0.

    # Assume Poisson stats for the noise.
    red_weight = np.ones(red_sed.flambda.size, dtype=float)  # 1./red_sed.flambda
    blue_weight = np.ones(red_sed.flambda.size, dtype=float)  # 1./blue_sed.flambda

    red_weight[np.where(red_sed.flambda == 0)] = 0.
    blue_weight[np.where(blue_sed.flambda == 0)] = 0.

    # red_cutoff = np.min(response.red_wavelen[np.where(response.red_response == 0.)])
    # blue_cutoff = np.max(response.blue_wavelen[np.where(response.blue_response == 0.)])

    #red_weight[np.where(red_sed.wavelen < red_cutoff+cutoff_pad)] = 0.
    #blue_weight[np.where(blue_sed.wavelen > blue_cutoff-cutoff_pad)] = 0.
    # Just make a simple switchover wavelength
    red_weight[np.where(red_sed.wavelen <= switch)] = 0.
    blue_weight[np.where(blue_sed.wavelen > switch)] = 0.

    # weight = np.zeros(red_sed.wavelen.size, dtype=float)
    # weight[np.where(red_sed.flambda > 0)] += 1
    # weight[np.where(blue_sed.flambda > 0)] += 1

    flambda = (red_sed.flambda*red_weight + blue_sed.flambda*blue_weight) / (red_weight + blue_weight)
    # flambda[np.where(weight == 0.)] = 0.

    finalSED = Sed(flambda=flambda, wavelen=red_sed.wavelen)
    return finalSED

# 需要导入模块: from lsst.sims.photUtils import Sed [as 别名]
# 或者: from lsst.sims.photUtils.Sed import resampleSED [as 别名]
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):
    galname = gallist[gal_name[i]]
    tmpgal = Sed(wavelen=gals[galname].wavelen, flambda=gals[galname].flambda)
    tmpgal.addCCMDust(a_int, b_int, ebv=ebv_int[i])
    tmpgal.redshiftSED(redshifts[i])
    tmpgal.resampleSED(wavelen_min=wavelen_min, wavelen_max=wavelen_max, wavelen_step=wavelen_step)
    tmpgal.addCCMDust(a_mw, b_mw, ebv=ebv_mw[i])
    tmpgal.multiplyFluxNorm(fluxnorm[i])
    #tmpgal.flambdaTofnu() #Not needed because multiplyFluxNorm calculates fnu
    tmpmags = tmpgal.manyMagCalc(phiarray, dlambda)
    j = 0
    for f in filterlist:
        mags2[f][i] = tmpmags[j]
        j = j+1
dt, t = dtime(t)
print "Calculating dust/redshift/dust/fluxnorm/%d magnitudes for %d galaxies optimized way took %f s" \
      %(len(filterlist), num_gal, dt)


# Check for differences in magnitudes.
import pylab