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

428 statements  

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

1''' 

2 main module for retreiving JPL-Horizons ephemeris data and 

3 convert it to CASA readable format (i.e. CASA table ) 

4''' 

5import os 

6import time 

7from math import sqrt, sin, cos, log 

8import numpy as np 

9 

10from casatools import table, quanta, measures 

11from casatasks import casalog 

12 

13_tb = table() 

14_qa = quanta() 

15_me = measures() 

16 

17# debug 

18_debug = False 

19 

20def gethorizonsephem(objectname, starttime, stoptime, incr, outtable, asis=False, rawdatafile=''): 

21 """ 

22 Main driver function for ephemeris data query from JPL-Horizons 

23  

24 arguments: 

25 objectname: ephemeris object name (case insensitive). Try to convert  

26 a common name to ID if asis = False. 

27 starttime: start time of the ephemeris data (expects YYYY/MM/DD/HH:MM  

28 format 

29 stoptime: stop (end) time of the ephemeris data in the same format as 

30 starttime 

31 incr: time increment (interval) of the ephemeris data 

32 outtable: output CASA table name 

33 asis: a toggle for additinal check to resolve object name  

34 rawdatafile: raw query result file name (optional) 

35 """ 

36 

37 # commented ones are not currently supported in setjy 

38 # Juno's table is exist in the data repo but not supported in setjy 

39 asteroids = {'ceres': '1;', 

40 'pallas': '2;', 

41 'juno': '3;', # large crater and temperature varies 

42 'vesta': '4;', 

43 'lutetia': '21;' 

44 } 

45 planets_and_moons = {'sun': '10', 

46 'mercury': '199', 

47 'venus': '299', 

48 'moon': '301', 

49 'mars': '499', 

50 'jupiter': '599', 

51 'io': '501', 

52 'europa': '502', 

53 'ganymede': '503', 

54 'callisto': '504', 

55 'saturn': '699', 

56 'titan': '606', 

57 'uranus': '799', 

58 'neptune': '899', 

59 'pluto': '999'} 

60 known_objects = planets_and_moons 

61 known_objects.update(asteroids) 

62 # default quantities (required by setjy) for CASA 

63 quantities = '1,14,15,17,19,20,24' 

64 ang_format = 'DEG' 

65 

66 # check objectname. If it matches with the known object list. If asis is specified 

67 # it will skip the check. 

68 # For small bodies such as comets, objectname must include 'DES=' following designated 

69 # object name or id (SPK ID). Based on the JPL-Horizons documentaiton, for IAU ID, 'DES=' 

70 # can be omitted. https://ssd.jpl.nasa.gov/horizons/app.html#command 

71 start_time = None 

72 stop_time = None 

73 step_size = None 

74 if not asis: 

75 if not objectname.lower() in known_objects: 

76 raise ValueError( 

77 f"{objectname} is not in the known object list for CASA. To skip this check set asis=True") 

78 else: 

79 target = known_objects[objectname.lower()] 

80 else: 

81 target = objectname 

82 try: 

83 if starttime.startswith('JD') or starttime.startswith('MJD'): 

84 start_time = starttime 

85 stop_time = stoptime 

86 else: 

87 start_time = _qa.time(starttime, form='ymd') 

88 stop_time = _qa.time(stoptime, form='ymd') 

89 except ValueError as e: 

90 casalog.post(e) 

91 

92 try: 

93 step_size = incr.replace(' ', '') 

94 except ValueError as e: 

95 casalog.post(e) 

96 

97 if _debug: 

98 print("target=",target) 

99 print(f"start_time={start_time}, stop_time={stop_time}, step_size={step_size}") 

100 print(f"quantities={quantities}, ang_format={ang_format}") 

101 ephemdata = queryhorizons(target, start_time, stop_time, step_size, quantities, ang_format, rawdatafile) 

102 # return ephemdata 

103 if ephemdata and 'result' in ephemdata: 

104 casalog.post('converting ephemeris data to a CASA table') 

105 tocasatb(ephemdata, outtable) 

106 if not os.path.exists(outtable): 

107 casalog.post('Failed to produce the CASA ephemeris table', 'WARN') 

108 if rawdatafile=='': 

109 casalog.post('You may want to re-run the task with rawdatafile'+\ 

110 ' parameter specified and check the content of the output raw data file for the further information on the error.','WARN') 

111 else: 

112 casalog.post(f'Please check the content of {rawdatafile} for the further information of the error','WARN') 

113 

114 

115def queryhorizons(target, starttime, stoptime, stepsize, quantities, ang_format, rawdatafile=''): 

