Python Sed.setupCCM_ab方法代碼示例

# 需要導入模塊: from lsst.sims.photUtils.Sed import Sed [as 別名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setupCCM_ab [as 別名]

#.........這裏部分代碼省略.........
            Using the dict assumes that the wavelength sampling and range
            is the same for all elements of the dict.

            if provided, overrides wavelen input and the SED is
            obtained at the wavelength values native to bandpass
            object.


        Returns
        -------
        `sims_photutils.sed` object containing the wavelengths and SED
        values from the SN at time time in units of ergs/cm^2/sec/nm


        .. note: If both wavelen and bandpassobject are `None` then exception,
                 will be raised.
        Examples
        --------
        >>> sed = SN.SNObjectSED(time=0., wavelen=wavenm)
        '''

        if wavelen is None and bandpass is None:
            raise ValueError('A non None input to either wavelen or\
                              bandpassobject must be provided')

        # if bandpassobject present, it overrides wavelen
        if bandpass is not None:
            if isinstance(bandpass, BandpassDict):
                firstfilter = bandpass.keys()[0]
                bp = bandpass[firstfilter]
            else:
                bp = bandpass
            # remember this is in nm
            wavelen = bp.wavelen

        flambda = np.zeros(len(wavelen))


        # self.mintime() and self.maxtime() are properties describing
        # the ranges of SNCosmo.Model in time. Behavior beyond this is 
        # determined by self.modelOutSideTemporalRange
        if (time >= self.mintime()) and (time <= self.maxtime()):
            # If SNCosmo is requested a SED value beyond the wavelength range
            # of model it will crash. Try to prevent that by returning np.nan for
            # such wavelengths. This will still not help band flux calculations
            # but helps us get past this stage.

            flambda = flambda * np.nan

            # Convert to Ang
            wave = wavelen * 10.0
            mask1 = wave >= self.minwave()
            mask2 = wave <= self.maxwave()
            mask = mask1 & mask2
            wave = wave[mask]

            # flux density dE/dlambda returned from SNCosmo in
            # ergs/cm^2/sec/Ang, convert to ergs/cm^2/sec/nm

            flambda[mask] = self.flux(time=time, wave=wave)
            flambda[mask] = flambda[mask] * 10.0
        else:
            # use prescription for modelOutSideTemporalRange
            if self.modelOutSideTemporalRange != 'zero':
                raise NotImplementedError('Model not implemented, change to zero\n')
                # Else Do nothing as flambda is already 0.
                # This takes precedence over being outside wavelength range
                
        if self.rectifySED:
            # Note that this converts nans into 0.
            flambda = np.where(flambda > 0., flambda, 0.)
        SEDfromSNcosmo = Sed(wavelen=wavelen, flambda=flambda)

        if not applyExtinction:
            return SEDfromSNcosmo

        # Apply LSST extinction
        global _sn_ax_cache
        global _sn_bx_cache
        global _sn_ax_bx_wavelen
        if _sn_ax_bx_wavelen is None \
        or len(wavelen)!=len(_sn_ax_bx_wavelen) \
        or (wavelen!=_sn_ax_bx_wavelen).any():

            ax, bx = SEDfromSNcosmo.setupCCM_ab()
            _sn_ax_cache = ax
            _sn_bx_cache = bx
            _sn_ax_bx_wavelen = np.copy(wavelen)
        else:
            ax = _sn_ax_cache
            bx = _sn_bx_cache

        if self.ebvofMW is None:
            raise ValueError('ebvofMW attribute cannot be None Type and must'
                             ' be set by hand using set_MWebv before this'
                             'stage, or by using setcoords followed by'
                             'mwEBVfromMaps\n')

        SEDfromSNcosmo.addDust(a_x=ax, b_x=bx, ebv=self.ebvofMW)
        return SEDfromSNcosmo

# 需要導入模塊: from lsst.sims.photUtils.Sed import Sed [as 別名]
# 或者: from lsst.sims.photUtils.Sed.Sed import setupCCM_ab [as 別名]
    def testAlternateBandpassesGalaxies(self):
        """
        the same as testAlternateBandpassesStars, but for galaxies
        """

        obs_metadata_pointed = ObservationMetaData(mjd=50000.0,
                                                   boundType='circle',
                                                   pointingRA=0.0, pointingDec=0.0,
                                                   boundLength=10.0)

        dtype = np.dtype([('galid', np.int),
                          ('ra', np.float),
                          ('dec', np.float),
                          ('uTotal', np.float),
                          ('gTotal', np.float),
                          ('rTotal', np.float),
                          ('iTotal', np.float),
                          ('zTotal', np.float),
                          ('uBulge', np.float),
                          ('gBulge', np.float),
                          ('rBulge', np.float),
                          ('iBulge', np.float),
                          ('zBulge', np.float),
                          ('uDisk', np.float),
                          ('gDisk', np.float),
                          ('rDisk', np.float),
                          ('iDisk', np.float),
                          ('zDisk', np.float),
                          ('uAgn', np.float),
                          ('gAgn', np.float),
                          ('rAgn', np.float),
                          ('iAgn', np.float),
                          ('zAgn', np.float),
                          ('bulgeName', str, 200),
                          ('bulgeNorm', np.float),
                          ('bulgeAv', np.float),
                          ('diskName', str, 200),
                          ('diskNorm', np.float),
                          ('diskAv', np.float),
                          ('agnName', str, 200),
                          ('agnNorm', np.float),
                          ('redshift', np.float)])

        test_cat = cartoonGalaxies(self.galaxy, obs_metadata=obs_metadata_pointed)
        with lsst.utils.tests.getTempFilePath('.txt') as catName:
            test_cat.write_catalog(catName)
            catData = np.genfromtxt(catName, dtype=dtype, delimiter=', ')

        self.assertGreater(len(catData), 0)

        cartoonDir = getPackageDir('sims_photUtils')
        cartoonDir = os.path.join(cartoonDir, 'tests', 'cartoonSedTestData')
        sedDir = getPackageDir('sims_sed_library')

        testBandpasses = {}
        keys = ['u', 'g', 'r', 'i', 'z']

        for kk in keys:
            testBandpasses[kk] = Bandpass()
            testBandpasses[kk].readThroughput(os.path.join(cartoonDir, "test_bandpass_%s.dat" % kk))

        imsimBand = Bandpass()
        imsimBand.imsimBandpass()

        specMap = defaultSpecMap

        ct = 0
        for line in catData:
            bulgeMagList = []
            diskMagList = []
            agnMagList = []
            if line['bulgeName'] == 'None':
                for bp in keys:
                    np.testing.assert_equal(line['%sBulge' % bp], np.NaN)
                    bulgeMagList.append(np.NaN)
            else:
                ct += 1
                dummySed = Sed()
                dummySed.readSED_flambda(os.path.join(sedDir, specMap[line['bulgeName']]))
                fnorm = dummySed.calcFluxNorm(line['bulgeNorm'], imsimBand)
                dummySed.multiplyFluxNorm(fnorm)
                a_int, b_int = dummySed.setupCCM_ab()
                dummySed.addDust(a_int, b_int, A_v=line['bulgeAv'])
                dummySed.redshiftSED(line['redshift'], dimming=True)
                dummySed.resampleSED(wavelen_match=testBandpasses['u'].wavelen)
                for bpName in keys:
                    mag = dummySed.calcMag(testBandpasses[bpName])
                    self.assertAlmostEqual(mag, line['%sBulge' % bpName], 10)
                    bulgeMagList.append(mag)

            if line['diskName'] == 'None':
                for bp in keys:
                    np.assert_equal(line['%sDisk' % bp], np.NaN)
                    diskMagList.append(np.NaN)
            else:
                ct += 1
                dummySed = Sed()
                dummySed.readSED_flambda(os.path.join(sedDir, specMap[line['diskName']]))
                fnorm = dummySed.calcFluxNorm(line['diskNorm'], imsimBand)
                dummySed.multiplyFluxNorm(fnorm)
#.........這裏部分代碼省略.........