-
Notifications
You must be signed in to change notification settings - Fork 0
/
ecsv_to_wstd.py
159 lines (120 loc) · 5.15 KB
/
ecsv_to_wstd.py
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
#!/usr/bin/env python
from datetime import datetime
import os
import sys
import numpy as np
from astropy.coordinates import SkyCoord
from astropy.table import Table
import astropy.units as u
from Wstd import Wstd
whirc_filter_lookup = {
'J': 'WHIRCJ',
'H': 'WHIRCH',
'KS': 'WHIRCK',
}
def get_redshift_for_target(target, target_info_file='observed_target_info.dr1.txt'):
"""Return redshift if known. Else return None."""
target_info = Table.read(target_info_file, format='ascii.commented_header')
w, = np.where(target_info['Name'] == target)
if len(w) < 1:
return None
target_info = target_info[w]
redshift = target_info['z']
return redshift
def build_header(table, target):
"""Return a Wstd header from information in table"""
redshift = get_redshift_for_target(target)
# Assume the coordinates of first row are correct:
col_ra, col_dec = table['coord_ra'], table['coord_dec']
coord = SkyCoord(col_ra[0], col_dec[0], unit=(col_ra.unit, col_dec.unit))
utc = datetime.utcnow()
time_stamp_str = "# Lightcurve generated at {} UTC".format(utc.isoformat())
header_str = Wstd().formatWstdHeader(
event=target,
RA=coord.ra.to_string(unit=u.hour, sep=':', pad=True),
Dec=coord.dec.to_string(unit=u.deg, sep=':', pad=True, alwayssign=True),
MW_E_BmV=0.,
z_hel=redshift,
lII=coord.galactic.l.to_string(unit=u.deg, decimal=True, precision=6),
bII=coord.galactic.b.to_string(unit=u.deg, decimal=True, precision=6),
extraHeaderInfo=time_stamp_str)
return header_str
class convert_ab_vega:
"""Convert AB magnitude and fluxes to Vega Magnitude and fluxes.
Supported filters are 'J', 'H', 'K'
http://www.astronomy.ohio-state.edu/~martini/usefuldata.html
Blanton and Roweis, 2007, AJ, 133, 734.
http://adsabs.harvard.edu/abs/2007AJ....133..734B
>>> convert = convert_ab_vega()
>>> ab_mag, vega_mag = 10, 9.09
>>> print('%5.2f' % convert.ab_to_vega_mag(10, 'J'))
9.09
>>> print('%5.2f' % convert.vega_to_ab_mag(10, 'H'))
11.39
>>> ab_flux, vega_flux = 10**(-0.4*10), 0.00023120648
>>> print('%6.2g' % convert.ab_to_vega_flux(ab_flux, 'J'))
0.00023
>>> print('%6.8g' % convert.vega_to_ab_flux(vega_flux, 'J'))
0.0001
"""
def __init__(self):
self.mAB_minus_mVega = {'J': 0.91, 'H': 1.39, 'K': 1.85}
def ab_to_vega_flux(self, ab_flux, filt):
"""Convert Vega flux to AB flux for given filter."""
return ab_flux * 10**(+0.4*self.mAB_minus_mVega[filt])
def vega_to_ab_flux(self, vega_flux, filt):
"""Convert AB flux to Vega flux for given filter."""
return vega_flux * 10**(-0.4*self.mAB_minus_mVega[filt])
def ab_to_vega_mag(self, ab_mag, filt):
"""Convert AB magnitude to Vega Magnitude for given filter.
"""
return ab_mag - self.mAB_minus_mVega[filt]
def vega_to_ab_mag(self, vega_mag, filt):
"""Convert Vega magnitude to AB Magnitude for given filter."""
return vega_mag + self.mAB_minus_mVega[filt]
def build_data(table, in_system='AB', out_system='Vega'):
"""Return a Wstd data str from information in table
table is expected to be calibrated in 'in_system' (default 'AB')
Output will be in 'out_system' (default 'Vega')
"""
filt = table['filter']
mjd = table['mjd']
flux = table['base_PsfFlux_flux_zp25']
flux_err = table['base_PsfFlux_fluxSigma_zp25']
if in_system == 'AB' and out_system == 'Vega':
convert = convert_ab_vega()
flux = [convert.ab_to_vega_flux(fl, fi) for fl, fi in zip(flux, filt)]
flux_err = [convert.ab_to_vega_flux(fle, fi) for fle, fi in zip(flux_err, filt)]
if in_system == 'Vega' and out_system == 'AB':
convert = convert_ab_vega()
flux = [convert.vega_to_ab_flux(fl, fi) for fl, fi in zip(flux, filt)]
flux_err = [convert.vega_to_ab_flux(fle, fi) for fle, fi in zip(flux_err, filt)]
filters = [whirc_filter_lookup[f] for f in filt]
observations = ['%s_%011.5f' % (out_system, m) for m in mjd]
dataStr = Wstd().formatWstdLC(observations, mjd, filters, flux, flux_err)
return dataStr
def write_header_data(header, data, file):
"""Write header, data to file."""
with open(file, 'w') as f:
f.write(header)
f.write(data)
def convert_ecsv_to_wstd(infile, outfile=None, in_system='AB', out_system='Vega'):
"""
Convert an ECSV file generated by sn_lightcurve.py to a Wstd format
"""
table = Table.read(infile, format='ascii.ecsv')
target = os.path.splitext(infile)[0]
field, target = str.split(target, '_')
if outfile is None:
outfile = '%s_%s.%s' % (field, target, Wstd().WstdLCsuffix)
header = build_header(table, target)
data = build_data(table, in_system=in_system, out_system=out_system)
write_header_data(header, data, outfile)
if __name__ == "__main__":
verbose = True
in_system, out_system = 'AB', 'Vega'
files_to_convert = sys.argv[1:]
for infile in files_to_convert:
if verbose:
print("Converting %s" % infile)
convert_ecsv_to_wstd(infile, in_system=in_system, out_system=out_system)