116 """ 

117 Submit a query to the JPL-Horizons API 

118  

119 arguments: 

120 target: a target solar system object name/id (str) 

121 starttime: ephemeris start time (e.g. '2021/12/01/00:00') 

122 stoptime: ephemeris stop time 

123 stepsize: ephemeris data time increment 

124 quantities: output data quantities given as a list of integers which represent the specific  

125 data as defined in https://ssd,jpl.nasa.gov/horizons/manual.html 

126 ang_format: RA, Dec output format ('DEG' or 'HMS')  

127 rawdatafile: a file name to save raw query results (Optional) 

128 

129 """ 

130 import ast 

131 import certifi 

132 import ssl 

133 from urllib.request import Request, urlopen 

134 from urllib.parse import urlencode 

135 from urllib.error import URLError 

136 import socket 

137 

138 context = ssl.create_default_context(cafile=certifi.where()) 

139 urlbase = 'https://ssd.jpl.nasa.gov/api/horizons.api' 

140 data = None 

141 

142 values = {'format': 'json', 

143 'EPHEM_TYPE': 'OBSERVER', 

144 'OBJ_DATA': 'YES', 

145 'COMMAND': f"'{target}'", 

146 'START_TIME': starttime, 

147 'STOP_TIME': stoptime, 

148 'STEP_SIZE': stepsize, 

149 'CENTER': '500@399', 

150 'QUANTITIES': f"'{quantities}'", 

151 'ANG_FORMAT': ang_format} 

152 

153 pardata = urlencode(values, doseq=True, encoding='utf-8') 

154 params = pardata.encode('ascii') 

155 # for debugging 

156 #print("params=", params) 

157 

158 req = Request(urlbase, params) 

159 datastr = None 

160 try: 

161 with urlopen(req, context=context, timeout=10) as response: 

162 datastr = response.read().decode() 

163 except URLError as e: 

164 if hasattr(e, 'reason'): 

165 msg = f'URLError: Failed to reach URL={urlbase} (reason: {e.reason})' 

166 casalog.post(msg, 'WARN') 

167 if hasattr(e, 'code'): 

168 msg = f'URLError: Couldn\'t fullfill the request to {urlbase} (code: {e.code})' 

169 casalog.post(msg, 'WARN') 

170 except socket.timeout as e: 

171 msg = f'Failed to reach URL={urlbase}. Socket timeout {e}' 

172 casalog.post(msg, 'WARN') 

173 

174 status = response.getcode() 

175 if datastr: 

176 try: 

177 data = ast.literal_eval(datastr) 

178 except RuntimeError as e: 

179 casalog.post(e, 'SEVERE') 

180 if status == 200: 

181 # 

182 # If the ephemeris data file was generated, write it to the output file: 

183 if rawdatafile != "": 

184 if os.path.exists(rawdatafile): 

185 os.remove(rawdatafile) 

186 if "result" in data: 

187 with open(rawdatafile, "w") as f: 

188 f.write(data["result"]) 

189 # Otherwise, the ephemeris file was not generated so output an error 

190 else: 

191 casalog.post("ERROR: No data found. Ephemeris data file not generated", 'WARN') 

192 else: 

193 raise RuntimeError(f'Could not retrieve the data. Error code:{status}:{response.msg}') 

194 else: 

195 data = None 

196 return data 

197 

198 

199def tocasatb(indata, outtable): 

200 """ 

201 convert a JPL-Horizons query results to a CASA table 

202 indata: either query data ('results') or file name that contains the result data 

203 outtable: output CASA table name 

204 """ 

205 # input data columns  

206 # Date__(UT)__HR:MN R.A.___(ICRF)___DEC ObsSub-LON ObsSub-LAT SunSub-LON SunSub-LAT ¥ 

207 # NP.ang NP.dist r rdot delta deldot S-T-O 

208 import re 

209 

210 # output CASA table columns  

211 # These are required columns for SetJy 

