rewrite with skyfield

pyphoon.py

@@ -94,7 +94,7 @@ def putmoon(pctphase, lines, atfiller, hemisphere):  # pylint: disable=too-many-
 
                 # rotate char upside-down if needed
                 char = char.translate(str.maketrans("().`_'",
-                                                    ")(`.^,"))
+                                                    "!!!!!!"))
 
             if char != '@':
                 putchar(char)

xaphoon.py

@@ -1,103 +1,43 @@
 #!/usr/bin/env python
-from datetime import datetime, timezone
+"""Xaphoon - Displays the phase of the moon as well as other related information."""
-import math
+
 import time
+
 from argparse import ArgumentParser
-import ephem
+from datetime import datetime, timezone
+from skyfield import almanac
+from skyfield_data import get_skyfield_data_path
+import skyfield.api
+
 from pyphoon import putmoon
 
-# Second resolution for culmination/illumination calculations
+# Initialize certain skyfield parameters globally
-DAY_INCREMENT=1/86400
+sf_load = skyfield.api.Loader(get_skyfield_data_path(), expire=False) # loader
+ts = sf_load.timescale(builtin=False) # timescale
+eph = sf_load('de421.bsp') # ephemerides
+earth, sun, moon = eph['Earth'], eph['Sun'], eph['Moon'] # moooon
 
-def to_deg(rad):
+def to_timestr(t, date=False, local=True):
-    """Convert radians to a displayable integer number of degrees."""
+    """Convert a skyfield time to a time string, optionally in the local time zone."""
-    return round(math.degrees(rad))
+    t = t.utc_datetime()
-
-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)
+        t = t.astimezone()
-    else:
+    if date:
-        date = date.datetime()
+        return t.strftime('%Y-%m-%d %H:%M:%S')
-    return date.strftime("%H:%M:%S")
+    return t.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():
+    """Main function
+
+    Parses arguments, calculates values, and displays them.
+    """
     parser = ArgumentParser()
     parser.add_argument("lat",
-                        help="Observer latitude")
+                        help="Observer latitude",
-    parser.add_argument("long",
+                        type=float)
-                        help="Observer longitude")
+    parser.add_argument("lon",
+                        help="Observer longitude",
+                        type=float)
     parser.add_argument("elevation",
                         help="Observer elevation in meters",
                         type=int)
@@ -111,47 +51,42 @@ def main():
                         type=int)
     args = parser.parse_args()
 
-    now = ephem.date(datetime.fromtimestamp(args.time, timezone.utc))
+    t = ts.from_datetime(datetime.fromtimestamp(args.time, timezone.utc)) # current time
-    print(f"Current time: {to_timestr(now)}")
 
-    me = ephem.Observer()
+    print(f"Current time: {to_timestr(t)}")
-    me.date = now
-    me.lat = args.lat
-    me.lon = args.long
-    me.elevation = args.elevation
 
-    moon = ephem.Moon(me)
+    obs_geo = skyfield.api.wgs84.latlon(args.lat, args.lon,
+                       elevation_m=args.elevation) # geographic position vector
+    obs = earth + obs_geo # barycentric position vector
 
-    az = to_deg(moon.az)
+    moon_apparent = obs.at(t).observe(moon).apparent()
-    el = to_deg(moon.alt)
+    el, az, _ = moon_apparent.altaz('standard')
-    print(f"Az: {az}° El: {el}°")
+    print(f"Az: {az.degrees:.0f}° El: {el.degrees:.0f}°")
 
-    rising = find_target_rising(moon, me)
+    # Find relevant moonrise. el is based on apparent location, so accounts
-    setting = me.next_setting(moon)
+    # for atmospheric refraction. y shouldn't be needed unless user is near
-
+    # one of the poles, so ignored for now. First [0] discards y (second
-    print (f"Rise: {to_timestr(rising)} Set: {to_timestr(setting)}")
+    # element of tuple); second [] selects from array of moonrises/moonsets
-
+    if el.degrees > 0:
-    culm = find_culmination(moon, me, rising, setting)
+        # Moon is up. Find last moonrise in the past 24 hours.
-
+        moonrise = almanac.find_risings(obs, moon, t-1, t)[0][-1]
-    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
+        # Moon is not up. Find first moonrise in the next 24 hours.
+        moonrise = almanac.find_risings(obs, moon, t, t+1)[0][0]
 
-    print(putmoon(full_cycle_phase, 20, '@', 'northern' if me.lat > 0 else 'southern'))
+    # Find first moonset in the next 24 hours after moonrise.
+    moonset = almanac.find_settings(obs, moon, moonrise, moonrise+1)[0][0]
+    print(f"Rise: {to_timestr(moonrise)} Set: {to_timestr(moonset)}")
+
+    transit = almanac.find_transits(obs, moon, moonrise, moonrise+1)[0]
+    print(f"Transit: {to_timestr(transit)}")
+
+    phase = almanac.moon_phase(eph, t)
+    print(f"Phase: {phase.degrees:.0f}°")
+
+    illum = moon_apparent.fraction_illuminated(sun)
+    print(f"Illumination: {illum*100:.0f}%")
+
+    print(putmoon(phase.degrees/360, 21, '@', 'north' if args.lat > 0 else 'south'))
 
 main()