conjunctions.py

#! /usr/bin/env python

# Predict planetary visibility in the early evening (sunset to midnight),
# and upcoming conjunctions between two or more planets.
#
# Specifically written to generate the Sky page on the
# PEEC Nature Guides; so it's hardwired to Los Alamos coordinates
# and timezone, and the SQL format uses URLs for images that the
# PEEC sky page will use.
#
# Copyright 2014 Akkana Peck -- share and enjoy under the GPLv2 or later.

from __future__ import print_function

import ephem
import math

verbose = False

# How low can a planet be at sunset or midnight before it's not interesting?
# We'll half it for the moon.
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.

# Half the moon's diameter, for checking occultations:
halfmoon = ephem.degrees('.25')

# How close does a bright planet need to be from the moon to be mentioned?
moon_sep = 25 * math.pi / 180.

# How little % illuminated do we need to consider an inner planet a crescent?
crescent_percent = 40

# Start and end times for seeing a crescent phase:
crescents = { "Mercury": [ None, None ], "Venus": [ None, None ] }

# 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.
timezone = 7

sun = ephem.Sun()

planets = [
    ephem.Moon(),
    ephem.Mercury(),
    ephem.Venus(),
    ephem.Mars(),
    ephem.Jupiter(),
    ephem.Saturn()
    ]

planets_by_name = ["Moon", "Mercury", "Venus", "Mars", "Jupiter", "Saturn"]

planets_up = {}
for planet in planets:
    planets_up[planet.name] = None

saw_conjunction = False
visible_planets = []


def datestr(d):
    # The date may be wrong because of time zones. Convert to our timezone.
    lt = ephem.localtime(d)
    # return lt.strftime("%m/%d/%Y")
    return lt.strftime("%Y-%m-%d")


# For the web page, it's a little more friendly to say something like
# "will be closest on Mar 3" vs. "will be closest on 2017-03-03".
def friendlydate(d):
    lt = ephem.localtime(d)
    return lt.strftime("%b %d")


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


