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-10-31 19:10 +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
10from casatools import table, quanta, measures
11from casatasks import casalog
13_tb = table()
14_qa = quanta()
15_me = measures()
17# debug
18_debug = False
20def gethorizonsephem(objectname, starttime, stoptime, incr, outtable, asis=False, rawdatafile=''):
21 """
22 Main driver function for ephemeris data query from JPL-Horizons
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 """
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'
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)
92 try:
93 step_size = incr.replace(' ', '')
94 except ValueError as e:
95 casalog.post(e)
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')
115def queryhorizons(target, starttime, stoptime, stepsize, quantities, ang_format, rawdatafile=''):
116 """
117 Submit a query to the JPL-Horizons API
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)
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
138 context = ssl.create_default_context(cafile=certifi.where())
139 urlbase = 'https://ssd.jpl.nasa.gov/api/horizons.api'
140 data = None
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}
153 pardata = urlencode(values, doseq=True, encoding='utf-8')
154 params = pardata.encode('ascii')
155 # for debugging
156 #print("params=", params)
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')
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
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
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 }
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
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
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
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
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")
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')
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)
615 else:
616 missingcols = len(cols) - foundncols
617 casalog.post(r"Missing {missingcols}")
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()
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.")
634 except RuntimeError as e:
635 casalog.post(str(e),"SEVERE")
637 finally:
638 tempfiles = [tempfname, tempconvfname]
639 _clean_up(tempfiles)
641# This needs to change to read the ephem data from the generated casa table.
642# Will be called in fill_keywords_from_dict()
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)
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]
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)
702 return (Rmean / (index + 1.0)) / 1000.0
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")
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)