Simulation

First, define a set of “observations”. These are the properties of our observations: the time, bandpass and depth.

import sncosmo
from astropy.table import Table
obs = Table({'time': [56176.19, 56188.254, 56207.172],
             'band': ['desg', 'desr', 'desi'],
             'gain': [1., 1., 1.],
             'skynoise': [191.27, 147.62, 160.40],
             'zp': [30., 30., 30.],
             'zpsys':['ab', 'ab', 'ab']})
print obs
skynoise zpsys band gain    time    zp
-------- ----- ---- ---- --------- ----
  191.27    ab desg  1.0  56176.19 30.0
  147.62    ab desr  1.0 56188.254 30.0
   160.4    ab desi  1.0 56207.172 30.0

Suppose we want to simulate a SN with the SALT2 model and the following parameters:

model = sncosmo.Model(source='salt2')
params = {'z': 0.4, 't0': 56200.0, 'x0':1.e-5, 'x1': 0.1, 'c': -0.1}

To get the light curve for this single SN, we’d do:

lcs = sncosmo.realize_lcs(obs, model, [params])
print lcs[0]
   time   band      flux        fluxerr     zp  zpsys
--------- ---- ------------- ------------- ---- -----
 56176.19 desg 96.0531272705  191.27537908 30.0    ab
56188.254 desr 456.360196623  149.22627064 30.0    ab
56207.172 desi  655.40885611 162.579572369 30.0    ab

Note that we’ve passed the function a one-element list, [params], and gotten back a one-element list in return. (The realize_lcs function is designed to operate on lists of SNe for convenience.)

Generating SN parameters

We see above that it is straightforward to simulate SNe once we already know the parameters of each one. But what if we want to pick SN parameters from some defined distribution?

Suppose we want to generate SN parameters for all the SNe we would find in a given search area over a defined period of time. We start by defining an area and time period, as well as a maximum redshift to consider:

area = 1.  # area in square degrees
tmin = 56175.  # minimum time
tmax = 56225.  # maximum time
zmax = 0.7

First, we’d like to get the number and redshifts of all SNe that occur over our 1 square degree and 50 day time period:

redshifts = list(sncosmo.zdist(0., zmax, time=(tmax-tmin), area=area))
print len(redshifts), "SNe"
print "redshifts:", redshifts
9 SNe
redshifts: [0.4199710008856507, 0.3500118339133868, 0.5915676316485601, 0.5857452631151785, 0.49024466410556855, 0.5732679644841575, 0.6224436826380927, 0.5853477892182203, 0.5522300320124105]

Generate a list of SN parameters using these redshifts, drawing x1 and c from normal distributions:

from numpy.random import uniform, normal
params = [{'x0':1.e-5, 'x1':normal(0., 1.), 'c':normal(0., 0.1),
           't0':uniform(tmin, tmax), 'z': z}
          for z in redshifts]
for p in params:
    print p
{'z': 0.4199710008856507, 'x0': 1e-05, 'x1': -0.9739877070754421, 'c': -0.1465835504611458, 't0': 56191.57686616353}
{'z': 0.3500118339133868, 'x0': 1e-05, 'x1': 0.04454878604727126, 'c': -0.04920811869083081, 't0': 56222.76963606611}
{'z': 0.5915676316485601, 'x0': 1e-05, 'x1': -0.26765265677262423, 'c': -0.06456008680932701, 't0': 56211.706219411404}
{'z': 0.5857452631151785, 'x0': 1e-05, 'x1': 0.8255953341731204, 'c': 0.08520083775049729, 't0': 56209.33583211229}
{'z': 0.49024466410556855, 'x0': 1e-05, 'x1': -0.12051827966517584, 'c': -0.09490756669333822, 't0': 56189.37571007927}
{'z': 0.5732679644841575, 'x0': 1e-05, 'x1': 0.3051310078192594, 'c': -0.10967604820261241, 't0': 56198.04368422346}
{'z': 0.6224436826380927, 'x0': 1e-05, 'x1': -0.6329407028587257, 'c': -0.009789183239376284, 't0': 56179.88133113836}
{'z': 0.5853477892182203, 'x0': 1e-05, 'x1': 0.6373371286596669, 'c': 0.05151693090038232, 't0': 56212.04579735962}
{'z': 0.5522300320124105, 'x0': 1e-05, 'x1': 0.04762095339856289, 'c': -0.005018877828783951, 't0': 56182.14827040906}

So far so good. The only problem is that x0 doesn’t vary. We’d like it to be randomly distributed with some scatter around the Hubble line, so it should depend on the redshift. Here’s an alternative:

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

for p in params:
    print p
{'c': -0.060104568346581566, 'x0': 2.9920355958896461e-05, 'z': 0.4199710008856507, 'x1': -0.677121283126299, 't0': 56217.93979718883}
{'c': 0.10405991801014292, 'x0': 2.134500759148091e-05, 'z': 0.3500118339133868, 'x1': 1.6034252041294512, 't0': 56218.008314206476}
{'c': -0.14777109151711296, 'x0': 7.9108889725043354e-06, 'z': 0.5915676316485601, 'x1': -2.2082282760850993, 't0': 56218.013686428785}
{'c': 0.056034777154805086, 'x0': 6.6457371815973038e-06, 'z': 0.5857452631151785, 'x1': 0.675413080007434, 't0': 56189.03517395757}
{'c': -0.0709158052635228, 'x0': 1.2228145655148946e-05, 'z': 0.49024466410556855, 'x1': 0.5449847454420981, 't0': 56198.02895700289}
{'c': -0.22101146234021096, 'x0': 7.4299221264917702e-06, 'z': 0.5732679644841575, 'x1': -1.543245858395605, 't0': 56189.04585414441}
{'c': 0.06964843664572477, 'x0': 9.7121906557832662e-06, 'z': 0.6224436826380927, 'x1': 1.7419604610283943, 't0': 56212.827270197355}
{'c': 0.07320513053870191, 'x0': 3.22205341646521e-06, 'z': 0.5853477892182203, 'x1': -0.39658066375434153, 't0': 56200.421464066916}
{'c': 0.18555773972769227, 'x0': 7.5955258508017471e-06, 'z': 0.5522300320124105, 'x1': -0.24463691193386283, 't0': 56190.492271332616}

Now we can generate the lightcurves for these parameters:

lcs = sncosmo.realize_lcs(obs, model, params)
print lcs[0]
   time   band      flux       fluxerr     zp  zpsys
--------- ---- ------------- ------------ ---- -----
 56176.19 desg 6.70520005464       191.27 30.0    ab
56188.254 desr 106.739113709       147.62 30.0    ab
56207.172 desi  1489.7521011 164.62420476 30.0    ab

Note that the “true” parameters are saved in the metadata of each SN:

lcs[0].meta
{'c': -0.060104568346581566,
 't0': 56217.93979718883,
 'x0': 2.9920355958896461e-05,
 'x1': -0.677121283126299,
 'z': 0.4199710008856507}