#!/usr/bin/env python

# Predict when the moon (or another body) will be at a specified

# altitude and azimuth during a specified time window

# over the course of a year.

# For instance: figure out when you can take a photo of the full

# moon rising over Lick Observatory; or figure out when the full

# moon shining in your skylight might keep you awake.

# Copyright 2015 by Akkana Peck. Share and enjoy under the GPL v2 or later.

import ephem

import sys, os

import datetime

# Define your own observer parameters here:

observer = ephem.Observer()

observer.name = "White Rock"

observer.lon = '-106.22'

observer.lat = '35.82'

observer.horizon = ephem.degrees(7.*ephem.pi/180.)

observer.elevation = 1980 # meters, though the docs don't actually say

DEGREES = 180. / ephem.pi

def discont_range(start, end, max):

'''Return a discontinuous range.

Like range(start, end) except that if start > end, will "loop around"

e.g. discont_range(22, 2, 24) will return [22, 23, 0, 1].

'''

start = start%max

end = end%max

if start <= end:

return range(start, end)

return range(start, max) + range(0, end)

def when_at_position(body, observer, targetalt, targetaz,

starttime, endtime,

slop=5., minphase=0, maxphase=100):

'''When will body be at the alt-az position during the time window,

in the year following the date of the observer passed in

(which defaults to today's date.)

Args:

body (ephem.Body) -- The celestial body to be calculated.

observer(ephem.Observer) -- the observing position.

targetalt (float) -- Altitude in decimal degrees.

targetaz (float) -- Azimuth in decimal degrees.

starttime(int) -- Earliest hour of day to check, GMT,

on 24-hour clock, e.g. 18 for 6pm.

endtime(int) -- Latest hour of day to check, GMT.

Will be ignored if starttime is a string.

slop (Optional) -- How much slop to allow in the alt/az

positions each way (float, decimal degrees).

minphase (Optional, int) -- What phase are we interested in (% illum)

maxphase (Optional, int) -- What phase are we interested in (% illum)

Returns:

List of [[date, alt, az, phase], ...]

'''

'''

Examples:

When will the moon transit at altitude 45 degrees and phase at least half,

and be bright in my skylight during MST nighttime hours?

when_at_position(ephem.Moon(), observer, 45., 180., 5, 12, 5, 75, 25)

When will the full moon rise exactly due east?

when_rise_set_at_position(ephem.Moon(), observer, 90., "rise", 5,

100, 20)

'''

start_triple = observer.date.triple()

start_triple = (start_triple[0], start_triple[1], int(start_triple[2]))

results = []

observer.date = ephem.Date(start_triple + (starttime, 0, 0))

end_date = ephem.Date((start_triple[0]+1, start_triple[1], start_triple[2],

starttime, 0, 0))

while observer.date <= end_date:

# At the beginning of each loop, observer.date is set to

# the starttime on a new day.

# We need to loop through to the endtime on the same day,

# then set observer.date to the starttime on the next day

# and continue.

# First get the time zone offset. We'll assume the offset

# is the same all day. That might make us be off by an hour

# on start and end times for a few hours twice a year.

local_hour = ephem.localtime(observer.date).hour

gmt_hour = observer.date.tuple()[3]

tzoffset = (gmt_hour - local_hour) % 24

# print "tzoffset:", tzoffset, "=", gmt_hour, "-", local_hour

daytriple = observer.date.triple()

for hour in discont_range(starttime+tzoffset, endtime+tzoffset+1, 24):

observer.date = ephem.Date(daytriple + (hour, 0, 0))

body.compute(observer)

alt = body.alt * DEGREES

az = body.az * DEGREES

if alt >= targetalt - slop and alt <= targetalt + slop \

and az >= targetaz - slop and az <= targetaz + slop:

if (body.phase > minphase and body.phase < maxphase):

results.append([observer.date, alt, az, body.phase])

return results

def when_rise_set_at_position(body, observer, targetaz, rise_set, slop=5.,

minphase=0, maxphase=100):