212 cols = { 

213 'MJD': {'header': 'Date__\(UT\)', 

214 'comment': 'date in MJD', 

215 'unit': 'd'}, 

216 'RA': {'header': 'R.A.', 

217 'comment': 'astrometric Right Ascension (ICRF)', 

218 'unit': 'deg'}, 

219 'DEC': {'header': 'DEC', 

220 'comment': 'astrometric Declination (ICRF)', 

221 'unit': 'deg'}, 

222 'Rho': {'header': 'delta', 

223 'comment': 'geocentric distance', 

224 'unit': 'AU'}, 

225 'RadVel': {'header': 'deldot', 

226 'comment': 'geocentric distance rate', 

227 'unit': 'AU/d'}, 

228 'NP_ang': {'header': 'NP.ang', 

229 'comment': 'North-Pole pos. angle', 

230 'unit': 'deg'}, 

231 'NP_dist': {'header': 'NP.dist', 

232 'comment': 'North-Pole ang. distance', 

233 'unit': 'arcsec'}, 

234 'DiskLong': {'header': 'ObsSub-LON', 

235 'comment': 'Sub-observer longitude', 

236 'unit': 'deg'}, 

237 'DiskLat': {'header': r'ObsSub-LAT', 

238 'comment': 'Sub-observer latitude', 

239 'unit': 'deg'}, 

240 'Sl_lon': {'header': 'SunSub-LON', 

241 'comment': 'Sub-Solar longitude', 

242 'unit': 'deg'}, 

243 'Sl_lat': {'header': r'SunSub-LAT', 

244 'comment': 'Sub-Solar longitude', 

245 'unit': 'deg'}, 

246 'r': {'header': 'r', 

247 'comment': 'heliocentric distance', 

248 'unit': 'AU'}, 

249 'rdot': {'header': 'rdot', 

250 'comment': 'heliocentric distance rate', 

251 'unit': 'km/s'}, 

252 'phang': {'header': 'S-T-O', 

253 'comment': 'phase angle', 

254 'unit': 'deg'} 

255 } 

256 

257 # do in a two-step 

258 # read the original query result dictionary containing ephemeris data and save the data part 

259 # to a temporary text file. Then re-read the temporary file to a casa table 

260 # with the approriate data format that setjy and measure expect. 

261 tempfname = 'temp_ephem_'+str(os.getpid())+'.dat' 

262 tempconvfname = 'temp_ephem_conv_'+str(os.getpid())+'.dat' 

263 try: 

264 exceedthelinelimit = None 

265 ambiguousname = None 

266 queryerrmsg = '' 

267 ephemdata = None 

268 # Scan the original data 

269 if isinstance(indata, dict) and 'result' in indata: 

270 #print("ephem data dict") 

271 ephemdata = indata['result'] 

272 elif isinstance(indata, str): 

273 if os.path.exists(indata): 

274 with open(indata, 'r') as infile: 

275 #ephemdata = infile.readlines() 

276 ephemdata = infile.read() 

277 # the relevant information in the main header section is extracted as 

278 else: 

279 raise IOError(f'{indata} does not exist') 

280 else: 

281 raise IOError(f'{indata} should be a dictionary contains ephemeris data under'+ 

282 ' the keyword, "results" or the text file name contains the raw JPL-Horizons'+ 

283 ' query results') 

284 # it is read line by line. The ephemeris data is stored in the separate text file 

285 # to be further processed in subsequent steps. 

286 with open(tempfname, 'w') as outfile: 

287 # Some initialization 

288 headerdict = {} 

289 # jplhorizonsdataIdFound = False 

290 datalines = 0 

291 readthedata = False 

292 startmjd = None 

293 endmjd = None 

294 incolnames = None 

295 # multiple entries (in different units) may exit for orb. per. 

296 foundorbper = False 

297 ### 

298 lcnt = 0 

299 # for lnum, line in enumerate(infile): 

300 for lnum, line in enumerate(ephemdata.split('\n')): 

301 # JPL-Horizons data should contain this line at the beginning 

302 if re.search(r'JPL/HORIZONS', line): 

303 # jplhorizondataIdFound = True 

304 casalog.post("Looks like JPL-Horizons data","INFO2") 

305 elif re.search(r'^\s*Ephemeris\s+', line): # date the data file was retrieved and created 

306 #m = re.search(r'(API_USER\S+\s+(\S+)\s+([0-9]+)\s+(\S+)\s+(\S+)') 

307 (_, _, _, wday, mon, day, tm, year, _) = re.split(' +', line, 8) 

308 # date info in 3-7th items 

309 

310 try: 

311 float(mon) 

312 nmon = mon 

313 except: 

314 tonummon=time.strptime(mon,"%b") 

315 nmon = f"{tonummon.tm_mon:02d}" 

316 day2 = f"{int(day):02d}" 

317 headerdict['VS_CREATE'] = year + '/' + nmon + '/' + day2 + '/' + tm[0:5] 

318 # VS_DATE - use the current time to indicate the time CASA table is created 

319 headerdict['VS_DATE'] = time.strftime('%Y/%m/%d/%H:%M',time.gmtime() ) 

320 headerdict['VS_TYPE'] = 'Table of comet/planetary positions' 

