|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 """ 00002 casapy functions for converting ASCII ephemerides from JPL-Horizons into 00003 CASA tables and installing them where casapy can find them. 00004 00005 jplfiles_to_repository() puts it all together, so it is most likely the 00006 function you want. 00007 00008 There are various utilities like convert_radec, datestr*, get_num_from_str, 00009 mean_radius*, and construct_tablepath defined in here as well. 00010 """ 00011 00012 from glob import glob 00013 import os 00014 import re 00015 import scipy.special 00016 import time # We can always use more time. 00017 00018 from taskinit import gentools, qa, casalog 00019 me = gentools(['me'])[0] 00020 00021 from dict_to_table import dict_to_table 00022 00023 # Possible columns, as announced by their column titles. 00024 # The data is scooped up by 'pat'. Either use ONE group named by the column 00025 # key, or mark it as unwanted. '-' is not valid in group names. 00026 # Leading and trailing whitespace will be taken care of later. 00027 # Sample lines: 00028 # Date__(UT)__HR:MN R.A.___(ICRF/J2000.0)___DEC Ob-lon Ob-lat Sl-lon Sl-lat NP.ang NP.dist r rdot delta deldot S-T-O L_s 00029 # 2010-May-01 00:00 09 01 43.1966 +19 04 28.673 286.52 18.22 246.99 25.34 358.6230 3.44 1.661167637023 -0.5303431 1.28664311447968 15.7195833 37.3033 84.50 00030 # 00031 # some mod to label names and comments so that they corresponds to 00032 # JPL-Horizons naming comvension 00033 cols = { 00034 'MJD': {'header': r'Date__\(UT\)__HR:MN', 00035 'comment': 'date', 00036 'pat': r'(?P<MJD>\d+-\w+-\d+ \d+:\d+)'}, 00037 'RA': {'header': r'R.A._+\([^)]+', 00038 'comment': 'Right Ascension', 00039 'pat': r'(?P<RA>(\d+ \d+ )?\d+\.\d+)'}, # require a . for safety 00040 'DEC': {'header': r'\)_+DEC.', 00041 'comment': 'Declination', 00042 'pat': r'(?P<DEC>([-+]?\d+ \d+ )?[-+]?\d+\.\d+)'}, 00043 'NP_ang': {'header': r'NP\.ang', 00044 'comment': 'North-Pole pos. angle', 00045 'pat': r'(?P<NP_ang>[0-9.]+|n\.a\.)', 00046 'unit': 'deg'}, 00047 'NP_dist': {'header': r'NP\.dist', 00048 'comment': 'North-Pole ang. distance', 00049 'pat': r'(?P<NP_dist>[-+0-9.]+|n\.a\.)', 00050 'unit': 'arcsec'}, 00051 'illu': {'header': r'Illu%', 00052 #'comment': 'Illumination', 00053 'comment': 'Illuminated fraction', 00054 'pat': r'(?P<illu>[0-9.]+)', 00055 'unit': r'%'}, 00056 # put back to original heading... 00057 'DiskLong': {'header': r'Ob-lon', 00058 #'Obs_Long': {'header': r'Ob-lon', 00059 'comment': 'Sub-observer longitude', 00060 # 'pat': r'(?P<Obs_Long>[0-9.]+|n\.a\.)', 00061 'pat': r'(?P<DiskLong>[0-9.]+|n\.a\.)', 00062 'unit': 'deg'}, 00063 'DiskLat': {'header': r'Ob-lat', 00064 #'Obs_Lat': {'header': r'Ob-lat', 00065 'comment': 'Sub-observer latitude', 00066 # 'pat': r'(?P<Obs_Lat>[-+0-9.]+|n\.a\.)', 00067 'pat': r'(?P<DiskLat>[-+0-9.]+|n\.a\.)', 00068 'unit': 'deg'}, 00069 'Sl_lon': {'header': r'Sl-lon', 00070 'comment': 'Sub-Solar longitude', 00071 'pat': r'(?P<Sl_lon>[0-9.]+|n\.a\.)', 00072 'unit': 'deg'}, 00073 'Sl_lat': {'header': r'Sl-lat', 00074 'comment': 'Sub-Solar longitude', 00075 'pat': r'(?P<Sl_lat>[-+0-9.]+|n\.a\.)', 00076 'unit': 'deg'}, 00077 00078 # These are J2000 whether or not ra and dec are apparent directions. 00079 'NP_RA': {'header': r'N\.Pole-RA', 00080 'comment': 'North Pole right ascension', 00081 'pat': r'(?P<NP_RA>(\d+ \d+ )?\d+\.\d+)'}, # require a . for safety 00082 'NP_DEC': {'header': r'N\.Pole-DC', 00083 'comment': 'North Pole declination', 00084 'pat': r'(?P<NP_DEC>([-+]?\d+ \d+ )?[-+]?\d+\.\d+)'}, 00085 00086 'r': {'header': 'r', 00087 'comment': 'heliocentric distance', 00088 'unit': 'AU', 00089 'pat': r'(?P<r>[0-9.]+)'}, 00090 'rdot': {'header': 'rdot', 00091 #'comment': 'heliocentric velocity', 00092 'comment': 'heliocentric distance rate', 00093 'unit': 'km/s', 00094 'pat': r'(?P<rdot>[-+0-9.]+)'}, 00095 # 'unwanted': True}, 00096 'Rho': {'header': 'delta', 00097 'comment': 'geocentric distance', 00098 'unit': 'AU', 00099 'pat': r'(?P<Rho>[0-9.]+)'}, 00100 'RadVel': {'header': 'deldot', 00101 #'comment': 'Radial velocity relative to the observer', 00102 'comment': 'geocentric distance rate', 00103 'pat': r'(?P<RadVel>[-+0-9.]+)', 00104 'unit': 'km/s'}, 00105 'phang': {'header': 'S-T-O', 00106 'comment': 'phase angle', 00107 'unit': 'deg', 00108 'pat': r'(?P<phang>[0-9.]+)'}, 00109 'ang_sep': {'header': 'ang-sep/v', 00110 'comment': 'Angular separation from primary', 00111 'pat': r'(?P<ang_sep>[0-9.]+/.)'}, # arcsec, "visibility code". 00112 # t: transiting primary 00113 # O: occulted by primary 00114 # p: partial umbral eclipse 00115 # P: occulted partial umbral eclipse 00116 # u: total umbral eclipse 00117 # U: occulted total umbral eclipse 00118 # *: none of the above 00119 # -: target is primary 00120 'L_s': {'header': 'L_s', # 08/2010: JPL does not supply this and 00121 'unit': 'deg', # says they cannot. Ask Bryan Butler. 00122 'comment': 'Season angle', 00123 'pat': r'(?P<L_s>[-+0-9.]+)'} 00124 } 00125 00126 def readJPLephem(fmfile): 00127 """ 00128 Reads a JPL Horizons text file (see 00129 http://ssd.jpl.nasa.gov/horizons.cgi#top ) for a solar system object and 00130 returns various quantities in a dictionary. The dict will be blank ({}) if 00131 there is a failure. 00132 """ 00133 retdict = {} 00134 casalog.origin('readJPLephem') 00135 00136 # Try opening fmfile now, because otherwise there's no point continuing. 00137 try: 00138 ephem = open(fmfile) 00139 except IOError: 00140 casalog.post("Could not open ephemeris file " + fmfile, 00141 priority="SEVERE") 00142 return {} 00143 00144 # Setup the regexps. 00145 00146 # Headers (one time only things) 00147 00148 # Dictionary of quantity label: regexp pattern pairs that will be searched 00149 # for once. The matching quantity will go in retdict[label]. Only a 00150 # single quantity (group) will be retrieved per line. 00151 headers = { 00152 'NAME': {'pat': r'^Target body name:\s+\d*\s*(\w+)'}, # object name, w.o. number 00153 'ephtype': {'pat': r'\?s_type=1#top>\]\s*:\s+\*(\w+)'}, # e.g. OBSERVER 00154 'obsloc': {'pat': r'^Center-site name:\s+(\w+)'}, # e.g. GEOCENTRIC 00155 # Catch either an explicit mean radius or a solitary target radius. 00156 'meanrad': {'pat': r'(?:Mean radius \(km\)\s*=|^Target radii\s*:)\s*([0-9.]+)(?:\s*km)?\s*$', 00157 'unit': 'km'}, 00158 # Triaxial target radii 00159 #'radii': {'pat': r'Target radii\s*:\s*([0-9.]+\s*x\s*[0-9.]+\s*x\s*[0-9.]+)\s*km.*Equator, meridian, pole', 00160 'radii': {'pat': r'Target radii\s*:\s*([0-9.]+\s*x\s*[0-9.]+\s*x\s*[0-9.]+)\s*km.*Equator, meridian, pole|Target radii\s*:\s*([0-9.]+)\s*km\s*', 00161 'unit': 'km'}, 00162 'T_mean': {'pat': r'Mean Temperature \(K\)\s*=\s*([0-9.]+)', 00163 'unit': 'K'}, 00164 00165 # Figure out the units later. 00166 'rot_per': {'pat': r'(?i)(?<!Inferred )\b(rot(ation(al)?|\.)?\s*per.*=\s*([-0-9.]+\s*[dhr]*|Synchronous))'}, 00167 'orb_per': {'pat': r'Orbital period((, days)?\s*=\s*[-0-9.]+\s*[dhr](\s*\(?R\)?)?)'}, 00168 00169 # MeasComet does not read units for these! E-lon(deg), Lat(deg), Alt(km) 00170 'GeoLong': {'pat': r'^Center geodetic\s*: ([-+0-9.]+,\s*[-+0-9.]+,\s*[-+0-9.]+)'}, 00171 'dMJD': {'pat': r'^Step-size\s*:\s*(.+)'}, 00172 00173 # request method v wday mth mday hh mm ss yyyy 00174 'VS_CREATE': {'pat': r'^Ephemeris / \w+ \w+ (\w+\s+\d+\s+\d+:\d+:\d+\s+\d+)'} 00175 } 00176 for hk in headers: 00177 headers[hk]['pat'] = re.compile(headers[hk]['pat']) 00178 00179 # Data ("the rows of the table") 00180 00181 # need date, r (heliocentric distance), delta (geocentric distance), and phang (phase angle). 00182 # (Could use the "dot" time derivatives for Doppler shifting, but it's 00183 # likely unnecessary.) 00184 datapat = r'^\s*' 00185 00186 stoppat = r'\$\$EOE$' # Signifies the end of data. 00187 00188 # Read fmfile into retdict. 00189 num_cols = 0 00190 in_data = False 00191 comp_mismatches = [] 00192 print_datapat = False 00193 # define interpretation of invalid values ('n.a.') 00194 invalid=-999. 00195 for line in ephem: 00196 if in_data: 00197 if re.match(stoppat, line): 00198 break 00199 matchobj = re.search(datapat, line) 00200 if matchobj: 00201 gdict = matchobj.groupdict() 00202 for col in gdict: 00203 if gdict[col]=='n.a.': 00204 gdict[col]=invalid 00205 retdict['data'][col]['data'].append(gdict[col]) 00206 if len(gdict) < num_cols: 00207 print "Partially mismatching line:" 00208 print line 00209 print "Found:" 00210 print gdict 00211 print_datapat = True 00212 else: 00213 print_datapat = True 00214 # Chomp trailing whitespace. 00215 comp_mismatches.append(re.sub(r'\s*$', '', line)) 00216 elif re.match(r'^\s*' + cols['MJD']['header'] + r'\s+' 00217 + cols['RA']['header'], line): 00218 # See what columns are present, and finish setting up datapat and 00219 # retdict. 00220 havecols = [] 00221 # extract coordinate ref info 00222 00223 m=re.match(r'(^\s*)(\S+)(\s+)('+cols['RA']['header']+')', line) 00224 coordref=m.group(4).split('(')[-1] 00225 cols['RA']['comment']+='('+coordref+')' 00226 cols['DEC']['comment']+='('+coordref+')' 00227 # Chomp trailing whitespace. 00228 myline = re.sub(r'\s*$', '', line) 00229 titleline = myline 00230 remaining_cols = cols.keys() 00231 found_col = True 00232 # This loop will terminate one way or another. 00233 while myline and remaining_cols and found_col: 00234 found_col = False 00235 #print "myline = '%s'" % myline 00236 #print "remaining_cols =", ', '.join(remaining_cols) 00237 for col in remaining_cols: 00238 if re.match(r'^\s*' + cols[col]['header'], myline): 00239 #print "Found", col 00240 havecols.append(col) 00241 remaining_cols.remove(col) 00242 myline = re.sub(r'^\s*' + cols[col]['header'], 00243 '', myline) 00244 found_col = True 00245 break 00246 datapat += r'\s+'.join([cols[col]['pat'] for col in havecols]) 00247 sdatapat = datapat 00248 casalog.post("Found columns: " + ', '.join(havecols)) 00249 datapat = re.compile(datapat) 00250 retdict['data'] = {} 00251 for col in havecols: 00252 if not cols[col].get('unwanted'): 00253 retdict['data'][col] = {'comment': cols[col]['comment'], 00254 'data': []} 00255 num_cols = len(retdict['data']) 00256 elif re.match(r'^\$\$SOE\s*$', line): # Start of ephemeris 00257 casalog.post("Starting to read data.", priority='INFO2') 00258 in_data = True 00259 else: 00260 #print "line =", line 00261 #print "looking for", 00262 for hk in headers: 00263 if not retdict.has_key(hk): 00264 matchobj = re.search(headers[hk]['pat'], line) 00265 if matchobj: 00266 if hk=='radii': 00267 mobjs=matchobj.groups() 00268 for gp in mobjs: 00269 if gp!=None: 00270 retdict[hk] = gp 00271 break 00272 break 00273 else: 00274 retdict[hk] = matchobj.group(1) # 0 is the whole line 00275 break 00276 ephem.close() 00277 00278 # If there were errors, provide debugging info. 00279 if comp_mismatches: 00280 print "Completely mismatching lines:" 00281 print "\n".join(comp_mismatches) 00282 if print_datapat: 00283 print "The apparent title line is:" 00284 print titleline 00285 print "datapat = r'%s'" % sdatapat 00286 00287 # Convert numerical strings into actual numbers. 00288 try: 00289 retdict['earliest'] = datestr_to_epoch(retdict['data']['MJD']['data'][0]) 00290 retdict['latest'] = datestr_to_epoch(retdict['data']['MJD']['data'][-1]) 00291 except Exception, e: 00292 print "Error!" 00293 if retdict.has_key('data'): 00294 if retdict['data'].has_key('MJD'): 00295 if retdict['data']['MJD'].has_key('data'): 00296 #print "retdict['data']['MJD']['data'] =", retdict['data']['MJD']['data'] 00297 print "retdict['data'] =", retdict['data'] 00298 else: 00299 print "retdict['data']['MJD'] has no 'data' key." 00300 print "retdict['data']['MJD'].keys() =", retdict['data']['MJD'].keys() 00301 else: 00302 print "retdict['data'] has no 'MJD' key." 00303 print "retdict['data'].keys() =", retdict['data'].keys() 00304 else: 00305 print "retdict has no 'data' key." 00306 raise e 00307 00308 for hk in headers: 00309 if retdict.has_key(hk): 00310 if headers[hk].has_key('unit'): 00311 if hk == 'radii': 00312 radii = retdict[hk].split('x') 00313 if len(radii)==1: 00314 a = float(radii[0]) 00315 retdict[hk] = {'unit': headers[hk]['unit'], 'value': (a,a,a)} 00316 retdict['meanrad'] = {'unit': headers[hk]['unit'], 00317 'value': a} 00318 else: 00319 a, b, c = [float(r) for r in radii] 00320 retdict[hk] = {'unit': headers[hk]['unit'], 00321 'value': (a, b, c)} 00322 retdict['meanrad'] = {'unit': headers[hk]['unit'], 00323 'value': mean_radius(a, b, c)} 00324 else: 00325 try: 00326 # meanrad might already have been converted. 00327 if type(retdict[hk]) != dict: 00328 retdict[hk] = {'unit': headers[hk]['unit'], 00329 'value': float(retdict[hk])} 00330 except Exception, e: 00331 print "Error converting header", hk, "to a Quantity." 00332 print "retdict[hk] =", retdict[hk] 00333 raise e 00334 elif hk == 'GeoLong': 00335 long_lat_alt = retdict[hk].split(',') 00336 retdict['GeoLong'] = float(long_lat_alt[0]) 00337 retdict['GeoLat'] = float(long_lat_alt[1]) 00338 retdict['GeoDist'] = float(long_lat_alt[2]) 00339 elif hk == 'dMJD': 00340 retdict[hk] = qa.convert(qa.totime(retdict[hk].replace('minutes', 'min')), 00341 'd')['value'] 00342 elif hk == 'orb_per': 00343 unit = 'h' 00344 retrograde = False 00345 if 'd' in retdict[hk].lower(): 00346 unit = 'd' # Actually this is most common. 00347 if 'r' in retdict[hk].lower(): 00348 retrograde = True 00349 value = get_num_from_str(retdict[hk], 'orbital period') 00350 if value != False: 00351 if retrograde and value > 0.0: 00352 value = -value 00353 retdict[hk] = {'unit': unit, 'value': value} 00354 else: 00355 del retdict[hk] 00356 00357 # The rotation period might depend on the orbital period ("Synchronous"), 00358 # so handle it after all the other headers have been done. 00359 if 'rot_per' in retdict: 00360 rpstr = retdict['rot_per'] 00361 if 'ROTPER' in rpstr: # Asteroid 00362 retdict['rot_per'] = {'unit': 'h', # Always seems to be for asteroids. 00363 'value': get_num_from_str(rpstr, 'rotation period')} 00364 elif 'Synchronous' in rpstr: 00365 retdict['rot_per'] = retdict['orb_per'] 00366 else: # Most likely a planet. 00367 match = re.search(r'(\d+)h\s*(\d+)m\s*([0-9.]+)s', rpstr) 00368 if match: 00369 hms = [float(match.group(i)) for i in range(1, 4)] 00370 retdict['rot_per'] = {'unit': 'h', 00371 'value': hms[0] + (hms[1] + hms[2] / 60.0) / 60.0} 00372 else: 00373 # DON'T include the optional r in hr! qa.totime can't handle it. 00374 try: 00375 match = re.search(r'([-0-9.]+)(?:\s*\+-[0-9.]+)?\s*([dh])', rpstr) 00376 if match: 00377 retdict['rot_per'] = {'unit': match.group(2), 00378 'value': float(match.group(1))} 00379 except: 00380 print "Error parsing the rotation period from" 00381 print rpstr 00382 00383 if retdict['data'].has_key('ang_sep'): 00384 retdict['data']['obs_code'] = {'comment': 'Obscuration code'} 00385 for dk in retdict['data']: 00386 if dk == 'obs_code': 00387 continue 00388 if cols[dk].has_key('unit'): 00389 retdict['data'][dk]['data'] = {'unit': cols[dk]['unit'], 00390 'value': scipy.array([float(s) for s in retdict['data'][dk]['data']])} 00391 if dk == 'RadVel': 00392 # Convert from km/s to AU/d. Blame MeasComet, not me. 00393 retdict['data'][dk]['data']['unit'] = 'AU/d' 00394 kmps_to_AUpd = qa.convert('1km/s', 'AU/d')['value'] 00395 retdict['data'][dk]['data']['value'] *= kmps_to_AUpd 00396 00397 if re.match(r'.*(RA|DEC)$', dk): 00398 retdict['data'][dk] = convert_radec(retdict['data'][dk]) 00399 elif dk == 'MJD': 00400 retdict['data']['MJD'] = datestrs_to_MJDs(retdict['data']['MJD']) 00401 elif dk == 'ang_sep': 00402 angseps = [] 00403 obscodes = [] 00404 for asoc in retdict['data'][dk]['data']: 00405 angsep, obscode = asoc.split('/') 00406 angseps.append(float(angsep)) 00407 obscodes.append(obscode) 00408 retdict['data'][dk]['data'] = {'unit': 'arcseconds', 00409 'value': angseps} 00410 retdict['data']['obs_code']['data'] = obscodes 00411 00412 if len(retdict.get('radii', {'value': []})['value']) == 3 \ 00413 and retdict['data'].has_key('NP_RA') and retdict['data'].has_key('NP_DEC'): 00414 # Do a better mean radius estimate using the actual theta. 00415 retdict['meanrad']['value'] = mean_radius_with_known_theta(retdict) 00416 00417 # To be eventually usable as a MeasComet table, a few more keywords are needed. 00418 retdict['VS_TYPE'] = 'Table of comet/planetary positions' 00419 retdict['VS_VERSION'] = '0003.0001' 00420 if retdict.has_key('VS_CREATE'): 00421 dt = time.strptime(retdict['VS_CREATE'], "%b %d %H:%M:%S %Y") 00422 else: 00423 casalog.post("The ephemeris creation date was not found. Using the current time.", 00424 priority="WARN") 00425 dt = time.gmtime() 00426 retdict['VS_CREATE'] = time.strftime('%Y/%m/%d/%H:%M', dt) 00427 00428 # VS_DATE is required by MeasComet, but it doesn't seem to be actually used. 00429 retdict['VS_DATE'] = time.strftime('%Y/%m/%d/%H:%M', time.gmtime()) 00430 00431 if retdict['data'].has_key('MJD'): 00432 retdict['MJD0'] = retdict['data']['MJD']['value'][0] - retdict['dMJD'] 00433 else: 00434 print "The table will not be usable with me.framecomet because it lacks MJD." 00435 00436 return retdict 00437 00438 def convert_radec(radec_col): 00439 """ 00440 Returns a column of RAs or declinations as strings, radec_col, as a 00441 quantity column. (Unfortunately MeasComet assumes the columns are 00442 Quantities instead of Measures, and uses GeoDist == 0.0 to toggle between 00443 APP and TOPO.) 00444 """ 00445 angstrlist = radec_col['data'] 00446 angq = {} 00447 nrows = len(angstrlist) 00448 00449 if len(angstrlist[0].split()) > 1: 00450 # Prep angstrlist for qa.toangle() 00451 if radec_col['comment'][:len("declination")].lower() == 'declination': 00452 for i in xrange(nrows): 00453 dms = angstrlist[i].replace(' ', 'd', 1) 00454 angstrlist[i] = dms.replace(' ', 'm') + 's' 00455 else: # R.A. 00456 for i in xrange(nrows): 00457 angstrlist[i] = angstrlist[i].replace(' ', ':') 00458 00459 # Do first conversion to get unit. 00460 try: 00461 firstang = qa.toangle(angstrlist[0]) 00462 except Exception, e: 00463 print "Error: Could not convert", angstrlist[0], "to an angle." 00464 raise e 00465 angq['unit'] = firstang['unit'] 00466 angq['value'] = [firstang['value']] 00467 00468 for angstr in angstrlist[1:]: 00469 angq['value'].append(qa.toangle(angstr)['value']) 00470 else: 00471 angq['unit'] = 'deg' # This is an assumption! 00472 angq['value'] = [float(a) for a in angstrlist] 00473 00474 return {'comment': radec_col['comment'], 00475 'data': {'unit': angq['unit'], 00476 'value': scipy.array(angq['value'])}} 00477 00478 def get_num_from_str(fstr, wanted="float"): 00479 """ 00480 Like float(fstr) on steroids, in that it ignores things in fstr that aren't 00481 numbers. Returns False on failure. 00482 00483 wanted: an optional label for the type of number you wanted. 00484 Only used for distinguishing error messages. 00485 00486 Example: 00487 >>> from JPLephem_reader import get_num_from_str 00488 >>> get_num_from_str(' Sidereal rot. period = 58.6462 d ') 00489 58.6462 00490 >>> get_num_from_str('Rotation period = 16.11+-0.01 hr', wanted='rotation period') 00491 16.109999999999999 00492 >>> get_num_from_str('Rotation period = Synchronous', wanted='rotation period') 00493 Could not convert "Rotation period = Synchronous" to a rotation period. 00494 False 00495 """ 00496 match = re.search(r'([-+]?(\d+(\.\d*)?|\d*\.\d+)([eEdD][-+]?\d+)?)', fstr) 00497 if match: 00498 value = float(match.group(1)) 00499 else: 00500 print "Could not convert \"%s\" to a %s." % (fstr, wanted) 00501 value = False 00502 return value 00503 00504 def mean_radius(a, b, c): 00505 """ 00506 Return the average apparent mean radius of an ellipsoid with semiaxes 00507 a >= b >= c. 00508 "average" means average over time naively assuming the pole orientation 00509 is uniformly distributed over the whole sphere, and "apparent mean radius" 00510 means a radius that would give the same area as the apparent disk. 00511 """ 00512 # This is an approximation, but it's not bad. 00513 # The exact equations for going from a, b, c, and the Euler angles to the 00514 # apparent ellipse are given in Drummond et al, Icarus, 1985a. 00515 # It's the integral over the spin phase that I have approximated, so the 00516 # approximation is exact for b == a, and appears to hold well for b << a. 00517 R = 0.5 * c**2 * (1.0 / b**2 + 1.0 / a**2) # The magic ratio. 00518 if R < 0.95: 00519 sqrt1mR = scipy.sqrt(1.0 - R) 00520 # There is fake singularity (RlnR) at R = 0, but it is unlikely to 00521 # be a problem. 00522 try: 00523 Rterm = 0.5 * R * scipy.log((1.0 + sqrt1mR) / (1.0 - sqrt1mR)) / sqrt1mR 00524 except: 00525 Rterm = 0.0 00526 else: 00527 # Use a (rapidly converging) series expansion to avoid a fake 00528 # singularity at R = 1. 00529 Rterm = 1.0 # 0th order 00530 onemR = 1.0 - R 00531 onemRtothei = 1.0 00532 for i in xrange(1, 5): # Start series at 1st order. 00533 onemRtothei *= onemR 00534 Rterm -= onemRtothei / (0.5 + 2.0 * i**2) 00535 avalfabeta = 0.5 * a * b * (1.0 + Rterm) 00536 return scipy.sqrt(avalfabeta) 00537 00538 def mean_radius_with_known_theta(retdict): 00539 """ 00540 Return the average apparent mean radius of an ellipsoid with semiaxes 00541 a >= b >= c (= retdict['radii']['value']). 00542 "average" means average over a rotation period, and "apparent mean radius" 00543 means the radius of a circle with the same area as the apparent disk. 00544 """ 00545 a = retdict['radii']['value'][0] 00546 b2 = retdict['radii']['value'][1]**2 00547 c2 = retdict['radii']['value'][2]**2 00548 onemboa2 = 1.0 - b2 / a**2 00549 units = {} 00550 values = {} 00551 for c in ['RA', 'DEC', 'NP_RA', 'NP_DEC']: 00552 units[c] = retdict['data'][c]['data']['unit'] 00553 values[c] = retdict['data'][c]['data']['value'] 00554 av = 0.0 00555 nrows = len(retdict['data']['RA']['data']['value']) 00556 for i in xrange(nrows): 00557 radec = me.direction('app', {'unit': units['RA'], 'value': values['RA'][i]}, 00558 {'unit': units['DEC'], 'value': values['DEC'][i]}) 00559 np = me.direction('j2000', {'unit': units['NP_RA'], 'value': values['NP_RA'][i]}, 00560 {'unit': units['NP_DEC'], 'value': values['NP_DEC'][i]}) 00561 szeta2 = scipy.sin(qa.convert(me.separation(radec, np), 'rad')['value'])**2 00562 csinz2 = c2 * szeta2 00563 bcosz2 = b2 * (1.0 - szeta2) 00564 bcz2pcsz2 = bcosz2 + csinz2 00565 m = csinz2 * onemboa2 / bcz2pcsz2 00566 av += (scipy.sqrt(bcz2pcsz2) * scipy.special.ellipe(m) - av) / (i + 1.0) 00567 return scipy.sqrt(2.0 * a * av / scipy.pi) 00568 00569 def datestr_to_epoch(datestr): 00570 """ 00571 Given a UT date like "2010-May-01 00:00", returns an epoch measure. 00572 """ 00573 return me.epoch(rf='UTC', v0=qa.totime(datestr)) 00574 00575 def datestrs_to_MJDs(cdsdict): 00576 """ 00577 All of the date strings must have the same reference frame (i.e. UT). 00578 """ 00579 datestrlist = cdsdict['data'] 00580 00581 # Convert to FITS format, otherwise qa.totime() will silently drop the hours. 00582 datestrlist = [d.replace(' ', 'T') for d in datestrlist] 00583 00584 timeq = {} 00585 # Do first conversion to get unit. 00586 firsttime = qa.totime(datestrlist[0]) 00587 timeq['unit'] = firsttime['unit'] 00588 timeq['value'] = [firsttime['value']] 00589 00590 for datestr in datestrlist[1:]: 00591 timeq['value'].append(qa.totime(datestr)['value']) 00592 00593 return {'unit': timeq['unit'], 00594 'value': scipy.array(timeq['value'])} 00595 00596 def construct_tablepath(fmdict, prefix=''): 00597 """ 00598 Construct a suitable pathname for a CASA table made from fmdict, 00599 starting with prefix. prefix can contain a /. 00600 00601 If prefix is not given, it will be set to 00602 "ephem_JPL-Horizons_%s" % fmdict['NAME'] 00603 """ 00604 if not prefix: 00605 prefix = "ephem_JPL-Horizons_%s" % fmdict['NAME'] 00606 return prefix + "_%.0f-%.0f%s%s.tab" % (fmdict['earliest']['m0']['value'], 00607 fmdict['latest']['m0']['value'], 00608 fmdict['latest']['m0']['unit'], 00609 fmdict['latest']['refer']) 00610 00611 def ephem_dict_to_table(fmdict, tablepath='', prefix=''): 00612 """ 00613 Converts a dictionary from readJPLephem() to a CASA table, and attempts to 00614 save it to either to tablepath or a constructed directory name. 00615 Returns whether or not it was successful. 00616 00617 If tablepath is blank and prefix is not given, the table will go to 00618 something like ephem_JPL-Horizons_NAME_EARLIEST-LATESTdUTC.tab. 00619 00620 If tablepath is blank and prefix is given, the table will go to 00621 something like prefix_EARLIEST-LATESTdUTC.tab. prefix can contain a /. 00622 """ 00623 if not tablepath: 00624 tablepath = construct_tablepath(fmdict, prefix) 00625 print "Writing to", tablepath 00626 00627 # safe gaurd from zapping current directory by dict_to_table() 00628 if tablepath=='.' or tablepath=='./' or tablepath.isspace(): 00629 raise Exception, "Invalid tablepath: "+tablepath 00630 retval = True 00631 # keepcolorder=T preserves column ordering in collist below 00632 keepcolorder=False 00633 try: 00634 outdict = fmdict.copy() # Yes, I want a shallow copy. 00635 kws = fmdict.keys() 00636 kws.remove('data') 00637 collist = outdict['data'].keys() 00638 00639 # For cosmetic reasons, encourage a certain order to the columns, i.e. 00640 # start with alphabetical order, 00641 collist.sort() 00642 # but put these ones first, in the listed order (ignore the reverse and 00643 # the pops) if they are present. 00644 #put_these_first = ['MJD', 'RA', 'DEC', 'Rho', 'RadVel', 'NP_RA', 'NP_DEC', 00645 # 'DiskLong', 'DiskLat', 'sl_lon', 'sl_lat', 'r', 00646 # 'ang_sep', 'obs_code'] 00647 put_these_first = ['MJD', 'RA', 'DEC', 'Rho', 'RadVel', 'NP_ang', 'NP_dist', 00648 # 'Obs_Long', 'Obs_Lat', 'Sl_lon', 'Sl_lat', 'r','rdot'] 00649 'DiskLong', 'DiskLat', 'Sl_lon', 'Sl_lat', 'r','rdot'] 00650 # Like l.sort(), reverse() acts on its instance instead of returning a value. 00651 put_these_first.reverse() 00652 for c in put_these_first: 00653 if c in collist: 00654 collist.remove(c) 00655 collist.insert(0, c) 00656 00657 clashing_cols = [c for c in collist if c in kws] 00658 if clashing_cols: 00659 raise ValueError, 'The input dictionary lists ' + ', '.join(clashing_cols) + ' as both keyword(s) and column(s)' 00660 00661 # This promotes the keys in outdict['data'] up one level, and removes 00662 # 'data' as a key of outdict. 00663 outdict.update(outdict.pop('data')) 00664 00665 # This is primarily because MeasComet insists on it, not because it 00666 # ever gets used. Maybe subType should be changed to 'Asteroid', 00667 # 'Moon', or 'Planet', but I'm leaving it at 'Comet' for now. 00668 info = {'readme': 'Derived by JPLephem_reader.py from a JPL-Horizons ephemeris (http://ssd.jpl.nasa.gov/horizons.cgi#top)', 00669 'subType': 'Comet', 'type': 'IERS'} 00670 00671 retval = dict_to_table(outdict, tablepath, kws, collist, info, keepcolorder) 00672 except Exception, e: 00673 print 'Error', e, 'trying to export an ephemeris dict to', tablepath 00674 retval = False 00675 00676 return retval 00677 00678 00679 def jplfiles_to_repository(objs, jpldir='.', jplext='.ephem', 00680 log='null'): 00681 """ 00682 For each Solar System object obj in the list objs, 00683 look for matching JPL-Horizons ASCII files with jplext in jpldir, 00684 read them into python dictionaries, 00685 write the dicts to CASA tables in $CASAROOT/data/ephemerides/JPL-Horizons/, 00686 and check that they can be read by me.framecomet(). 00687 Returns the number of ephemerides processed + readable by me.framecomet. 00688 00689 jpldir and jplext can be glob patterns. 00690 00691 $CASAROOT is derived from $CASAPATH. 00692 00693 Log messages will be directed to log for the duration of this function. 00694 Note that 'null' makes a NullLogSink, so it might be better than /dev/null. 00695 00696 Example: 00697 import recipes.ephemerides.request as jplreq 00698 objs = jplreq.asteroids.keys() + jplreq.planets_and_moons.keys() 00699 jplfiles_to_repository(objs, os.getenv('CASAPATH').split()[0]) 00700 """ 00701 neph = 0 00702 casapath = os.getenv('CASAPATH') 00703 if not casapath: 00704 print "CASAPATH is not set." 00705 return 0 00706 datadir = casapath.split()[0] + '/data/ephemerides/JPL-Horizons' 00707 if not os.path.isdir(datadir): 00708 try: 00709 os.mkdir(datadir) 00710 print "Created", datadir 00711 print "You should probably svn add it." 00712 except Exception, e: 00713 "Error", e, "creating", datadir 00714 return 0 00715 datadir += '/' 00716 00717 #oldlog = casalog.logfile() 00718 # This is needed to stop WARN and above from printing to the console, 00719 # but it permanently severs the logger window. 00720 #casalog.setglobal(True) 00721 #casalog.setlogfile(log) 00722 00723 if jpldir[-1] != '/': 00724 jpldir += '/' 00725 for sob in objs: 00726 capob = sob.capitalize() 00727 lob = sob.lower() 00728 jplfiles = glob(jpldir + lob + jplext) + glob(jpldir + capob + jplext) 00729 for jplfile in jplfiles: 00730 casalog.post('Reading ' + jplfile) 00731 fmdict = readJPLephem(jplfile) 00732 tabpath = construct_tablepath(fmdict, datadir + capob) 00733 ephem_dict_to_table(fmdict, tabpath) 00734 00735 # Check if it is readable by me.framecomet. 00736 epoch = fmdict['earliest'] 00737 epoch['m0']['value'] += 0.5 * (fmdict['latest']['m0']['value'] - 00738 epoch['m0']['value']) 00739 me.doframe(epoch) 00740 if me.framecomet(tabpath): 00741 neph += 1 00742 else: 00743 casalog.post(tabpath + " was not readable by me.framecomet.", 00744 'WARN') 00745 00746 #casalog.setlogfile(oldlog) 00747 00748 return neph
1.8.0
