1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
180
181
|
Simulation
==========
First, define a set of "observations". These are the properties of our
observations: the time, bandpass and depth.
.. code:: python
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
.. parsed-literal::
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:
.. code:: python
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:
.. code:: python
lcs = sncosmo.realize_lcs(obs, model, [params])
print lcs[0]
.. parsed-literal::
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:
.. code:: python
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:
.. code:: python
redshifts = list(sncosmo.zdist(0., zmax, time=(tmax-tmin), area=area))
print len(redshifts), "SNe"
print "redshifts:", redshifts
.. parsed-literal::
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:
.. code:: python
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
.. parsed-literal::
{'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:
.. code:: python
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
.. parsed-literal::
{'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:
.. code:: python
lcs = sncosmo.realize_lcs(obs, model, params)
print lcs[0]
.. parsed-literal::
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:
.. code:: python
lcs[0].meta
.. parsed-literal::
{'c': -0.060104568346581566,
't0': 56217.93979718883,
'x0': 2.9920355958896461e-05,
'x1': -0.677121283126299,
'z': 0.4199710008856507}
|
