|
casa
$Rev:20696$
|
|
casa
$Rev:20696$
|
00001 """ 00002 Utilities for having JPL-Horizons ephemerides mailed to you. 00003 See JPL_ephem_reader.py for doing something with them. 00004 00005 Examples: 00006 00007 import recipes.ephemerides.request as jplreq 00008 00009 # I recommend you not ask for more than ~18 months of anything with 00010 # date_incr ~ 1h, because the result would be split into multiple 00011 # emails which you would have to stitch together. 00012 00013 for thing in jplreq.asteroids.keys() + jplreq.planets_and_moons.keys(): 00014 jplreq.request_from_JPL(thing, '2012-12-31') 00015 00016 # A trick to avoid fast moving objects: 00017 for thing in jplreq.asteroids.keys() + jplreq.planets_and_moons.keys(): 00018 if thing not in jplreq.default_date_incrs: 00019 jplreq.request_from_JPL(thing, '2012-12-31') 00020 """ 00021 00022 import os 00023 import re 00024 import smtplib 00025 import socket 00026 import time 00027 from email.mime.text import MIMEText 00028 00029 # Maps from object names to numbers that JPL-Horizons will recognize without 00030 # fussing around with choosing between barycenters, substrings, etc.. 00031 # Use lower case keys. 00032 # Do not use keys for asteroids that are also in planets_and_moons. (The IAU 00033 # might enforce this anyway.) 00034 # 00035 # They are numbered by order of discovery, which is not quite the same as size. 00036 00037 asteroids = {'ceres': 1, 00038 'pallas': 2, 00039 'juno': 3, # Large crater and temperature changes. 00040 'vesta': 4, 00041 'astraea': 5, 00042 'hygiea': 10, # Careful with the spelling. It used to be 00043 # Hygeia, and it is named after the Greek 00044 # goddess Hygieia (or Hygeia). Also, it is 00045 # fairly oblate and eccentric. 00046 'parthenope': 11, 00047 'victoria': 12, 00048 'eunomia': 15, 00049 'euphrosyne': 31, 00050 '52 europa': 52, 00051 'cybele': 65, 00052 'sylvia': 87, # Has two moons. 00053 'davida': 511, 00054 'interamnia': 704} 00055 00056 planets_and_moons = {'mercury': 199, 00057 'venus': 299, 00058 'mars': 499, 00059 'phobos': 401, 00060 'deimos': 402, 00061 'jupiter': 599, 00062 'io': 501, 00063 'europa': 502, 00064 'ganymede': 503, 00065 'callisto': 504, 00066 'saturn': 699, 00067 'mimas': 601, 00068 'enceladus': 602, 00069 'tethys': 603, 00070 'dione': 604, 00071 'rhea': 605, 00072 'titan': 606, 00073 'hyperion': 607, 00074 'iapetus': 608, 00075 'phoebe': 609, 00076 'janus': 610, 00077 'epimetheus': 611, 00078 'helene': 612, 00079 'telesto': 613, 00080 'calypso': 614, 00081 'atlas': 615, 00082 'prometheus': 616, 00083 'pandora': 617, 00084 'pan': 618, 00085 # I've been to Ymir, so I have a soft spot for it, but it 00086 # has an unknown radius (2010). 00087 #'ymir': 619, 00088 'uranus': 799, 00089 'ariel': 701, 00090 'umbriel': 702, 00091 'titania': 703, 00092 'oberon': 704, 00093 'miranda': 705, 00094 'cordelia': 706, 00095 'ophelia': 707, 00096 'bianca': 708, 00097 'cressida': 709, 00098 'desdemona': 710, 00099 'juliet': 711, 00100 'portia': 712, 00101 'rosalind': 713, 00102 'belinda': 714, 00103 'puck': 715, 00104 # 'caliban': 716, Uncertain radius, 2010. 00105 # 'sycorax': 717, Uncertain radius, 2010. 00106 # 'prospero': 718, Unknown radius, 2010 00107 # 'setebos': 719, Unknown radius, 2010 00108 # 'stephano': 720, Unknown radius, 2010 00109 # 'trinculo': 721, Unknown radius, 2010 00110 # 'francisco': 722, " 00111 # 'margaret': 723, Unknown radius, 2010 00112 # 'ferdinand': 724, Unknown radius, 2010 00113 # 'perdita': 725, Unknown radius, 2010 00114 # 'mab': 726, Unknown radius, 2010 00115 # 'cupid': 727, " 00116 'neptune': 899, 00117 'triton': 801, 00118 'nereid': 802, 00119 'naiad': 803, 00120 'thalassa': 804, 00121 'despina': 805, 00122 'galatea': 806, 00123 'larissa': 807, 00124 'proteus': 808, 00125 'pluto': 999, # It's still a planet in this sense. 00126 'charon': 901 00127 # 'nix': 902 Unknown radius, 2010 00128 # 'hydra': 903 Unknown radius, 2010 00129 } 00130 00131 should_have_orientation = ['mars', 'deimos', 'phobos', 'vesta', 'jupiter', 'io', 00132 'janus', 'enceladus', 'mimas', 'iapetus', 00133 'phoebe', 'tethys', 'uranus', 'ariel', 'miranda', 00134 'neptune'] 00135 should_have_sublong = ['mars', 'deimos', 'phobos', 'jupiter', 'io', 00136 'janus', 'enceladus', 'phoebe', 'mimas', 'tethys', 00137 'neptune'] 00138 00139 # Getting positions once a day is not enough for many moons, if the position of 00140 # the moon relative to its primary will be needed. Note that a maximum 00141 # suitable increment is imposed by Earth's motion. 00142 default_date_incrs = { 00143 'default': "1 d", # The default default. 00144 'ariel': '0.5d', 00145 'cordelia': '0.05d', 00146 'deimos': '0.25d', 00147 'dione': '0.5d', 00148 'enceladus': '0.25d', 00149 'io': '0.25d', 00150 'janus': '0.1d', 00151 'mimas': '0.2d', 00152 'miranda': '0.25 d', 00153 'phobos': '0.05d', 00154 'tethys': '0.4d' 00155 } 00156 00157 def request_from_JPL(objnam, enddate, 00158 startdate=None, 00159 date_incr=None, 00160 get_axis_orientation=None, 00161 get_axis_ang_orientation=None, 00162 get_sub_long=None, 00163 obsloc="", 00164 return_address=None, 00165 mailserver=None, 00166 use_apparent=True, 00167 get_sep=None): 00168 """ 00169 Request an ASCII ephemeris table from JPL-Horizons for a Solar System 00170 object. If all goes well it should arrive by email in a few minutes to 00171 an hour. (The return value from this function is whether or not it sent 00172 the request.) 00173 00174 All but the first two parameters have hopefully sensible defaults: 00175 objnam: 00176 The name of the object (case-insensitive). It will be used to refer to 00177 specifically its center, as opposed to other possible locations in the 00178 vicinity of the object. For example, if objnam ="Mars", it will choose 00179 Mars, not the Mars barycenter or the Mars Reconnaissance Orbiter. 00180 enddate: 00181 The date that the ephemeris should end on. 00182 It can be an epoch measure or string (yyyy-mm-dd, assumes UT). 00183 startdate: 00184 Defaults to today, but it can be specified like enddate. 00185 date_incr: 00186 The increment between dates in the ephemeris. casapy's setjy 00187 task and me tool automatically interpolate. It can be a (time) quantity 00188 or a string (which will be interpreted as if it were a quantity). 00189 00190 Unlike the JPL email interface, this does not need it to be an integer 00191 number of time units. request_from_JPL() will do its best to convert 00192 it to fit JPL's required format. 00193 00194 Default: 1 Earth day. 00195 get_axis_orientation: 00196 Request the orientation of the object's polar axis relative to the line 00197 of sight. This is needed (along with the flattening) if treating the 00198 disk as an ellipse, but it is often unavailable. 00199 True or False 00200 Defaults to whether or not objnam is in should_have_orientation. 00201 get_axis_ang_orientation: 00202 Request the angular orientation (position angle and angular distance from 00203 sub-observer point) of the object's polar axis relative to the line 00204 of sight. 00205 True or False (by default it is included) 00206 get_sub_long: 00207 Request the planetographic (geodetic) longitudes and latitudes of the 00208 subobserver and sub-Solar points. Only needed if the object has 00209 significant known surface features. 00210 True or False 00211 Defaults to whether or not objnam is in should_have_sublong. 00212 obsloc: 00213 Observatory name, used to get topocentric coordinates. 00214 Obviously not all observatories are recognized. 00215 Default: "" (geocentric) 00216 return_address: 00217 The email address that the ephemeris will be sent to. 00218 Default: <username>@<domainname>. 00219 mailserver: 00220 The computer at _your_ end to send the mail from. 00221 Default: a semi-intelligent guess. 00222 use_apparent: 00223 Get the apparent instead of J2000 RA and Dec. No refraction by Earth's 00224 atmosphere will be applied; MeasComet assumes apparent directions and 00225 JPL_ephem_reader would be confused if both apparent and J2000 00226 directions were present. 00227 Default: True 00228 get_sep: 00229 Get the angular separation from the primary, and whether it is 00230 transiting, in eclipse, etc.. This only makes sense for moons and does 00231 not guarantee that nothing else (like Earth, Luna, a bright extrasolar 00232 object) is in the line of sight! 00233 Default: True if it is in the moons list, False otherwise. 00234 """ 00235 lobjnam = objnam.lower() 00236 00237 # Handle defaults 00238 if get_axis_orientation == None: # remember False is valid. 00239 if lobjnam in should_have_orientation: 00240 get_axis_orientation = True 00241 else: 00242 get_axis_orientation = False 00243 00244 if get_sub_long == None: # remember False is valid. 00245 if lobjnam in should_have_sublong: 00246 get_sub_long = True 00247 else: 00248 get_sub_long = False 00249 00250 if not return_address: 00251 fqdn = socket.getfqdn() 00252 00253 # Only use the top two levels, i.e. eso.org and nrao.edu, not 00254 # (faraday.)cv.nrao.edu. 00255 domain = '.'.join(fqdn.split('.')[-2:]) 00256 00257 return_address = os.getlogin() + '@' + domain 00258 00259 if not mailserver: 00260 try: 00261 #mailserver = socket.getfqdn(socket.gethostbyname('mail')) 00262 mailserver = socket.getfqdn(socket.gethostbyname('smtp')) 00263 except socket.gaierror: 00264 print "Could not find a mailserver." 00265 return False 00266 00267 if not startdate: 00268 syr, smon, s_d, s_h, smin, s_s, swday, syday, sisdst = time.gmtime() 00269 startdate = "%d-%02d-%02d" % (syr, smon, s_d) 00270 00271 if not date_incr: 00272 date_incr = default_date_incrs.get(lobjnam, 00273 default_date_incrs['default']) 00274 00275 if get_sep == None: 00276 get_sep = (planets_and_moons.get(lobjnam, 99) % 100) < 99 00277 00278 # Get to work. 00279 if lobjnam in asteroids: 00280 objnum = str(asteroids[lobjnam]) + ';' 00281 elif lobjnam in planets_and_moons: 00282 objnum = str(planets_and_moons[lobjnam]) 00283 else: 00284 print "The JPL object number for", objnam, "is not known. Try looking it up at" 00285 print 'http://ssd.jpl.nasa.gov/horizons.cgi?s_body=1#top and adding it.' 00286 return False 00287 00288 if obsloc and obsloc.lower() != 'geocentric': 00289 print "Topocentric coordinates are not yet supported by this script." 00290 print "Defaulting to geocentric." 00291 # to set site coordinates, 00292 # CENTER=coord@500 00293 # COORD_TYPE='GEODETIC' 00294 # SITE_COORD='E.lon,lat,height' (in deg and km) 00295 # e.g for ALMA 00296 # SITE_COORD='-67.7549290,-23.022886,5.05680' 00297 00298 # 500@399: 399=earth 500=body center (==g@399==geo) 00299 center = '500@399' 00300 00301 # quantities = [2, 10, 12, 19, 20, 24] 00302 # Original default set of quantities: apparent RA/DEC, Illum. frac, ang. separation between 00303 # non-lunar target and the center of primary body, helio range and range rate, observer 00304 # range and range rate, S-T-O 00305 # Bryan's request 00306 # [1,14,15,17,19,20,24] 00307 # extra: 10 00308 # need : 17 00309 # new default quantities 00310 quantities = [2, 12, 17, 19, 20, 24] 00311 if not use_apparent: 00312 quantities[0] = 1 00313 if get_axis_orientation: 00314 quantities.append(32) 00315 if get_sub_long: 00316 quantities.append(14) 00317 quantities.append(15) 00318 if not get_sep: 00319 quantities.remove(12) 00320 if not get_axis_ang_orientation: 00321 quantities.remove(17) 00322 print "Retrieved quantity code list=",quantities 00323 00324 # It seems that STEP_SIZE must be an integer, but the unit can be changed 00325 # to hours or minutes. 00326 match = re.match(r'([0-9.]+)\s*([dhm])', date_incr) 00327 n_time_units = float(match.group(1)) 00328 time_unit = match.group(2) 00329 if n_time_units < 1.0: 00330 if time_unit == 'd': 00331 n_time_units *= 24.0 00332 time_unit = 'h' 00333 if time_unit == 'h' and n_time_units < 1.0: # Note fallthrough. 00334 n_time_units *= 60.0 00335 time_unit = 'm' 00336 if n_time_units < 1.0: # Uh oh. 00337 print date_incr, "is an odd request for a date increment." 00338 print "Please change it or make your request manually." 00339 return False 00340 print "Translating date_incr from", date_incr, 00341 date_incr = "%.0f %s" % (n_time_units, time_unit) 00342 print "to", date_incr 00343 00344 instructions = "\n".join(["!$$SOF", 00345 "COMMAND= '%s'" % objnum, 00346 'CENTER= ' + center, 00347 "MAKE_EPHEM= 'YES'", 00348 "TABLE_TYPE= 'OBSERVER'", 00349 "START_TIME= '%s'" % startdate, 00350 "STOP_TIME= '%s'" % enddate, 00351 "STEP_SIZE= '%s'" % date_incr, 00352 "CAL_FORMAT= 'CAL'", #date format(CAL,JD,or BOTH) 00353 "TIME_DIGITS= 'MINUTES'", 00354 "ANG_FORMAT= 'DEG'", 00355 "OUT_UNITS= 'KM-S'", 00356 "RANGE_UNITS= 'AU'", 00357 "APPARENT= 'AIRLESS'", # for apparent position 00358 "SOLAR_ELONG= '0,180'", 00359 "SUPPRESS_RANGE_RATE= 'NO'", 00360 "SKIP_DAYLT= 'NO'", 00361 "EXTRA_PREC= 'NO'", 00362 "R_T_S_ONLY= 'NO'", 00363 "REF_SYSTEM= 'J2000'", 00364 "CSV_FORMAT= 'NO'", 00365 "OBJ_DATA= 'YES'", 00366 "QUANTITIES= '%s'" % ','.join([str(q) for q in quantities]), 00367 '!$$EOF']) 00368 00369 # Set up a MIMEText object (it's a dictionary) 00370 msg = MIMEText(instructions) 00371 00372 msg['To'] = "horizons@ssd.jpl.nasa.gov" 00373 msg['Subject'] = 'JOB' 00374 msg['From'] = return_address 00375 msg['Reply-to'] = return_address 00376 00377 # Establish an SMTP object and connect to the mail server 00378 s = smtplib.SMTP() 00379 s.connect(mailserver) 00380 00381 # Send the email - real from, real to, extra headers and content ... 00382 s.sendmail(return_address, msg['To'], msg.as_string()) 00383 s.close() 00384 00385 return True 00386 00387 def list_moons(): 00388 """ 00389 List planets_and_moons in a more organized way. 00390 """ 00391 # Gather the moons by planet number. 00392 planets = {} 00393 moons = {} 00394 for lcname in planets_and_moons: 00395 num = planets_and_moons[lcname] 00396 planet = num / 100 00397 if num % 100 == 99: 00398 planets[planet] = lcname.title() 00399 else: 00400 if not moons.has_key(planet): 00401 moons[planet] = {} 00402 moons[planet][num % 100] = lcname.title() 00403 00404 #print "planets =", planets 00405 #print "moons:" 00406 #for p in planets: 00407 # print planets[p] 00408 # print " ", moons.get(p, "None") 00409 00410 # For formatting the output table, find the column widths, 00411 # and maximum number of moons per planet. 00412 maxmoons = max([len(moons.get(p, '')) for p in planets]) 00413 maxwidths = {} 00414 for planet in planets: 00415 if moons.has_key(planet): 00416 maxwidths[planet] = max([len(m) for m in moons[planet].values()]) 00417 else: 00418 maxwidths[planet] = 0 00419 if len(planets[planet]) > maxwidths[planet]: 00420 maxwidths[planet] = len(planets[planet]) 00421 00422 # Set up the table columns. 00423 plannums = planets.keys() 00424 plannums.sort() 00425 sortedmoons = {} 00426 formstr = '' 00427 hrule = '' 00428 for p in plannums: 00429 formstr += '| %-' + str(maxwidths[p]) + 's ' 00430 if p == 1: 00431 hrule += '|' 00432 else: 00433 hrule += '+' 00434 hrule += '-' * (maxwidths[p] + 2) 00435 moonkeys = moons.get(p, {}).keys() 00436 moonkeys.sort() 00437 sortedmoons[p] = {} 00438 for row in xrange(len(moonkeys)): 00439 sortedmoons[p][row] = moons[p][moonkeys[row]] 00440 formstr += '|' 00441 hrule += '|' 00442 00443 print formstr % tuple([planets[p] for p in plannums]) 00444 print hrule 00445 for row in xrange(maxmoons): 00446 print formstr % tuple([sortedmoons[p].get(row, '') for p in plannums]) 00447 00448 def list_asteroids(): 00449 """ 00450 Like list_moons, but list the asteroids by their numbers 00451 (= order of discovery, ~ albedo * size) 00452 """ 00453 astnums = asteroids.values() 00454 astnums.sort() 00455 invast = {} 00456 for a in asteroids: 00457 invast[asteroids[a]] = a.title() 00458 for n in astnums: 00459 print "%3d %s" % (n, invast[n])
1.8.0
