xaphoon/xaphoon.py

#!/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()