def check_occultation(b1name, b2name, minsep, closest_date, observer):
    """Use a more fine-grained time step to check for a lunar occultation.
    """
    # Unfortunately, by the time we get here we only have a name,
    # not an ephem body. Get the body back:
    b1 = planets[planets_by_name.index(b1name)]
    b2 = planets[planets_by_name.index(b2name)]

    timestep = ephem.hour / 20

    # Go from a day earlier to a day later:
    enddate = closest_date + oneday
    observer.date = closest_date - oneday
    while observer.date < enddate:
        b1.compute(observer)
        b2.compute(observer)
        sep = ephem.separation(b1, b2)
        if sep < minsep:
            minsep = sep
            closest_date = observer.date
        observer.date += timestep

    return minsep, closest_date


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 __repr__(self):
        return "Conjection: " + ', '.join(map(str, 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]
        else:
            return ', '.join(names[:-1]) + ' and ' + names[-1]

    def closeout(self, observer):
        """Time to figure out what we have and print it."""

        if verbose:
            print("closeout", self.start_date(), "-", self.end_date())
            print("  bodies", self.bodies)
            print("  pairs", self.pairs)

        # Find the list of minimum separations between each pair.
        startdate = ephem.date('3000/1/1')
        enddate = ephem.date('0001/1/1')
        minseps = []
        moonclose = ephem.degrees('1')
        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:
                    # First check for lunar occultations.
                    if (b1 == 'Moon' or b2 == 'Moon') and minsep < moonclose:
                        minsep, closest_date = check_occultation(b1, b2, minsep,
                                                                 closest_date,
                                                                 observer)
                    minseps.append((closest_date, minsep, b1, b2))

        minseps.sort()

        def conjstr(b1, b2, sepdate, minsep):
            if b1 == 'Moon' and minsep < halfmoon:
                # It's probably an occultation.
                # Rather than try to predict the time,
                # just narrow it down to the nearest hour.
                sepdate_tuple = list(sepdate.tuple())
                hour = sepdate_tuple[3] - timezone
                if hour > sepdate_tuple[3]:
                    sepdate_tuple[2] -= 1
                sepdate_tuple[3] = hour
                if hour > 12:
                    hourstr = str(hour-12) + " pm"
                else:
                    hourstr = str(hour) + " am"

                seps = "POSSIBLE OCCULTATION around " + hourstr
            else:
                seps = sepstr(minsep)
            return " %s and %s will be closest on %s (%s)." % \
                     (b1, b2, friendlydate(sepdate), seps)

        if output_format == "csv":
            s = '"Conjunction of ' + self.andjoin(self.bodies) + '",'
            s += datestr(startdate) + "," + datestr(enddate) + ",,"
            s += "\""
            for m in minseps:
                s += conjstr(m[2], m[3], m[0], m[1])

            s += "\",,astronomy/starry_moon.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)

        elif output_format == "sql":
            s = "('Conjunction of " + self.andjoin(self.bodies) + "', "
            s += "'astronomy', 'naked eye', "
            s += "'" + datestr(startdate) + "', '" + datestr(enddate) + "', "
            s += "'"
            for m in minseps:
                s += conjstr(m[2], m[3], m[0], m[1])
            s += "', "
            s += "'astronomy/starry_moon.jpg', "
            s += "240, 169, "
            s += "'<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), end=' ')
            print("lasts from %s to %s." % (datestr(startdate), datestr(enddate)))
            for m in minseps:
                print(" ", conjstr(m[2], m[3], m[0], 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 __repr__(self):
        s = "ConjunctionList:"
        for c in self.clist:
            s += "\n    " + str(c)
        return s

    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, observer):
        """When we have a day with no conjunctions, check the list
           and close out any pending conjunctions.
        """
        for c in self.clist:
            c.closeout(observer)
        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),
    "Moon" : ("astronomy/Phase-088.jpg", '''"<a href='http://commons.wikimedia.org/wiki/User:JayTanner/gallery'>Jay Tanner</a>"''', 240, 240),
    "Mercury" : ("astronomy/mercury.jpg", "", 240, 182),
    "Venus" : ("astronomy/venus.jpg", "", 240, 192),
    # "Mars" : ("http://imgsrc.hubblesite.org/hu/db/images/hs-2001-24-a-small_web.jpg", "Hubble Space Telescope", 200, 200),
    "Mars" : ("astronomy/mars.jpg", "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),
    "Jupiter" : ("astronomy/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)
    "Saturn" : ("astronomy/saturn.jpg", "Voyager 2", 182, 240)
}

descriptions = {
    "Mars": " as a bright, reddish \"star\".",
    "Saturn": ". A small telescope will show its rings.",
    "Jupiter": ". With binoculars you can see its four brightest moons."
    }

def quotecsv(s):
    if ',' in s or '"' in s:
        return '"' + s.replace('"', '""') + '"'
    return s

def escape_singlequotes(s):
    return s.replace("'", "\\'")

def finish_planet(p, d, observer, output_format):
    if not planets_up[p]:
        return

    # Morning or evening? For that, we need to get a handle on
    # the ephem body object again.
    body = planets[planets_by_name.index(p)]
    observer.date = d
    transit = observer.previous_transit(body) - ephem.hour * timezone
    transit = list(ephem.Date(transit).tuple())
    if transit[3] < 3 or transit[3] > 12:
        when = "evening"
    else:
        when = "morning"

    if p == "Venus" or p == "Mercury":
        if when == "evening":
            isvis = p + " is visible in the early evening sky."
        else:
            isvis = p + " is visible in the morning sky."

    elif p in list(descriptions.keys()):
        desc = p + " is visible in the " + when + descriptions[p]
        if output_format == "csv":
            isvis = quotecsv(desc)
        elif output_format == "sql":
            isvis = escape_singlequotes(desc)
        else:
            isvis = desc

    elif p == "Moon":
        isvis = "The moon is visible in the " + when + "."

    else:
        isvis = p + " is visible in the " + when + "."

    # How about crescent info?
    if p in list(crescents.keys()):
        if crescents[p][0]:
            isvis += " A telescope will show a crescent from " \
                     + friendlydate(crescents[p][0])
            if crescents[p][1]:
                isvis += " to " + friendlydate(crescents[p][1])
            isvis += '.'
        crescents[p] = [ None, None ]

    if output_format == "csv" or output_format == "sql":
        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 = ""

            if output_format == "csv":
                print("%s,%s,%s,,%s,,%s,%s,%s,%s" % \
                      (p, datestr(planets_up[p]), datestr(d), isvis,
                       img, w, h, cred))
            else:
                print("('%s', 'astronomy', 'naked eye', '%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, output_format):
    """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.
       toolate==None means look for anything between sunset and sunries.
    """
    global visible_planets, saw_conjunction

    d = start
    conjunctions = ConjunctionList()

    if output_format == "csv":
        print('name,start,end,time,longname,URL,'
              'image,image width,image height,image credit')
    elif output_format == "sql":
        print("INSERT INTO peecnatu_guides_dev.astronomy(common_name, "
              "guide_group, visibility, start, end, comments, image, "
              "image_width, image_height, image_credit) VALUES")
    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 up on the given date, do housekeeping to remember
           that status, then return True if it's up, False otherwise.
           The date passed in is just a date; we will try different times
           on that date, including the immediately preceding sunset
           and a "toolate" hour of the night.
        """
        global crescents, planets_up, visible_planets, saw_conjunction

        # The moon is easy to see, so allow it half the alt of anything else.
        if planet.name == "Moon":
            if planet.alt < min_alt/2:   # moon isn't up
                return False
        elif planet.alt < min_alt:       # planet is not up
            return False

        # Planet is up.
        if not planets_up[planet.name]:
            planets_up[planet.name] = d;
        visible_planets.append(planet)

        if planet.name not in list(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

    # Loop over the days in the time range:
    while d < end:

        # Set d to mid-day sometime, definitely before sunset
        midday = list(d.tuple())
        midday[3:6] = [12 - timezone, 0, 0]
        if midday[3] < 0:
            midday[3] += 12
        observer.date = ephem.date(tuple(midday))
        sunset = observer.next_setting(sun)

        # Stop at a fixed hour of the evening?
        if toolate:
            latenight = list(observer.date.tuple())
            latenight[3:6] = [toolate + 24 - timezone, 0, 0]
            latenight = ephem.date(tuple(latenight))
            if latenight < sunset:
                latenight += oneday

        # Stop at sunrise
        else:
            observer.date = sunset
            latenight = observer.next_rising(sun)

        if verbose:
            print("\n***", d, "from", sunset, "to", latenight)

        # 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 late-night
            # is greater than a visible_threshold
            observer.date = sunset
            planet.compute(observer)
            if verbose:
                print(observer.date, planet, "alt", planet.alt)
            if not check_if_planet_up(planet, observer.date):
                # If it's not up at sunset, try latenight
                observer.date = latenight
                if observer.date < sunset:
                    observer.date += oneday
                planet.compute(observer)
                if verbose:
                    print("  ", observer.date, planet, "alt", planet.alt)
                if not check_if_planet_up(planet, observer.date):
                    # Planet is not up. Was it up yesterday?
                    if planets_up[planet.name] and planet.name != "Moon":
                        finish_planet(planet.name, observer.date,
                                      observer, output_format)

        # print()
        # 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 latenight.
        saw_conjunction = False
        observer.date = ephem.date((sunset + latenight)/2)
        moon = planets[0]
        if len(visible_planets) > 1:
            for p, planet in enumerate(visible_planets):
                for planet2 in visible_planets[p+1:]:
                    sep = ephem.separation(planet, planet2)
                    # print(observer.date, "moon -", planet2.name, sep)
                    if sep <= max_sep:
                        # print (datestr(observer.date), planet.name,
                        #        planet2.name, sepstr(sep))
                        if verbose:
                            print("adding sep", planet.name, planet2.name,
                                  observer.date, sep)
                        conjunctions.add(planet.name, planet2.name,
                                         observer.date, sep)
                        saw_conjunction = True
                    elif planet == moon and sep <= moon_sep:
                        if verbose:
                            print("adding moon sep", planet.name, planet2.name,
                                  observer.date, sep)
                        conjunctions.add(planet.name, planet2.name,
                                         observer.date, sep)
                        saw_conjunction = True

        if not saw_conjunction:
            conjunctions.closeout(observer)

        # Add a day:
        d = ephem.date(d + oneday)

    if saw_conjunction:
        conjunctions.closeout()
    for p in visible_planets:
        if p.name != "Moon":
            finish_planet(p.name, d, observer, output_format)


def moon_phases(start, end, output_format):
    d = ephem.previous_full_moon(start)
    d = ephem.previous_full_moon(d)
    print("Starting from", d)

    def output_moon_phase(d, phasename, img, attr):
        if output_format == "sql":
            print("('%s', 'astronomy', 'naked eye', '%s', '%s', '%s', "
                  "'%s', 240, 240, '%s')," % (phasename + " moon",
                                              datestr(d), datestr(d),
                                              phasename + " moon",
                                              img, attr))
        else:
            print(datestr(d), ":", phasename + " moon")

    while d <= end:
        d = ephem.next_first_quarter_moon(d)
        output_moon_phase(d, "First quarter", 'astronomy/Phase-088.jpg',
                          '<a href=\"http://commons.wikimedia.org/wiki/'
                          'File:Phase-088.jpg\">Jay Tanner</a>')
        d = ephem.next_full_moon(d)
        output_moon_phase(d, "Full", 'astronomy/Phase-180.jpg',
                          '<a href=\"http://commons.wikimedia.org/wiki/'
                          'File:Phase-180.jpg\">Jay Tanner</a>')
        d = ephem.next_last_quarter_moon(d)
        output_moon_phase(d, "Last quarter", 'astronomy/Phase-270.jpg',
                          '<a href=\"http://commons.wikimedia.org/wiki/'
                          'File:Phase-270.jpg\">Jay Tanner</a>')
        d = ephem.next_new_moon(d)
        output_moon_phase(d, "New", 'astronomy/New_Moon.jpg',
                          '<a href="https://commons.wikimedia.org/wiki/'
                          'File:New_Moon.jpg">QuimGil</a>')


if __name__ == '__main__':
    import sys, os

    if len(sys.argv) > 1 and (sys.argv[1] == "-h" or sys.argv[1] == "--help"):
        print("Usage: %s [-c|-s] start_date end_date"
              % os.path.basename(sys.argv[0]))
        print("  -c: CSV output")
        print("  -s: SQL output")
        print("  Otherwise output will be text")
        print("Dates in format yyyy/mm/dd [hh:mm] (PyEphem's preference)")
        sys.exit(0)

    if len(sys.argv) > 1 and sys.argv[1] == "-v":
        verbose = True
        sys.argv = sys.argv[1:]

    if len(sys.argv) > 1 and sys.argv[1] == "-c":
        output_format = "csv"
        sys.argv = sys.argv[1:]
    elif len(sys.argv) > 1 and sys.argv[1] == "-s":
        output_format = "sql"
        sys.argv = sys.argv[1:]
    else:
        output_format = "text"

    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

    # How late is too late to observe? E.g. 1 for 1 am, 23 for 11 pm.
    # None means we're interested in all events from sunset to sunrise.
    toolate = None

    try:
        run(start, end, observer, toolate, output_format)
        moon_phases(start, end, output_format)

    except KeyboardInterrupt:
        print("Interrupted")