Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/eop.py: 8%
299 statements
« prev ^ index » next coverage.py v7.6.4, created at 2024-10-31 18:48 +0000
1import math
2import re
3import shutil
5from casatasks import casalog
6from casatools import ctsys
7from casatools import table
8from casatools import msmetadata
10import numpy as np
11import erfa
13c = 299792458.0
14arcsec = 0.00000484813681109535993589914098765
16def ut1_utc(eops, jd1, jd2):
17 mjd1, mjd2 = jd1 - 2400000.5, jd2
18 i = np.searchsorted(eops['MJD'], mjd1 + mjd2, side='right')
19 i1 = np.clip(i, 1, len(eops) - 1)
20 i0 = i1 - 1
21 mjd_0, mjd_1 = eops['MJD'][i0], eops['MJD'][i1]
22 val_0, val_1 = eops['UT1_UTC'][i0], eops['UT1_UTC'][i1]
23 d_val = val_1 - val_0
24 d_val -= round(d_val)
25 val = val_0 + (mjd1 - mjd_0 + mjd2) / (mjd_1 - mjd_0) * d_val
26 return val
28def pm_xy(eops, jd1, jd2):
29 mjd1, mjd2 = jd1 - 2400000.5, jd2
30 i = np.searchsorted(eops['MJD'], mjd1 + mjd2, side='right')
31 i1 = np.clip(i, 1, len(eops) - 1)
32 i0 = i1 - 1
33 mjd_0, mjd_1 = eops['MJD'][i0], eops['MJD'][i1]
34 val_0, val_1 = eops['PM_x'][i0], eops['PM_x'][i1]
35 xp = val_0 + (mjd1 - mjd_0 + mjd2) / (mjd_1 - mjd_0) * (val_1 - val_0)
36 val_0, val_1 = eops['PM_y'][i0], eops['PM_y'][i1]
37 yp = val_0 + (mjd1 - mjd_0 + mjd2) / (mjd_1 - mjd_0) * (val_1 - val_0)
38 return xp * arcsec, yp * arcsec
40def utctt(jd1, jd2):
41 tai = erfa.utctai(jd1, jd2)
42 return erfa.taitt(tai[0], tai[1])
44def calc_delay(eops, pos, p, mjd):
45 xp, yp = pm_xy(eops, 2400000.5, mjd)
46 dut1 = ut1_utc(eops, 2400000.5, mjd) / 86400
47 tt = utctt(2400000.5, mjd)
48 r = erfa.c2t00a(tt[0], tt[1], 2400000.5 + mjd, dut1, xp, yp)
49 p = erfa.rxp(r, p)
50 return erfa.pdp(pos, p) / c
52def add_row(tb, obsid, field, antenna, spws, reffreq, t, pos, pc,
53 old_eops, new_eops, delta):
54 dtm = delta / 86400
55 tm = t / 86400
56 old_delay = [0, 0, 0]
57 old_delay[0] = calc_delay(old_eops, pos, pc, tm - dtm)
58 old_delay[1] = calc_delay(old_eops, pos, pc, tm)
59 old_delay[2] = calc_delay(old_eops, pos, pc, tm + dtm)
60 new_delay = [0, 0, 0]
61 new_delay[0] = calc_delay(new_eops, pos, pc, tm - dtm)
62 new_delay[1] = calc_delay(new_eops, pos, pc, tm)
63 new_delay[2] = calc_delay(new_eops, pos, pc, tm + dtm)
64 delay = new_delay[1] - old_delay[1]
65 new_rate = (new_delay[2] - new_delay[0]) / (2 * delta)
66 old_rate = (old_delay[2] - old_delay[0]) / (2 * delta)
67 rate = new_rate - old_rate
69 param = np.zeros(shape=(8,1), dtype='float32')
70 paramerr = -np.ones(shape=(8,1), dtype='float32')
71 flag = np.zeros(shape=(8,1), dtype='float32')
72 snr = np.ones(shape=(8,1), dtype='float32')
73 weight = np.ones(shape=(8,1), dtype='float32')
75 param[1] = delay * 1e9
76 param[2] = rate
77 param[5] = delay * 1e9
78 param[6] = rate
80 for spw in spws:
81 param[0] = 2 * math.pi * delay * reffreq[spw]
82 param[4] = 2 * math.pi * delay * reffreq[spw]
84 row = tb.nrows()
85 tb.addrows(1)
86 tb.putcell('TIME', row, t)
87 tb.putcell('INTERVAL', row, 0)
88 tb.putcell('ANTENNA1', row, antenna)
89 tb.putcell('ANTENNA2', row, -1)
90 tb.putcell('FIELD_ID', row, field)
91 tb.putcell('SCAN_NUMBER', row, -1)
92 tb.putcell('OBSERVATION_ID', row, obsid)
93 tb.putcell('SPECTRAL_WINDOW_ID', row, spw)
94 tb.putcell('FPARAM', row, param)
95 tb.putcell('PARAMERR', row, paramerr)
96 tb.putcell('FLAG', row, flag)
97 tb.putcell('SNR', row, snr)
98 tb.putcell('WEIGHT', row, weight)
99 continue
100 pass
102def get_eops_from_ms(vis, obsid):
103 tb = table()
105 eops = {}
106 eops['MJD'] = []
107 eops['UT1_UTC'] = []
108 eops['PM_x'] = []
109 eops['PM_y'] = []
111 try:
112 tb.open(vis + '/EARTH_ORIENTATION')
113 res = tb.query('OBSERVATION_ID == %d' % obsid)
114 for row in range(res.nrows()):
115 pm = res.getcell('PM', row)
116 eops['MJD'].append(res.getcell('TIME', row) / 86400)
117 eops['PM_x'].append(pm[0] / arcsec)
118 eops['PM_y'].append(pm[1] / arcsec)
119 eops['UT1_UTC'].append(res.getcell('UT1_UTC', row))
120 continue
121 tb.close()
122 except:
123 return None
125 if len(eops['MJD']) == 0:
126 return None
127 return eops
129def get_eops_from_casadata(mjd_min, mjd_max):
130 tb = table()
132 eops = {}
133 eops['MJD'] = []
134 eops['UT1_UTC'] = []
135 eops['PM_x'] = []
136 eops['PM_y'] = []
138 try:
139 path = ctsys.resolve('geodetic/IERSeop2000')
140 tb.open(path)
141 res = tb.query('MJD > %f && MJD < %f' % (mjd_min - 0.5, mjd_max + 0.5))
142 for row in range(res.nrows()):
143 eops['MJD'].append(res.getcell('MJD', row))
144 eops['UT1_UTC'].append(res.getcell('dUT1', row))
145 eops['PM_x'].append(res.getcell('x', row))
146 eops['PM_y'].append(res.getcell('y', row))
147 continue
148 tb.close()
149 except:
150 return None
152 if len(eops['MJD']) == 0:
153 return None
154 return eops
156# Parse UNSO finals version 2.0 format
157def parse_usno_2_0(fp, mjd_min, mjd_max):
158 eops = {}
159 eops['MJD'] = []
160 eops['UT1_UTC'] = []
161 eops['PM_x'] = []
162 eops['PM_y'] = []
164 for line in fp:
165 if re.match(r'#', line):
166 continue
167 jd_tai = float(line[0:9])
168 mjd = jd_tai - 2400000.5
169 if mjd < mjd_min:
170 continue
171 if mjd > mjd_max:
172 continue
173 pm_x = float(line[10:17]) / 10
174 pm_y = float(line[18:25]) / 10
175 ut1_tai = float(line[26:35]) / 1e6
176 jd_utc = erfa.taiutc(jd_tai, 0)
177 utc_tai = (jd_utc[1] + (jd_utc[0] - jd_tai)) * 86400
178 ut1_utc = ut1_tai - utc_tai
179 mjd = jd_utc[1] + (jd_utc[0] - 2400000.5)
180 eops['MJD'].append(mjd)
181 eops['UT1_UTC'].append(ut1_utc)
182 eops['PM_x'].append(pm_x)
183 eops['PM_y'].append(pm_y)
184 continue
186 if len(eops['MJD']) == 0:
187 return None
188 return eops
190# Parse UNSO finals version 2.1 format
191def parse_usno_2_1(fp, mjd_min, mjd_max):
192 eops = {}
193 eops['MJD'] = []
194 eops['UT1_UTC'] = []
195 eops['PM_x'] = []
196 eops['PM_y'] = []
198 for line in fp:
199 if re.match(r'#', line):
200 continue
201 jd_tai = float(line[0:9])
202 mjd = jd_tai - 2400000.5
203 if mjd < mjd_min:
204 continue
205 if mjd > mjd_max:
206 continue
207 pm_x = float(line[10:18]) / 10
208 pm_y = float(line[19:27]) / 10
209 ut1_tai = float(line[28:39]) / 1e6
210 jd_utc = erfa.taiutc(jd_tai, 0)
211 utc_tai = (jd_utc[1] + (jd_utc[0] - jd_tai)) * 86400
212 ut1_utc = ut1_tai - utc_tai
213 mjd = jd_utc[1] + (jd_utc[0] - 2400000.5)
214 eops['MJD'].append(mjd)
215 eops['UT1_UTC'].append(ut1_utc)
216 eops['PM_x'].append(pm_x)
217 eops['PM_y'].append(pm_y)
218 continue
220 if len(eops['MJD']) == 0:
221 return None
222 return eops
224# Parse old IERS format
225def parse_eopc04(fp, mjd_min, mjd_max):
226 eops = {}
227 eops['MJD'] = []
228 eops['UT1_UTC'] = []
229 eops['PM_x'] = []
230 eops['PM_y'] = []
232 for line in fp:
233 if re.match(r'\s', line) or re.match(r'#', line):
234 continue
235 mjd = int(line[12:19])
236 if mjd < mjd_min:
237 continue
238 if mjd > mjd_max:
239 continue
240 pm_x = float(line[20:30])
241 pm_y = float(line[31:41])
242 ut1_utc = float(line[42:53])
243 eops['MJD'].append(mjd)
244 eops['UT1_UTC'].append(ut1_utc)
245 eops['PM_x'].append(pm_x)
246 eops['PM_y'].append(pm_y)
247 continue
249 if len(eops['MJD']) == 0:
250 return None
251 return eops
253def get_eops_from_file(infile, mjd_min, mjd_max):
254 fp = open(infile)
255 if not fp:
256 msg = 'Cannot open ' + infile
257 casalog.post(msg, 'SEVERE')
258 return None
259 try:
260 line = fp.readline()
261 if line.startswith("EOP-MOD Ver 2.0"):
262 return parse_usno_2_0(fp, mjd_min, mjd_max)
263 elif line.startswith("EOP-MOD Ver 2.1"):
264 return parse_usno_2_1(fp, mjd_min, mjd_max)
265 else:
266 return parse_eopc04(fp, mjd_min, mjd_max)
267 except:
268 return None
270def do_generate_eop(vis, caltable, infile):
271 msmd = msmetadata()
272 tb = table()
274 # Use a 1 minute calibration interval. This should be sufficient
275 # as we'll use both the delay and delay rate in the interpolation.
276 delta = 60
278 tb.open(caltable + '/SPECTRAL_WINDOW')
279 reffreq = tb.getcol('CHAN_FREQ')[0]
280 tb.close()
282 msmd.open(vis)
283 tb.open(caltable, nomodify=False)
285 # Prefetch antenna positions.
286 positions = {}
287 for antenna in msmd.antennaids():
288 pos = msmd.antennaposition(antenna)
289 pos = erfa.s2p(pos['m0']['value'], pos['m1']['value'], pos['m2']['value'])
290 positions[antenna] = pos
291 continue
293 # Prefetch phase centers.
294 phasecenters = {}
295 for field in msmd.fieldsforname():
296 pc = msmd.phasecenter(field)
297 pc = erfa.s2c(pc['m0']['value'], pc['m1']['value'])
298 phasecenters[field] = pc
299 continue
301 # The original EOPs may include differen values for the same day
302 # for different observations. Therefore we iterate over
303 # observations and create calibration table entries that include
304 # the OBSERVATION_ID.
305 for obsid in range(msmd.nobservations()):
306 # Get the original EOPs for this obervation from the MS.
307 old_eops = get_eops_from_ms(vis, obsid)
308 if not old_eops:
309 msg = 'No EOPS found in ' + vis
310 casalog.post(msg, 'SEVERE')
311 continue
313 # Get the updated EOPs for the timerange covered by the
314 # original EOPs.
315 mjd_min = np.min(old_eops['MJD'])
316 mjd_max = np.max(old_eops['MJD'])
317 if infile and infile != 'None':
318 new_eops = get_eops_from_file(infile, mjd_min, mjd_max)
319 else:
320 new_eops = get_eops_from_casadata(mjd_min, mjd_max)
321 if not new_eops:
322 msg = 'No updated EOPS found for MJD ' + str(mjd_min) + '-' \
323 + str(mjd_max) + ' in ' + infile
324 casalog.post(msg, 'SEVERE')
325 continue
327 # Iterate by scan over all the antennas and fields in this
328 # observation.
329 for scan in msmd.scannumbers(obsid=obsid):
330 for antenna in msmd.antennasforscan(scan=scan, obsid=obsid):
331 for field in msmd.fieldsforscan(scan=scan, obsid=obsid):
332 pos = positions[antenna]
333 pc = phasecenters[field]
335 spws = msmd.spwsforscan(scan=scan, obsid=obsid)
336 times = msmd.timesforscan(scan=scan, obsid=obsid)
337 t0 = 0
338 for t in times:
339 if t - t0 > delta:
340 t0 = t
341 add_row(tb, obsid, field, antenna, spws, reffreq,
342 t, pos, pc, old_eops, new_eops, delta)
343 pass
344 continue
345 if t != t0:
346 add_row(tb, obsid, field, antenna, spws, reffreq,
347 t, pos, pc, old_eops, new_eops, delta)
348 pass
350 continue
351 continue
352 continue
353 continue
355 nrows = tb.nrows()
356 tb.close()
358 return nrows
360def generate_eop(vis, caltable, infile):
361 nrows = 0
362 try:
363 nrows = do_generate_eop(vis, caltable, infile)
364 finally:
365 if nrows == 0:
366 shutil.rmtree(caltable)