Permalink
Find file
Fetching contributors…
Cannot retrieve contributors at this time
executable file 405 lines (341 sloc) 14.1 KB
#! /usr/bin/env python
# Predict planetary visibility in the early evening (sunset to midnight),
# and upcoming conjunctions between two or more planets.
# Copyright 2014 Akkana Peck -- share and enjoy under the GPLv2 or later.
import ephem
import math
verbose = False
output_csv = True
# How low can a planet be at sunset or midnight before it's not interesting?
min_alt = 10. * math.pi / 180.
# How close do two bodies have to be to consider it a conjunction?
max_sep = 3.5 * math.pi / 180.
# How little percent illuminated do we need to consider something a crescent?
crescent_percent = 40
# Start and end times for seeing a crescent phase:
crescents = { "Mercury": [ None, None ], "Venus": [ None, None ] }
sun = ephem.Sun()
planets = [
ephem.Moon(),
ephem.Mercury(),
ephem.Venus(),
ephem.Mars(),
ephem.Jupiter(),
ephem.Saturn()
]
planets_up = {}
for planet in planets:
planets_up[planet.name] = None
def datestr(d):
tup = d.tuple()
return "%d/%d/%d" % (tup[0], tup[1], tup[2])
def sepstr(sep):
deg = float(sep) * 180. / math.pi
# if deg < .5:
# return "less than a half a degree (%.2f)" % deg
# if deg < 1.:
# return "less than a degree (%.2f)" % deg
return "%.1f deg" % deg
class ConjunctionPair:
'''A conjunction between a pair of objects'''
def __init__(self, b1, b2, date, sep):
self.bodies = [b1, b2]
self.date = date
self.sep = sep
def __repr__(self):
return "%s: %s and %s, sep %s" % (datestr(self.date), self.bodies[0],
self.bodies[1], sepstr(self.sep))
def __contains__(self, body):
return body in self.bodies
class Conjunction:
'''A collection of ConjunctionPairs which may encompass more
than two bodies and several days.
The list is not guaranteed to be in date (or any other) order.
'''
def __init__(self):
self.bodies = []
self.pairs = []
def __contains__(self, body):
return body in self.bodies
def add(self, body1, body2, date, sep):
self.pairs.append(ConjunctionPair(body1, body2, date, sep))
if body1 not in self.bodies:
self.bodies.append(body1)
if body2 not in self.bodies:
self.bodies.append(body2)
def start_date(self):
date = ephem.date('3000/1/1')
for pair in self.pairs:
if pair.date < date:
date = pair.date
return date
def end_date(self):
date = ephem.date('0001/1/1')
for pair in self.pairs:
if pair.date > date:
date = pair.date
return date
def find_min_seps(self):
return mindate, maxdate, minseps
def andjoin(self, names):
'''Join a list together like a, b, c and d'''
if len(names) == 1:
return names[0]
elif len(names) < 4:
return ', '.join(names[:-1]) + ' and ' + names[-1]
def closeout(self):
'''Time to figure out what we have and print it.'''
# Find the list of minimum separations between each pair.
startdate = ephem.date('3000/1/1')
enddate = ephem.date('0001/1/1')
minseps = []
for i, b1 in enumerate(self.bodies):
for b2 in self.bodies[i+1:]:
minsep = 360 # degrees
closest_date = None
for pair in self.pairs:
if pair.date < startdate:
startdate = pair.date
if pair.date > enddate:
enddate = pair.date
if b1 in pair and b2 in pair:
if pair.sep < minsep:
minsep = pair.sep
closest_date = pair.date
# Not all pairs will be represented. In a triple conjunction,
# the two outer bodies may never get close enough to register
# as a conjunction in their own right.
if minsep < max_sep:
minseps.append((closest_date, minsep, b1, b2))
minseps.sort()
if output_csv:
s = '"Conjunction of ' + self.andjoin(self.bodies) + '",'
s += datestr(startdate) + "," + datestr(enddate) + ",,"
s += "\""
for m in minseps:
s += " %s and %s will be closest on %s (%s)." % \
(m[2], m[3], datestr(m[0]), sepstr(m[1]))
s += "\",,http://upload.wikimedia.org/wikipedia/commons/thumb/4/47/Sachin_Nigam_-_starry_moon_%28by-sa%29.jpg/320px-Sachin_Nigam_-_starry_moon_%28by-sa%29.jpg,240,169,\"<a href='http://commons.wikimedia.org/wiki/File:Sachin_Nigam_-_starry_moon_%28by-sa%29.jpg'>starry moon on Wikimedia Commons</a>\""
print s
else:
print "Conjunction of", self.andjoin(self.bodies),
print "lasts from %s to %s." % (datestr(startdate), datestr(enddate))
for m in minseps:
print " %s and %s are closest on %s (%s)." % \
(m[2], m[3], datestr(m[0]), sepstr(m[1]))
def merge(self, conj):
'''Merge in another Conjunction -- it must be that the two
sets of pairs have bodies in common.
'''
for p in conj.pairs:
self.pairs.append(p)
for body in conj.bodies:
if body not in self.bodies:
self.bodies.append(body)
class ConjunctionList:
'''A collection of Conjunctions -- no bodies should be shared
between any of the conjunctions we contain.
'''
def __init__(self):
self.clist = []
def add(self, b1, b2, date, sep):
for i, c in enumerate(self.clist):
if b1 in c or b2 in c:
c.add(b1, b2, date, sep)
# But what if one of the bodies is already in one of our
# other Conjunctions? In that case, we have to merge.
for cc in self.clist[i+1:]:
if b1 in cc or b2 in cc:
c.merge(cc)
self.clist.delete(cc)
return
# It's new, so just add it
c = Conjunction()
c.add(b1, b2, date, sep)
self.clist.append(c)
def closeout(self):
'''When we have a day with no conjunctions, check the list
and close out any pending conjunctions.
'''
for c in self.clist:
c.closeout()
self.clist = []
oneday = ephem.hour * 24
web_image = {
"Moon" : ("http://upload.wikimedia.org/wikipedia/commons/thumb/5/54/Phase-088.jpg/240px-Phase-088.jpg", '''"<a href='http://commons.wikimedia.org/wiki/User:JayTanner/gallery'>Jay Tanner</a>"''', 240, 240),
"Mercury" : ("../resources/astronomy/mercury.jpg", "", 240, 182),
"Venus" : ("../resources/astronomy/venus.jpg", "", 240, 192),
"Mars" : ("http://imgsrc.hubblesite.org/hu/db/images/hs-2001-24-a-small_web.jpg,200,200", "Hubble Space Telescope", 200, 200),
"Jupiter" : ("http://upload.wikimedia.org/wikipedia/commons/thumb/e/e2/Jupiter.jpg/240px-Jupiter.jpg", '"USGS, JPL and NASA"', 240, 240),
"Saturn" : ("http://upload.wikimedia.org/wikipedia/commons/thumb/b/b4/Saturn_%28planet%29_large.jpg/384px-Saturn_%28planet%29_large.jpg", "Voyager 2", 192, 240)
}
descriptions = {
"Mars": "Mars is visible as a bright, reddish \"star\".",
"Saturn": "Saturn is visible. A small telescope will show its rings.",
"Jupiter": "Jupiter is visible. With binoculars you can see its four brightest moons."
}
def quotecsv(s):
if ',' in s or '"' in s:
return '"' + s.replace('"', '""') + '"'
return s
def finish_planet(p, d):
if not planets_up[p]:
return
if p in descriptions.keys():
if output_csv:
isvis = quotecsv(descriptions[p])
else:
isvis = descriptions[p]
elif p == "Venus" or p == "Mercury":
isvis = p + " is visible in the early evening sky."
else:
isvis = p + " is visible."
# How about crescent info?
if p in crescents.keys():
if crescents[p][0]:
isvis += " A telescope will show a crescent from " \
+ datestr(crescents[p][0])
if crescents[p][1]:
isvis += " to " + datestr(crescents[p][1])
isvis += '.'
crescents[p] = [ None, None ]
if output_csv:
if p != 'Moon':
if web_image[p]:
img = web_image[p][0]
cred = web_image[p][1]
w = web_image[p][2]
h = web_image[p][3]
else:
img = ""
cred = ""
w = ""
h = ""
print "%s,%s,%s,,%s,,%s,%s,%s,%s" % \
(p, datestr(planets_up[p]), datestr(d), isvis,
img, w, h, cred)
else:
print datestr(planets_up[p]), "to", datestr(d), ":", isvis
planets_up[p] = None
def run(start, end, observer, toolate):
'''Find planetary visibility between dates start and end,
for an observer whose location has been set,
between sunset and "toolate" on each date, where toolate is a GMT hour,
e.g. toolate=7 means we'll stop at 0700 GMT or midnight MDT.
'''
d = start
conjunctions = ConjunctionList()
if output_csv:
print 'name,start,end,time,longname,URL,image,image width,image height,image credit'
else:
print "Looking for planetary events between %s and %s:\n" % \
(datestr(d), datestr(end))
def check_if_planet_up(planet, d):
'''If the planet is currently up, do housekeeping to remember
that status, then return True if it's up, False otherwise.
'''
global crescents, planets_up
if planet.alt < min_alt: # planet is not up
return False
if not planets_up[planet.name]:
planets_up[planet.name] = d;
visible_planets.append(planet)
if planet.name not in crescents.keys():
return True
# Is it a crescent? Update its crescent dates.
if planet.phase <= crescent_percent: # It's a crescent now
if not crescents[planet.name][0]:
crescents[planet.name][0] = d
else:
crescents[planet.name][1] = d
return True
while d < end:
observer.date = d
sunset = observer.previous_setting(sun)
# sunrise = observer.next_rising(sun)
# print "Sunset:", sunset, " Sunrise:", sunrise
midnight = list(observer.date.tuple())
midnight[3:6] = [toolate, 0, 0]
midnight = ephem.date(tuple(midnight))
# We have two lists of planets: planets_up and visible_planets.
# planets_up is a dictionary of the time we first saw each planet
# in its current apparition. It's global, and used by finish_planet.
# visible_planets is a list of planets currently visible.
visible_planets = []
for planet in planets:
# A planet is observable this evening (not morning)
# if its altitude at sunset OR its altitude at midnight
# is greater than a threshold, which we'll set at 10 degrees.
observer.date = sunset
planet.compute(observer)
# print planet.name, "alt at sunset:", planet.alt
if not check_if_planet_up(planet, observer.date):
# If it's not up at sunset, try midnight
observer.date = midnight
if observer.date < sunset:
observer.date += oneday
planet.compute(observer)
if not check_if_planet_up(planet, observer.date):
# Planet is not up. Was it up yesterday?
if planets_up[planet.name]:
finish_planet(planet.name, observer.date)
# print datestr(d), "visible planets:", \
# ' '.join([p.name for p in visible_planets])
# print "planets_up:", planets_up
# Done with computing visible_planets.
# Now look for conjunctions, anything closer than 5 degrees.
# Split the difference, use a time halfway between sunset and midnight.
saw_conjunction = False
observer.date = ephem.date((sunset + midnight)/2)
if len(visible_planets) > 1:
for p, planet in enumerate(visible_planets):
for planet2 in visible_planets[p+1:]:
sep = ephem.separation(planet, planet2)
if sep <= max_sep:
# print datestr(observer.date), planet.name, \
# planet2.name, sepstr(sep)
conjunctions.add(planet.name, planet2.name,
observer.date, sep)
saw_conjunction = True
if not saw_conjunction:
conjunctions.closeout()
# Add a day:
d = ephem.date(d + oneday)
for p in visible_planets:
finish_planet(p.name, d)
if __name__ == '__main__':
import sys
if len(sys.argv) > 1 and sys.argv[1] == "-c":
output_csv = True
sys.argv = sys.argv[1:]
else:
output_csv = False
if len(sys.argv) > 1:
start = ephem.date(sys.argv[1])
else:
start = ephem.date('2014/8/15 04:00')
if len(sys.argv) > 2:
end = ephem.date(sys.argv[2])
else:
end = ephem.date('2017/1/1')
# Loop from start date to end date,
# using a time of 10pm MST, which is 4am GMT the following day.
# end = ephem.date('2016/1/1')
# For testing, this spans a Mars/Moon/Venus conjunction:
# d = ephem.date('2015/2/10 04:00')
# end = ephem.date('2015/3/10')
observer = ephem.Observer()
observer.name = "Los Alamos"
observer.lon = '-106.2978'
observer.lat = '35.8911'
observer.elevation = 2286 # meters, though the docs don't actually say
# What hour GMT corresponds to midnight here?
# Note: we're not smart about time zones. This will calculate
# a time based on the time zone offset right now, whether we're
# currently in DST or not.
# And for now we don't even calculate it, just hardwire it.
midnight = 7
try:
run(start, end, observer, midnight)
except KeyboardInterrupt:
print "Interrupted"