import sniacatalogs as snia import numpy as np %matplotlib inline import matplotlib.pyplot as plt import seaborn as sns sns.set() # Give it an ra, dec # sims.photUtils uses O'Donnell 94 to calculate extinction, and does not use any other fits. So, the fit function is not # a parameter at present # Following catsim conventions the entries are in degrees SNObjectfull = snia.snObject.SNObject(ra=30., dec=60.) SNObjectfull.setCoords(ra=30., dec=60.) # The coordinates are stored in radians, but these are not 'public' variables print 'ra, dec: ', SNObjectfull._ra , SNObjectfull._dec # Publicly accessed variables having ra and dec in an array, in the form required by catsim SNObjectfull.skycoord # The color excess has been calculated print 'E(B-V) = ', SNObjectfull.ebvofMW sed = SNObjectfull.SNObjectSED(time=0., wavelen=np.arange(300., 700.,10.)) print sed.wavelen print sed.flambda SN = snia.snObject.SNObject() print 'ra, dec: ', SN._ra , SN._dec # The color excess has been calculated print 'E(B-V) = ', SN.ebvofMW print SN.SNObjectSED(time=0., wavelen=np.arange(300., 800., 10.)) help(SN.set_MWebv) SN.set_MWebv(0.) # in this case, you could have changed this by setting SN.ebvofMW = 0., # but will usually do not allow that sort of thing print SN.SNObjectSED(time=[0.], wavelen=np.arange(300., 800., 10.)) print 'ra, dec: ', SN._ra , SN._dec SN.setCoords(ra=30., dec=60) SN.ebvofMW SN.mwEBVfromMaps() SN.ebvofMW SN.set(z=0.5, x1=0., c =0.) SN.set_source_peakabsmag(absmag=-19.5, band='bessellb', magsys='AB') print SN.summary() from lsst.sims.photUtils import Bandpass from lsst.sims.photUtils.Photometry import PhotometryBase fig, ax = plt.subplots() sed = SN.SNObjectSED(time=0., wavelen=np.arange(300., 800., 1.)) sed2 = SN.SNObjectSED(time=10., wavelen=np.arange(300., 800., 10.)) ax.plot(sed.wavelen, sed.flambda, 'k-', label='t=0.') ax.plot(sed2.wavelen, sed2.flambda, 'r-', label='t=10.') ax.legend(loc='best') ax.set_xlabel(r'$\lambda (nm)$') ax.set_ylabel(r'flambda (ergs/cm^2/s/nm)') # Setup a photUtils class that holds all the LSST bandpasses as a list pbase = PhotometryBase() pbase.loadBandpassesFromFiles() pbase.setupPhiArray_dict() # pbase.phiArray pbase.bandpassDict pbase.loadBandpassesFromFiles() # The following were tried with feature/SIM-1037-only-calculate-requested-magnitudes in photUtils # Store light curve in fluxes l = [] for time in np.arange(-20., 50., 1.): l.append([time]+ SN.catsimBandFluxes(time=time, bandpassobject=pbase.bandpassDict, phiarray=pbase.phiArray, observedBandPassInds=None).tolist()) # Store light curve in mags m = [] for time in np.arange(-20., 50., 1.): m.append([time]+ SN.catsimBandMags(time=time, bandpassobject=pbase.bandpassDict, phiarray=pbase.phiArray).tolist()) fig, ax = plt.subplots() sed = SN.SNObjectSED(time=0., bandpassobject=pbase.bandpassDict) sed2 = SN.SNObjectSED(time=10., wavelen=np.arange(300., 800., 10.)) ax.plot(sed.wavelen, sed.flambda, 'k-', label='t=0.') ax.plot(sed2.wavelen, sed2.flambda, 'r-', label='t=10.') ax.legend(loc='best') ax.set_xlabel(r'$\lambda (nm)$') ax.set_ylabel(r'flambda (ergs/cm^2/s/nm)') sed.flambdaTofnu() sed.manyFluxCalc(phiarray=pbase.phiArray, wavelen_step=pbase.waveLenStep, observedBandPassInd=None) pbase.waveLenStep l = np.asarray(l) fig, ax = plt.subplots() ax.plot(l[:, 0], l[:, 1], label='u') ax.plot(l[:, 0], l[:, 2], label='g') ax.plot(l[:, 0], l[:, 3], label='r') ax.plot(l[:, 0], l[:, 4], label='i') ax.plot(l[:, 0], l[:, 5], label='z') ax.plot(l[:, 0], l[:, 6], label='y') ax.set_xlabel('time -peak') ax.set_ylabel('Flux') ax.legend(loc='best') m = np.asarray(m) fig, ax = plt.subplots() ax.plot(m[:, 0], m[:, 1], 'o', label='u') ax.plot(m[:, 0], m[:, 2], label='g') ax.plot(m[:, 0], m[:, 3], label='r') ax.plot(m[:, 0], m[:, 4], label='i') ax.plot(m[:, 0], m[:, 5], label='z') ax.plot(m[:, 0], m[:, 6], label='y') ax.invert_yaxis() ax.set_xlabel('time -peak') ax.set_ylabel('Magnitude') ax.legend(loc='best') import sncosmo sdssu = sncosmo.get_bandpass('sdssu') print 'u band at given redshift has an effective wavelength of ', sdssu.wave_eff/(1 + SN.get('z')),\ ' which is too close to the lower wavelength limit of the model' min(np.asarray(l)[:,1]) from lsst.sims.catalogs.generation.db import ObservationMetaData from sniacatalogs import sncat from lsst.sims.catalogs.generation.db import CatalogDBObject galDB = CatalogDBObject.from_objid('galaxyTiled') from lsst.sims.catalogs.measures.instance import InstanceCatalog class galCopy(InstanceCatalog): column_outputs = ['id', 'raJ2000', 'decJ2000', 'redshift'] override_formats = {'raJ2000': '%8e', 'decJ2000': '%8e'} mjds = [571203.15] for i, myMJD in enumerate(mjds): myObsMD = ObservationMetaData(boundType='circle', boundLength=0.015, unrefractedRA=5.0, unrefractedDec=15.0, site='LSST', bandpassName=['u', 'g', 'r', 'i', 'z', 'y'], mjd=myMJD) catalog = sncat.SNIaCatalog(db_obj=galDB, obs_metadata=myObsMD) gals = galCopy(db_obj=galDB, obs_metadata=myObsMD) gals.write_catalog('gals.dat') catalog.suppressDimSN catalog.maxTimeSNVisible catalog.suppressDimSN = False # Check print catalog.suppressDimSN print catalog.suppressHighzSN print catalog.maxz catalog.write_catalog('catalog_day.txt', write_metadata=True, meta_key='#@') import numpy as np galcat = np.loadtxt('gals.dat', delimiter=',') sn = np.loadtxt('catalog_day.txt', delimiter=',') %matplotlib inline import matplotlib.pyplot as plt fig, ax = plt.subplots(2, figsize=(6, 8), sharex=True) countssn, bins, patches = ax[0].hist(sn[:, 3], histtype='step', label='SN', lw=2., color='k') countsgals, bins, patches = ax[0].hist(galcat[:, 3], bins=bins, histtype='step', linestyle='dashed', lw=2.0, color='r', label='gals') ax[0].set_ylabel('number of objects') ax[0].legend(loc='upper left') ax[1].plot((bins[:-1] + bins[1:])/2.0, countsgals - countssn, 'o') ax[1].set_xlabel('z') ax[1].set_ylabel('num (SN - gals)') plt.tight_layout() fig.set_label('Comparison of galaxy catalog with SN catalog') catalog.suppressDimSN = True catalog.write_catalog('SN_TimeWindow.txt') !head SN_TimeWindow.txt import lightcurve_utils as lu #This creates a list of obsMetaData with a starting point, a cadence and a total number of epochs def obsMetaDataList(startdate=570180, cadence=3.0, numepochs=20): """ create a list of obsMetaData variables for a list of pointings. Parameters ---------- startdate: float, optional, defaults to 570180 starting date of sequence cadence: float, optional, defaults to 3.0 interval between epochs of observations numepochs: int, optional, defaults to 20 number of epochs of observations Returns ------- list of `~lsst.sims.catalogs.generation.db.ObservationMetaData` variables .. note: This will finally come from OpSims output, this is just a mock up. """ # List of MJDs with a three day cadence with 20 epochs myMJDS = [startdate + cadence*i for i in range(numepochs)] filters = ['u', 'g', 'r', 'i', 'z', 'y'] unrefRA = 5.0 unrefDec = 15.0 boundLen = 0.015 boundType = 'circle' obsMetaDataList = [] for mjd in myMJDS: obsMetaDataList.append(ObservationMetaData(boundType=boundType, unrefractedRA=unrefRA, unrefractedDec=unrefDec, boundLength=boundLen, bandpassName=filters, mjd=mjd)) return obsMetaDataList obsMDList = obsMetaDataList() obsMDList[0].summary obsMDList[1].summary lu.writeCatalogtoDB(dbfile='data/sncat_new.db', dbtable='mysncat', ascii_root='data/SNIa_new_', galdb=galDB, obsMetaDataList=obsMDList) lcs = lu.getLCsFromDB(dbfile='data/sncat_new.db', dbtable='mysncat', lc_root='data/LightCurves/SN_new_') len(lcs) lcs[1][1] fig, ax = plt.subplots(figsize=(24,8)) lu.plotlc(lcs[0]) lu.plotlc(lcs[1], marker='s') ax.invert_yaxis()