#!/usr/bin/env python from datetime import datetime, timezone import math import time from argparse import ArgumentParser import ephem from pyphoon import putmoon # Second resolution for culmination/illumination calculations DAY_INCREMENT=1/86400 def to_deg(rad): """Convert radians to a displayable integer number of degrees.""" return round(math.degrees(rad)) def to_timestr(date, local=True): """Convert a pyephem date to a time string in the local time zone.""" if local: date = ephem.localtime(date) else: date = date.datetime() return date.strftime("%H:%M:%S") def find_target_rising(moon, me): """Return the relevant moonrise to base display and calculations off of.""" if moon.alt == 0: # i would love a better way to do this me = me.copy() me.date = me.previous_rising(moon) return me.next_rising(moon) if moon.alt > 0: return me.previous_rising(moon) # moon.alt < 0 return me.next_rising(moon) def cmp_culmination(moon, me, t): """Determine whether the culmination is before, after, or at t. Returns 0 if t is the culmination, -1 if t if culmination is before t, or 1 if culmination is after t. Assumes there is exactly one peak elevation, which seems to cause error of up to about 7 seconds due to float precision. """ me.date = t - DAY_INCREMENT moon.compute(me) e1 = moon.alt me.date = t moon.compute(me) e2 = moon.alt me.date = t + DAY_INCREMENT moon.compute(me) e3 = moon.alt if e1 > e2: return -1 if e3 > e2: return 1 return 0 def find_culmination(moon, me, rising, setting): """Finds culmination via binary search. Assumes rising and setting are from same pass. """ moon = moon.copy() me = me.copy() t1 = rising t3 = setting while True: t2 = (t1 + t3) / 2 match cmp_culmination(moon,me,t2): case 0: return ephem.date(t2) case -1: t3 = t2 case 1: t1 = t2 def cmp_illumination(moon, me, t): """Determine whether the moon is waxing, waning, or either full or new. Returns 0 if the moon is either full or new, -1 if moon is waning, or 1 if moon is waxing. """ moon = moon.copy() me = me.copy() me.date = t - DAY_INCREMENT moon.compute(me) i1 = moon.moon_phase me.date = t + DAY_INCREMENT moon.compute(me) i2 = moon.moon_phase if i1 > i2: return -1 if i1 < i2: return 1 return 0 def main(): parser = ArgumentParser() parser.add_argument("lat", help="Observer latitude") parser.add_argument("long", help="Observer longitude") parser.add_argument("elevation", help="Observer elevation in meters", type=int) parser.add_argument("-c", "--columns", help="Number of columns for the output to use (default 70)", default=70, type=int) parser.add_argument("-t", "--time", help="Unix epoch time to perform calculations for", default=time.time(), type=int) args = parser.parse_args() now = ephem.date(datetime.fromtimestamp(args.time, timezone.utc)) print(f"Current time: {to_timestr(now)}") me = ephem.Observer() me.date = now me.lat = args.lat me.lon = args.long me.elevation = args.elevation moon = ephem.Moon(me) az = to_deg(moon.az) el = to_deg(moon.alt) print(f"Az: {az}° El: {el}°") rising = find_target_rising(moon, me) setting = me.next_setting(moon) print (f"Rise: {to_timestr(rising)} Set: {to_timestr(setting)}") culm = find_culmination(moon, me, rising, setting) print(f"Culmination: {to_timestr(culm)}") direction = cmp_illumination(moon, me, now) match direction: case -1: direction_indicator = '-' case 0: direction_indicator = '' case 1: direction_indicator = '+' print(f"Phase: {moon.moon_phase:.0%}{direction_indicator}") # Convert illumination percentage and waxing/waning status to percent through full cycle if direction < 0: # waning full_cycle_phase = 1 - (moon.moon_phase / 2) else: full_cycle_phase = moon.moon_phase / 2 print(putmoon(full_cycle_phase, 20, '@', 'northern' if me.lat > 0 else 'southern')) main()