ison.py

#!/usr/bin/env python2.6
# -*- coding: iso-8859-15 -*-
"""Calculate position for comet

Author: Kelsey Jordahl
Time-stamp: <Sun Feb 26 12:28:17 EST 2012>
License: GPLv3
"""

import ephem
from numpy import pi
from datetime import date, time, datetime


def time_string(et):
    """Output a formatted string for an ephem.Date object
    """
    return ephem.localtime(et).ctime()
    et = et.tuple()
    d = date(et[0], et[1], et[2])
    t = time(et[3], et[4], int(et[5]))
    dt = datetime.combine(d, t)
    return dt.strftime("%d %b %Y %H:%M:%S UTC")


def deg(angle):
    """Return an integer number of degrees, given an angle in radians
    """
    a = round(angle * 180 / pi)
    return int(a)


def main():
    # http://www.minorplanetcenter.net/iau/Ephemerides/Comets/Soft03Cmt.txt
    garradd = ephem.readdb("C/2012 S1 (ISON),h,11/28.7672/2013,61.9802,295.7166,345.5222,1.000005,0.012487,2000,7.0,3.2")

    sun = ephem.Sun()
    moon = ephem.Moon()

    # home location
    home = ephem.Observer()
    # West Orange
    home.lat, home.lon = '40.71', '-74.00'
    home.date = ephem.date(datetime.utcnow())
    garradd.compute(home)
    sun.compute(home)
    sunrise, sunset = (home.next_rising(sun),
                       home.next_setting(sun))
    moon.compute(home)
    moonrise, moonset = (home.next_rising(moon),
                       home.next_setting(moon))
    home.horizon = '-6'
    civil_start, civil_end = (home.next_rising(sun, use_center=True),
                              home.next_setting(sun, use_center=True))
    home.horizon = '-12'
    nautical_start, nautical_end = (home.next_rising(sun, use_center=True),
                                    home.next_setting(sun, use_center=True))
    home.horizon = '-18'
    astro_start, astro_end = (home.next_rising(sun, use_center=True),
                              home.next_setting(sun, use_center=True))
    print time_string(home.date)
    print "Moon phase: %d%% full" % moon.phase
    d = {}
    d[sunset] = " Sunset:                       " + time_string(sunset)
    d[civil_end] = " Civil twilight ends:          " + time_string(civil_end)
    d[nautical_end] = " Nautical twilight ends:       " + time_string(nautical_end)
    d[astro_end] = " Astronomical twilight ends:   " + time_string(astro_end)
    d[astro_start] = " Astronomical twilight starts: " + time_string(astro_start)
    d[nautical_start] = " Nautical twilight starts:     " + time_string(nautical_start)
    d[civil_start] = " Civil twilight starts:        " + time_string(civil_start)
    d[sunrise] = " Sunrise:                      " + time_string(sunrise)
    d[moonrise] = " Moonrise:                     " + time_string(moonrise)
    d[moonset] = " Moonset:                      " + time_string(moonset)

    transit_time = home.next_transit(garradd)
    print u"\nC/2012 S1 (ISON) mag %3.1f in %s: RA %s dec %s" \
        % (garradd.mag, ephem.constellation(garradd)[1], garradd.ra, garradd.dec)
    d[transit_time] = " Comet Garradd transit:        " + time_string(transit_time)
    
    for time in sorted(d.keys()):
        print d[time]
    print ''
    home.date = astro_end
    garradd.compute(home)
    print u"alt %d°, az %d° at end of astronomical twilight" \
        % (deg(garradd.alt), deg(garradd.az))
    home.date = astro_start
    garradd.compute(home)
    print u"alt %d°, az %d° at start of astronomical twilight" \
        % (deg(garradd.alt), deg(garradd.az))
    home.date = transit_time
    garradd.compute(home)
    print u"altitude at transit: %d°" % deg(garradd.alt)

if __name__ == "__main__":
    main()