321 # VERSION stored in the output table may be incremented in the future. 

322 # For now, it is fixed, but it may be incremented from 0003 to 0004 to indiate 

323 # this new code is used to convert the jpl horizons data to a table. 

324 # ver 0004.0001 - new keyword: ephemeris_source, posrefsys 'ICRF/J2000' -> "ICRS'  

325 headerdict['VS_VERSION'] = '0004.0001' 

326 # target object name 

327 elif re.match(r'^[>\s]*Target body name', line): 

328 m = re.match(r'^[>\s]*Target body name:\s+(\S+)\s+(\S*)\s+\{(\w+)\:\s*(\S+)\}', line) 

329 if m: 

330 headerdict['NAME'] = m[1] 

331 if len(m.groups())==4 and m[3]=='source': 

332 headerdict['ephemeris_source'] = m[4] 

333 # start time (of the requested time range) 

334 elif re.search(r'Start time', line): 

335 m = re.match(r'^[>\s]*Start time\s+\S+\s+\S+\s+(\S+)\s+(\S+)\s+(\w+)', line) 

336 if m: 

337 startmjd = _qa.totime(m[1] + '/' + m[2]) 

338 #--This info will not be used but left here since it might be useful fo debugging. 

339 # # end time (of the requested time range) 

340 elif re.search(r'Stop time', line): 

341 m = re.match(r'^[>\s]*Stop time\s+\S+\s+\S+\s+(\S+)\s+(\S+)\s+(\w+)', line) 

342 if m: 

343 endmjd = _qa.totime(m[1] + '/' + m[2]) 

344 # date increment 

345 elif re.search(r'Step-size', line): 

346 m = re.match(r'^[>\s]*Step-size\s+\S+\s+(\S+)\s+(\w+)', line) 

347 if m: 

348 unit = m[2] 

349 if unit == 'minutes': 

350 # this is the default unit returned by the JPL-Horizons 

351 theunit = 'min' 

352 elif unit == 'hours': 

353 theunit = 'h' 

354 elif unit == 'days': 

355 theunit = 'd' 

356 elif unit == 'steps': 

357 theunit = 'steps' 

358 elif unit == 'calendar': 

359 raise RuntimeError('Unit of Step-size in calendar month or year is not supported.') 

360 else: 

361 raise RuntimeError(f'Unit of Step-size, {unit} is not recognized.') 

362 if theunit == 'd': 

363 dmjd = m[1] 

364 elif theunit == 'steps': 

365 # print('endtime=',endmjd['value']) 

366 dmjd = (endmjd['value'] - startmjd['value'])/int(m[1]) 

367 else: 

368 dmjd = _qa.convert(_qa.totime(m[1] + theunit), 'd')['value'] 

369 headerdict['dMJD'] = dmjd 

370 if startmjd is not None: # start mjd should be available before step-size line 

371 # MJD0 = firstMJD - dMJD (as defined casacore MeasComet documentation) 

372 headerdict['MJD0'] = startmjd['value'] - dmjd 

373 elif re.search(r'Center geodetic', line): 

374 m = re.match(r'^[>\s]*Center geodetic\s*: ([-+0-9.]+,\s*[-+0-9.]+,\s*[-+0-9.]+)', line) 

375 if m: 

376 long_lat_alt = m[1].split(',') 

377 headerdict['GeoLong'] = float(long_lat_alt[0]) 

378 headerdict['GeoLat'] = float(long_lat_alt[1]) 

379 headerdict['GeoDist'] = float(long_lat_alt[2]) 

380 # obs location 

381 elif re.search(r'Center-site name', line): 

382 m = re.match(r'^[>\s]*Center-site name:(\s*)([a-zA-Z\s*\/\-\(\)\,]+)', line) 

383 if m: 

384 # For the topocentric location, currently only ALMA, VLA and GBT are translated  

385 # to the proper observatory name which recongnized by Measures. 

386 # The key part is the name used by JPL-Horizons. 

387 observatories = {'ALMA':'ALMA', 'VLA':'VLA', 'Green Bank':'GBT'} 

388 if m[2]: 

389 headerdict['obsloc'] = m[2] 

390 for obs in observatories: 

391 if obs in m[2]: 

392 headerdict['obsloc'] = observatories[obs] 

393 break 

394 # Assume the query is made in ICRF reference frame 

395 # and for casacore measures, this will be 'ICRS' 

396 headerdict['posrefsys'] = 'ICRS' 

397 elif re.search(r'Target radii', line): 