'''When will body rise at the target azimuth

in the year following the date of the observer passed in

(which defaults to today's date.)

Args:

body (ephem.Body) -- The celestial body to be calculated.

observer(ephem.Observer) -- the observing position.

targetaz (float) -- Azimuth in decimal degrees.

rise_set(string) -- "rise" or "set"

slop (Optional) -- How much slop to allow in the alt/az

positions each way (float, decimal degrees).

minphase (Optional, int) -- What phase are we interested in (% illum)

maxphase (Optional, int) -- What phase are we interested in (% illum)

Returns:

List of [[date, alt, az, phase], ...]

'''

start_triple = observer.date.triple()

start_triple = (start_triple[0], start_triple[1], int(start_triple[2]))

observer.date = ephem.Date((start_triple[0], start_triple[1],

start_triple[2], 0, 0, 0))

end_date = ephem.Date((start_triple[0]+1, start_triple[1], start_triple[2],

0, 0, 0))

results = []

while observer.date <= end_date:

if rise_set == "rise":

observer.date = observer.next_rising(body)

else:

observer.date = observer.next_setting(body)

body.compute(observer)

az = body.az * DEGREES

if az >= targetaz - slop and az <= targetaz + slop:

if (body.phase > minphase and body.phase < maxphase):

results.append([observer.date, observer.horizon*DEGREES, az,

body.phase])

# Push the date forward a little bit to make sure we get the

# next rise/set, not this one again.

observer.date += ephem.hour

return results

def Usage():

progname = os.path.basename(sys.argv[0])

print """Predict when the moon will rise, set or be at a specific position.

Usage: %s [-f] alt az start_hour end_hour

%s [-f] [rise|set] azimuth

If -f is specified, only full or nearly-full moons will be considered;

-p sets minimum percent illuminated, +p sets maximum.

Hours are specified in local time.

""" % (progname, progname)

sys.exit(0)

if __name__ == '__main__':

if len(sys.argv) < 3 or sys.argv[1] == '-h' or sys.argv[1] == '--help':

Usage()

args = sys.argv[1:]

minphase = 0

maxphase = 100

phasestr = "moon"

MIN_FULL = 90

while True:

if args[0] == '-f':

minphase = MIN_FULL

phasestr = "full moon"

args = args[1:]

elif args[0] == '-p':

minphase = int(args[1])

phasestr = "moon with at least %d%% illuminated" % minphase

args = args[2:]

elif args[0] == '+p':

maxphase = int(args[1])

phasestr = "moon with at most %d%% illuminated" % maxphase

args = args[2:]

else:

break

if minphase > 0:

if minphase == MIN_FULL:

phasestr = "full moon"

else:

if maxphase < 100:

phasestr = "moon between %d%% and %d%% illuminated" % (minphase,

maxphase)

else:

phasestr = "moon with at least %d%% illuminated" % minphase

elif maxphase < 100:

phasestr = "moon with at most %d%% illuminated" % maxphase

if len(args) == 2:

if args[0] != "rise" and args[0] != "set":

Usage()

try:

az = float(args[1])

except ValueError:

print "Error: azimuth must be a number, not %s" % args[1]

print

Usage()

results = when_rise_set_at_position(ephem.Moon(), observer,

az, args[0], 5, minphase, maxphase)

print """The %s will %s at azimuth %.1f during the next year at these times:

""" % (phasestr, args[0], az)

elif len(args) == 4:

try:

alt = float(args[0])

az = float(args[1])

except ValueError:

print "Error: alt/az must be numbers"

print

Usage()

try:

start_hour = int(args[2])

end_hour = int(args[3])

except ValueError:

print "Error: start/end times must be integers"

print

Usage()

results = when_at_position(ephem.Moon(), observer,

alt, az, start_hour, end_hour,

5, minphase, maxphase)

print """The %s will be at %.1f, %.1f

between %dh and %dh during the next year at these times:

""" % (phasestr, alt, az, start_hour, end_hour)

else:

Usage()

print "%26s %5s %5s %s" % ("Date/time", "alt", "az", "phase")

for r in results:

if len(r) > 3 and r[3]:

phase = "Illuminated: %2d%%" % r[3]

else:

phase = 'x'

# pargs = tuple(r[0:3] + [phase])

pargs = (ephem.localtime(r[0]).ctime(), r[1], r[2], phase)

print "%26s %5.1f %5.1f %s" % pargs