Python script to compute optimal Milky Way photography dates and times

compute_milky_way.py

#!/usr/bin/env python2.7

from datetime import date, datetime
import ephem

mylocation = ephem.Observer()
mylocation.lat, mylocation.lon = '39.039502', '-77.486371'

min_alt = ephem.degrees('10:00')

milkyway_visual = ephem.readdb("Milky Way Visual Center,f|M|F7,18:55:03.06,-30:28:42.30,3.00,2000")
milkyway_sagastar = ephem.readdb("Milky Way Black Hole (Sgr A*),f|M|F7,17:45:40.0409,-29:00:28.118,25.00,2000")
milk = milkyway_sagastar
sun = ephem.Sun()
moon = ephem.Moon()

astrosunset_offset = 2 * ephem.hour

# This is around noon EST
start_date = ephem.Date('2016/1/1 16:00:00')
end_date = ephem.Date('2017/1/1 16:00:00')

cur_date = start_date

while cur_date < end_date:
    cur_date = ephem.Date(cur_date + (ephem.hour * 24))
    next_date = cur_date + (ephem.hour * 24)
    mylocation.date = cur_date

    sunset = mylocation.next_setting(sun)
    sunrise = mylocation.next_rising(sun)
    astrosunset = sunset + astrosunset_offset
    milkset = mylocation.next_setting(milk)
    milkrise = mylocation.next_rising(milk)

    if milkset < astrosunset or milkrise > sunrise:
        #print "Milky Way is down before sunset %s" % (cur_date)
        continue

    moonset = mylocation.next_setting(moon)
    moonrise = mylocation.next_rising(moon)

    # Check 10-minute intervals
    interval = ephem.hour / 6
    fine_date = sunset - interval
    target_dates = []
    #print "=>> Starting day at %s [%s]" % (ephem.Date(fine_date), cur_date)

    while fine_date < sunrise:
        fine_date = fine_date + interval
        mylocation.date = ephem.Date(fine_date)
        moon.compute(mylocation)
        if moon.alt > 0:
            #print "Moon is still up %s" % (ephem.Date(fine_date))
            continue

        milk.compute(mylocation)
        if milk.alt < min_alt:
            #print "Milky Way is still below horizon (%s) %s [moon: %s]" % (str(milk.alt), ephem.Date(fine_date), str(moon.alt).split(":")[0])
            continue

        nice_date = ephem.localtime(ephem.Date(fine_date)).strftime('%Y-%m-%d %H:%M:%S')
        nice_time = ephem.localtime(ephem.Date(fine_date)).strftime('%H:%M:%S')
        nice_alt = str(milk.alt).split(":")[0]
        nice_az = str(milk.az).split(":")[0]
        #print "Adding date %s" % ephem.Date(fine_date)
        target_dates.append( {'date': ephem.localtime(ephem.Date(fine_date)), 'nice_date': nice_date, 'nice_time': nice_time, 'alt': nice_alt, 'az': nice_az, 'moon': str(moon.alt).split(":")[0]} )

    if len(target_dates) == 0:
        continue

    if len(target_dates) == 1:
        first, last, apex = target_dates[0], target_dates[0], target_dates[0]
    else:
        first, last, apex = target_dates[0], target_dates[len(target_dates)-1], target_dates[(len(target_dates)/2)]


    print "%s - %s   [apex: Alt: %s deg, Az: %s deg at %s]" % (first['nice_date'], last['nice_time'], apex['alt'], apex['az'], apex['nice_time'])