398 m = re.match(r'^[>\s]*Target radii\s*:\s*([0-9.]+\s*x\s*[0-9.]+\s*x\s*[0-9.]+)\s*km.*|' 

399 '^[>/s]*Target radii\s*:\s*([0-9.]+)\s*km', line) 

400 if m: 

401 matchlist = m.groups() 

402 radiiarr = np.zeros(3) 

403 if len(matchlist)==2: 

404 if m[2] is None: 

405 radiiarr = np.asarray(np.array(m[1].split(' x ')), dtype=np.float64) 

406 headerdict['radii'] = {'unit': 'km', 'value': radiiarr} 

407 elif m[1] is None: 

408 radiiarr = np.array([m[2],m[2],m[2]], dtype=np.float64) 

409 headerdict['radii'] = {'unit': 'km', 'value': radiiarr} 

410 else: 

411 casalog.post(f"Unexpected number or matches for Target radii:{m.groups} (expected 2)", 'WARN') 

412 #rotational period (few pattens seem to exist) 

413 elif re.search(r'rot. period|Rotational period', line): 

414 m = re.search(r'rot. period\s+\S*=\s*([0-9.]+h\s*[0-9.]+m\s*[0-9.]+\s*s)|' 

415 'rot. period\s+\S*=\s*([0-9.]+)(?:\s*\+-[0-9.]+)?\s*([dh])|' 

416 'Rotational period\s*=\s+Synchronous', line) 

417 if m: 

418 if m[0].find('Synchronous') != -1: 

419 headerdict['rot_per'] = 'Synchronous' 

420 else: 

421 if len(m.groups()) == 3: 

422 if m[1] is None: 

423 headerdict['rot_per'] = _qa.quantity(m[2] + m[3]) 

424 elif m[2] is None and m[3] is None: 

425 #subm = re.search(r'([0-9]+)h\s*([0-9]+)m\s*([0-9.]+)\s*s',m[1]) 

426 headerdict['rot_per'] = _qa.convert(re.sub(r'\s+','',m[1]),'h') 

427 # another variation of rotational period entry 

428 elif re.search(r'rot.\s+per.', line): 

429 m = re.search(r'rot.\s+per.\s*=\s*([0-9.]+)\s*([dh])',line) 

430 if m: 

431 headerdict['rot_per'] = _qa.quantity(m[1]+m[2]) 

432 # rotational period for asteroids 

433 elif re.search(r'ROTPER', line): 

434 m = re.search(r'ROTPER=\s*([0-9.]+)\s*', line) 

435 if m: 

436 headerdict['rot_per'] = {'unit': 'h', 'value': float(m[1])} 

437 elif re.search(r'orbit period|orb per|orb. per.|orbital period', line.lower()) \ 

438 and not foundorbper: 

439 m = re.search(r'Orbital period\s*[=~]\s*([-0-9.]+)\s*(\w+)\s*|' 

440 'orb. per., (\w)\s+=\s+([0-9.]+)\s+|' 

441 'orb per\s+=\s+([0-9.]+)\s+(\w+)\s+|' 

442 'orbit period\s+=\s+([0-9.]+)\s+(\w+)\s+', line) 

443 if m: 

444 if m[0].find('Orbital period') != -1: 

445 headerdict['orb_per'] = {'unit': m[2], 'value': float(m[1])} 

446 elif m[0].find('orb. per.') != -1: 

447 headerdict['orb_per'] = {'unit': m[3], 'value': float(m[4])} 

448 elif m[0].find('orb per') != -1: 

449 headerdict['orb_per'] = {'unit': m[6], 'value': float(m[5])} 

450 elif m[0].find('orbit period') != -1: 

451 headerdict['orb_per'] = {'unit': m[8], 'value': float(m[7])} 

452 if 'orb_per' in headerdict: 

453 foundorbper = True 

454 elif re.search(r'Mean [Tt]emperature', line): 

455 m = re.search(r'Mean [Tt]emperature\s+\(K\)\s+=\s+([0-9.]+)\s+',line) 

456 if m: 

457 headerdict['T_mean'] = {'unit':'K', 'value':float(m[1])} 

458 # start reading data 

459 elif re.search(r'\s*Date__\(UT\)', line): 

460 incolnames = line.split() 

461 elif re.search(r'\$\$SOE', line): 

462 readthedata = True 

463 elif re.search(r'\$\$EOE', line): 

464 readthedata = False 

465 elif readthedata: 

466 datalines += 1 

467 outfile.write(line + '\n') 

468 elif re.search(r'Projected output length', line): 

469 # this message occurs when requested data is too large 

470 exceedthelinelimit = line 

471 elif re.search(r'Multiple major-bodies match', line): 

