In [1]:
import numpy as np
%matplotlib inline
import matplotlib.pyplot as plt
import sncosmo
import os,sys
from astropy.units import Unit
from numpy.random import normal, uniform

The main structure in SNCosmo is a model (as seen from the earth) described as a time series of spectra.


In [48]:
model  =  sncosmo.Model(source="Salt2-extended")

The SALT model has parameters (x1 (related to 'stretch', ie. temporal width of light curve' , c (related slope of spectra against wavelength), x0 (close to apparent magnitude), t0 (time of peak), obviously redshift z). The model parameters may be set through a dictionary as below, or individually


In [3]:
redshifts = [0.3, 0.5, 0.8, 1.1] 
params = []
for z in redshifts:
    #mabs = normal(-19.3, 0.3)
    mabs = -19.3
    model.set(z=z)
    model.set_source_peakabsmag(mabs, 'bessellb', 'ab')
    x0 = model.get('x0')
    p = {'z':z, 't0':0, 'x0':x0, 'x1': 0., 'c':0.}
    #p = {'z':z, 't0':uniform(tmin, tmax), 'x0':x0, 'x1': normal(0., 1.), 'c': normal(0., 0.1)}
    params.append(p)

#params is a list of dictionaries. Each dictionary stores the SALT2 model parameters of a SN. 
for p in params:
    print p


{'c': 0.0, 'x0': 3.3789156007734127e-05, 'z': 0.3, 'x1': 0.0, 't0': 0}
{'c': 0.0, 'x0': 1.0115257191922824e-05, 'z': 0.5, 'x1': 0.0, 't0': 0}
{'c': 0.0, 'x0': 3.2102222568952169e-06, 'z': 0.8, 'x1': 0.0, 't0': 0}
{'c': 0.0, 'x0': 1.457503209988081e-06, 'z': 1.1, 'x1': 0.0, 't0': 0}

Let us take a grid of points from -20 to 70 days on the earth. Note that the extremes are too dim for a SNIa to be observed with decent signal in usual surveys, and so the SALT model is only provided at -20, 50 in the rest frame; even at the extremeties of this, it is quite unconstrained since the number of SNIa epochs used in constraining the model is small.


In [4]:
days =np.arange(-20, 50,1.)
wavelens = np.arange(3000.,8000.,10.)
times = [-20., 0., 20., 40.]

In [49]:
#Setting all the model parameters, can also be done individually as model.set(z= 0.5, c=0.1, x1=1.0)

In [15]:
model.set(**params[0])

In [14]:
data = model.flux(time=times, wave=wavelens)
colors = ['b','k','g','r']
labels = ['-20', '0', '20', '40']
fig, ax = plt.subplots()
for i,dat in enumerate(data):
    ax.plot(wavelens, dat, color=colors[i], label=labels[i])
    ax.legend(loc='best')
ax.set_ylabel("flux")
ax.set_xlabel('wavelength (Ang)')


Out[14]:
<matplotlib.text.Text at 0x10a3d2450>

In [19]:
#magnitudes in sdss bands
g = model.bandmag(band='sdssg', magsys='ab', time=days)
r = model.bandmag(band='sdssr', magsys='ab', time=days)
i = model.bandmag(band='sdssi', magsys='ab', time=days)
z = model.bandmag(band='sdssz', magsys='ab', time=days)

In [29]:
#magnitudes for model with higher c
model.set(c=0.1)
ghc = model.bandmag(band='sdssg', magsys='ab', time=days)
rhc = model.bandmag(band='sdssr', magsys='ab', time=days)
ihc = model.bandmag(band='sdssi', magsys='ab', time=days)
zhc = model.bandmag(band='sdssz', magsys='ab', time=days)

In [37]:
#magnitudes for model with lower c
model.set(c=-0.1)
glc = model.bandmag(band='sdssg', magsys='ab', time=days)
rlc = model.bandmag(band='sdssr', magsys='ab', time=days)
ilc = model.bandmag(band='sdssi', magsys='ab', time=days)
zlc = model.bandmag(band='sdssz', magsys='ab', time=days)

In [50]:
fig, ax = plt.subplots(3,1, figsize=(4,10))
ax[0].plot(i, g - r, lw=2, color='k')
ax[1].plot(i, r - i, lw=2, color='k')
ax[2].plot(i, z - r, lw=2, color='k', label='c=0.')
ax[0].set_ylabel('g-r')
ax[1].set_ylabel('r-i')
ax[2].set_ylabel('i-z')
ax[2].set_xlabel('i')
#plot values with c = 0.1
ax[0].plot(i, ghc - rhc, lw=2, color= 'b')
ax[1].plot(i, rhc - ihc, lw=2, color='b')
ax[2].plot(i, zhc - rhc, lw=2, color='b', label='c=0.1')
#plot values with c = -0.1
ax[0].plot(i, glc - rlc, lw=2, color= 'r')
ax[1].plot(i, rlc - ilc, lw=2, color='r')
ax[2].plot(i, zlc - rlc, lw=2, color='r', label='c=-0.1')


Out[50]:
[<matplotlib.lines.Line2D at 0x10c89b410>]

In [44]:
fig2, ax2 = plt.subplots(1,3, figsize=(12,3))
ax2[0].plot(days, g- r, lw=2.,color='k')
ax2[1].plot(days, r -i, lw=2., color='k')
ax2[2].plot(days, z - r, lw=2, color='k')
#plot values with c = 0.1
ax2[0].plot(days, ghc - rhc, lw=2.,color='b')
ax2[1].plot(days, rhc - ihc, lw=2., color='b')
ax2[2].plot(days, zlc - rhc, lw=2, color='b')
#plot values with c = -0.1
ax2[0].plot(days, glc - rlc, lw=2.,color='r')
ax2[1].plot(days, rlc - ilc, lw=2., color='r')
ax2[2].plot(days, zlc - rlc, lw=2, color='r')


Out[44]:
[<matplotlib.lines.Line2D at 0x10c5cdbd0>]

In [45]:
arr = np.array([days, g,r,i,z,ghc, rhc, ihc, zhc, glc, rlc, ilc, zlc]).T

In [46]:
np.savetxt('mags.dat', arr)