Coverage for /wheeldirectory/casa-6.7.0-12-py3.10.el8/lib/py/lib/python3.10/site-packages/casatasks/private/eop.py: 85%

299 statements  

« prev     ^ index     » next       coverage.py v7.6.4, created at 2024-11-01 07:19 +0000

1import math 

2import re 

3import shutil 

4 

5from casatasks import casalog 

6from casatools import ctsys 

7from casatools import table 

8from casatools import msmetadata 

9 

10import numpy as np 

11import erfa 

12 

13c = 299792458.0 

14arcsec = 0.00000484813681109535993589914098765 

15 

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 

27 

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 

39 

40def utctt(jd1, jd2): 

41 tai = erfa.utctai(jd1, jd2) 

42 return erfa.taitt(tai[0], tai[1]) 

43 

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 

51 

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 

68 

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') 

74 

75 param[1] = delay * 1e9 

76 param[2] = rate 

77 param[5] = delay * 1e9 

78 param[6] = rate 

79 

80 for spw in spws: 

81 param[0] = 2 * math.pi * delay * reffreq[spw] 

82 param[4] = 2 * math.pi * delay * reffreq[spw] 

83 

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 

101 

102def get_eops_from_ms(vis, obsid): 

103 tb = table() 

104 

105 eops = {} 

106 eops['MJD'] = [] 

107 eops['UT1_UTC'] = [] 

108 eops['PM_x'] = [] 

109 eops['PM_y'] = [] 

110 

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 

124 

125 if len(eops['MJD']) == 0: 

126 return None 

127 return eops 

128 

129def get_eops_from_casadata(mjd_min, mjd_max): 

130 tb = table() 

131 

132 eops = {} 

133 eops['MJD'] = [] 

134 eops['UT1_UTC'] = [] 

135 eops['PM_x'] = [] 

136 eops['PM_y'] = [] 

137 

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 

151 

152 if len(eops['MJD']) == 0: 

153 return None 

154 return eops 

155 

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'] = [] 

163 

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 

185 

186 if len(eops['MJD']) == 0: 

187 return None 

188 return eops 

189 

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'] = [] 

197 

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 

219 

220 if len(eops['MJD']) == 0: 

221 return None 

222 return eops 

223 

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'] = [] 

231 

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 

248 

249 if len(eops['MJD']) == 0: 

250 return None 

251 return eops 

252 

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 

269 

270def do_generate_eop(vis, caltable, infile): 

271 msmd = msmetadata() 

272 tb = table() 

273 

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 

277 

278 tb.open(caltable + '/SPECTRAL_WINDOW') 

279 reffreq = tb.getcol('CHAN_FREQ')[0] 

280 tb.close() 

281 

282 msmd.open(vis) 

283 tb.open(caltable, nomodify=False) 

284 

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 

292 

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 

300 

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 

312 

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 

326 

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] 

334 

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 

349 

350 continue 

351 continue 

352 continue 

353 continue 

354 

355 nrows = tb.nrows() 

356 tb.close() 

357 

358 return nrows 

359 

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)