472 ambiguousname = line 

473 else: 

474 pass 

475 lcnt += 1 

476 if 'radii' in headerdict: 

477 radiival = headerdict['radii']['value'] 

478 meanrad = _mean_radius(radiival[0], radiival[1], radiival[2]) 

479 headerdict['meanrad'] = {'unit': 'km', 'value': meanrad} 

480 if datalines == 0: 

481 casalog.post("No ephemeris data was found", "WARN") 

482 queryerrmsg = exceedthelinelimit if exceedthelinelimit != None else ambiguousname 

483 

484 if queryerrmsg: 

485 raise RuntimeError(f'Error occurred at the query:{queryerrmsg}') 

486 else: 

487 raise RuntimeError(f'Error occurred at the query') 

488 else: 

489 casalog.post(f'Number of data lines={datalines}') 

490 casalog.post(f'Number of all lines in the file={lcnt}') 

491 #print("headerdict=", headerdict) 

492 # output to a casa table 

493 

494 # check the input data columns and stored the order as indices 

495 foundncols = 0 

496 indexoffset = 0 

497 colkeys = {} 

498 #print('incolnames=',incolnames) 

499 if incolnames is not None: 

500 for outcolname in cols: 

501 # all colnames in cols should have unit defined. 

502 if 'unit' in cols[outcolname]: 

503 colkeys[outcolname] = np.array([cols[outcolname]['unit']]) 

504 else: 

505 raise KeyError(f'Missing unit for {outcolname}.') 

506 inheadername = cols[outcolname]['header'] 

507 # expect date is in the first column (date and mm:hh seperated by spaces) 

508 if outcolname == 'MJD': 

509 for incol in incolnames: 

510 if re.search(inheadername + '.+', incol): 

511 cols[outcolname]['index'] = incolnames.index(incol) 

512 foundncols += 1 

513 indexoffset = 1 

514 if 'index' not in cols[outcolname]: 

515 casalog.post("Cannot find the Date column", "WARN") 

516 elif outcolname == 'RA': 

517 for incol in incolnames: 

518 if re.search(inheadername + '.+(ICRF).+', incol): 

519 cols[outcolname]['index'] = incolnames.index(incol) + indexoffset 

520 foundncols += 1 

521 

522 if 'index' not in cols[outcolname]: 

523 casalog.post("Cannot find the astrometric RA and Dec column", "WARN") 

524 elif outcolname == 'DEC': 

525 if 'index' in cols['RA']: 

526 # Dec data col is next to RA data col 

527 cols[outcolname]['index'] = cols['RA']['index'] + 1 

528 # add additional offset for index (4 data columns at this point) 

529 indexoffset +=1 

530 foundncols += 1 

531 else: 

532 if inheadername in incolnames: 

533 cols[outcolname]['index'] = incolnames.index(inheadername) + indexoffset 

534 foundncols += 1 

535 else: 

536 casalog.post(f"Cannot find {inheadername}", "WARN") 

537 

538 #print(cols) 

539 casalog.post(f"expected n cols = {len(cols)} ") 

540 casalog.post(f"foundncols = {foundncols}") 

541 if foundncols == len(cols): 

542 # Format the data to comply with measure/setjy 

543 casalog.post("Found all the required columns") 

544 with open(tempconvfname, 'w') as outf, open(tempfname, 'r') as inf: 

545 ndata = 0 

546 earliestmjd = None 

547 mjd = None 

548 prevmjd = None 

549 for line in inf: 

550 outline = '' 

551 sep = ' ' 

552 extraindexoffset = 0 

553 rawtempdata = line.split() 

554 tempdata = ['-999.0' if x == 'n.a.' else x for x in rawtempdata] 

555 # Check if there are extra codes (solar/lunar presence codes) after the date. 

556 # The extra codes may be present for topocentric observer location. Even in that case 

557 # not every data line has the code so it needs to be checked line by line. 

558 ntempdata = len(tempdata) 

559 nincolnames = len(incolnames) + 2 # 2 for date_time and ra_dec splits  

560 if ntempdata > nincolnames : 

561 # assume there is the extra codes (no space between solar and lunar presence code) 

562 if ntempdata - nincolnames == 1: 

563 extraindexoffset = 1 

564 else: 

565 raise RuntimeError(f'Unexpected number of data items was detected:'+ 

566 f' {ntempdata} while expected {nincolnames}') 

567 # construct mjd from col 1 (calendar date) + col 2 (time) 

568 caldatestr = tempdata[0] + ' ' + tempdata[1] 

569 mjd = _qa.totime(caldatestr) 

570 if mjd == prevmjd: 

571 raise RuntimeError(f'Duplicated timestamp, {mjd}, is detected. This may occur when '+ 

572 'a time range is specified in MJD with a short time interval. If that is the case, '+ 

573 'try calendar date+time string for the time range.') 

574 outline += str(mjd['value']) + sep 

575 # position 

576 rad = tempdata[cols['RA']['index'] + extraindexoffset] 

577 decd = tempdata[cols['DEC']['index'] + extraindexoffset] 

578 outline += rad + sep + decd + sep 

579 # geocentric dist. (Rho) 

580 delta = tempdata[cols['Rho']['index'] + extraindexoffset] 

581 outline += delta + sep 

582 # geocentric range rate (RadVel) 

583 valinkmps = tempdata[cols['RadVel']['index'] + extraindexoffset] 

584 deldot = _qa.convert(_qa.quantity(valinkmps+'km/s'), 'AU/d' )['value'] 

585 outline += str(deldot) + sep 

586 # NP_ang & NP_dist 

587 npang = tempdata[cols['NP_ang']['index'] + extraindexoffset] 

588 npdist = tempdata[cols['NP_dist']['index'] + extraindexoffset] 

589 outline += npang + sep + npdist + sep 

590 # DiskLong & DiskLat 

591 disklong = tempdata[cols['DiskLong']['index'] + extraindexoffset] 

592 disklat = tempdata[cols['DiskLat']['index'] + extraindexoffset] 

593 outline += disklong + sep + disklat + sep 

594 # sub-long & sub-lat 

595 sllon = tempdata[cols['Sl_lon']['index'] + extraindexoffset] 

596 sllat = tempdata[cols['Sl_lat']['index'] + extraindexoffset] 

597 outline += sllon + sep + sllat + sep 

598 # r, rot 

599 r = tempdata[cols['r']['index'] + extraindexoffset] 

600 rdot = tempdata[cols['rdot']['index'] + extraindexoffset] 

601 outline += r + sep + rdot + sep 

602 # S-T-O 

603 phang = tempdata[cols['phang']['index'] + extraindexoffset] 

604 outline += phang 

605 outf.write(outline + '\n') 

606 

607 ndata += 1 

608 if ndata == 1: 

609 earliestmjd = mjd 

610 prevmjd = mjd 

611 # record first and last mjd in the data 

612 headerdict['earliest'] = _me.epoch('UTC',earliestmjd) 

613 headerdict['latest'] = _me.epoch('UTC', mjd) 

614 

615 else: 

616 missingcols = len(cols) - foundncols 

617 casalog.post(r"Missing {missingcols}") 

618 

619 # final step: convert to a CASA table 

620 dtypes = np.array(['D' for _ in range(len(cols))]) 

621 _tb.fromascii(outtable, tempconvfname, sep=' ', columnnames=list(cols.keys()), 

622 datatypes=dtypes.tolist()) 

623 _tb.done() 

624 

625 # fill keyword values in the ephem table 

626 if os.path.exists(outtable): 

627 _fill_keywords_from_dict(headerdict, colkeys, outtable) 

628 casalog.post(f"Output is written to a CASA table, {outtable}") 

629 else: 

630 raise RuntimeError("Error occured. The output table, " + outtable + "is not generated") 

631 else: #incolnames is None 

632 raise RuntimeError("No data header was found.") 

633 

634 except RuntimeError as e: 

635 casalog.post(str(e),"SEVERE") 

636 

637 finally: 

638 tempfiles = [tempfname, tempconvfname] 

639 _clean_up(tempfiles) 

640 

641# This needs to change to read the ephem data from the generated casa table. 

642# Will be called in fill_keywords_from_dict() 

643 

644# Copied from JPLephem_reader2.py for mean_radius calc (when no NP-ang, NP_dist) 

645# Note: this is not fully verified for the correctness but has been used  

646# through JPLephem_reader2 for processing the ephemeris data to a CASA table. 

647# For setjy, solar_system_setjy does own mean radius calculation and  

648# meanrad is the table is not used for that. 2022.01.12 TT 

649def _mean_radius(a, b, c): 

650 """ 

651 Return the average apparent mean radius of an ellipsoid with semiaxes 

652 a >= b >= c. 

653 "average" means average over time naively assuming the pole orientation 

654 is uniformly distributed over the whole sphere, and "apparent mean radius" 

655 means a radius that would give the same area as the apparent disk. 

656 """ 

657 # This is an approximation, but it's not bad. 

658 # The exact equations for going from a, b, c, and the Euler angles to the 

659 # apparent ellipse are given in Drummond et al., Icarus, 1985a. 

660 # It's the integral over the spin phase that I have approximated, so the 

661 # approximation is exact for b == a and appears to hold well for b << a. 

662 R = 0.5 * c ** 2 * (1.0 / b ** 2 + 1.0 / a ** 2) # The magic ratio. 

663 if R < 0.95: 

664 sqrt1mR = sqrt(1.0 - R) 

665 # There is fake singularity (RlnR) at R = 0, but it is unlikely to 

666 # be a problem. 

667 try: 

668 Rterm = 0.5 * R * log((1.0 + sqrt1mR) / (1.0 - sqrt1mR)) / sqrt1mR 

669 except RuntimeError: 

670 Rterm = 0.0 

671 else: 

672 # Use a (rapidly converging) series expansion to avoid a fake 

673 # singularity at R = 1. 

674 Rterm = 1.0 # 0th order 

675 onemR = 1.0 - R 

676 onemRtothei = 1.0 

677 for i in range(1, 5): # Start series at 1st order. 

678 onemRtothei *= onemR 

679 Rterm -= onemRtothei / (0.5 + 2.0 * i ** 2) 

680 avalfabeta = 0.5 * a * b * (1.0 + Rterm) 

681 return sqrt(avalfabeta) 

682 

683 

684def _mean_radius_with_known_theta(disklat, radii): 

685 """ 

686 Mean radius calculation extracted from solar_system_setjy 

687 """ 

688 Req = 1000.0 * (radii[0] + radii[1]) / 2 

689 Rp = 1000.0 * radii[2] 

690 

691 Rmean = 0.0 

692 index = 0 

693 for index, lat in enumerate(disklat): 

694 if lat == -999.0: 

695 lat = 0.0 

696 rlat = lat * (180. / np.pi) 

697 Rpap = sqrt(Req * Req * sin(rlat) ** 2.0 + 

698 Rp * Rp * cos(rlat) ** 2.0) 

699 # Rmean += ( sqrt (Rpap * Req) - Rmean )/ (index + 1.0) 

700 Rmean += sqrt(Rpap * Req) 

701 

702 return (Rmean / (index + 1.0)) / 1000.0 

703 

704 

705def _fill_keywords_from_dict(keydict, colkeys, tablename): 

706 # open tb 

707 # get ra,dec,np_ra, np_dec, and radii in the keyword 

708 # call mod version of mean_radius_with_known_theta 

709 orderedmainkeys = ['VS_CREATE','VS_DATE','VS_TYPE','VS_VERSION','NAME', 

710 'MJD0','dMJD','GeoDist','GeoLat','GeoLong','obsloc', 

711 'posrefsys','ephemeris_source','earliest','latest','radii', 

712 'meanrad','orb_per','rot_per','T_mean'] 

713 try: 

714 _tb.open(tablename, nomodify=False) 

715 disklat = _tb.getcol('DiskLat') 

716 if 'radii' in keydict: 

717 radiival = keydict['radii']['value'] 

718 calc_meanrad = _mean_radius_with_known_theta(disklat, radiival) 

719 if calc_meanrad and ('meanrad' in keydict): 

720 keydict['meanrad']['value'] = calc_meanrad 

721 if 'rot_per' in keydict and 'orb_per' in keydict and keydict['rot_per'] == 'Synchronous': 

722 keydict['rot_per'] = keydict['orb_per'] 

723 sortedkeydict = {k: keydict[k] for k in orderedmainkeys if k in keydict} 

724 #print('sorteddict=',sortedkeydict) 

725 for k in sortedkeydict: 

726 _tb.putkeyword(k, sortedkeydict[k]) 

727 datacolnames = _tb.colnames() 

728 for col in datacolnames: 

729 if col in colkeys: 

730 _tb.putcolkeyword(col, 'QuantumUnits', colkeys[col]) 

731 # add table info required by measComet 

732 maintbinfo = {'readme': 'Derivied by jplhorizons-query from JPL-Horizons API ' 

733 '(https://ssd.jpl.nasa.gov/api/horizons.api)', 

734 'subType':'Comet', 

735 'type': 'IERS'} 

736 _tb.putinfo(maintbinfo) 

737 _tb.flush() 

738 _tb.done() 

739 except RuntimeError: 

740 casalog.post('Internal error: Cannot add the data in keywords', "WARN") 

741 

742def _clean_up(filelist): 

743 """ 

744 Clean up the temporary files  

745 """ 

746 for f in filelist: 

747 if os.path.exists(f): 

748 os.remove(f)