astronomy/sun_moon.py: Beta version and docs.

astronomy/README.md

@@ -0,0 +1,65 @@
+# Astronomical calculations in MicroPython
+
+This module enables sun and moon rise and set times to be determined at any
+geographical location. Times are in seconds from midnight and refer to any
+event in a 24 hour period starting at midnight. The midnight datum is defined in
+local time. The start is a day being the current day plus an offset in days.
+
+A `moonphase` function is also provided enabling the moon phase to be determined
+for any date.
+
+The code was ported from C/C++ as presented in "Astronomy on the Personal
+Computer" by Montenbruck and Pfleger, with mathematical improvements contributed
+by Raul Kompaß and Marcus Mendenhall.
+
+Caveat. I am not an astronomer. If there are errors in the fundamental
+algorithms I am unlikely to be able to offer an opinion, still less a fix.
+
+The code is currently under development: the API may change.
+
+# The RiSet class
+
+## Constructor
+
+Args (float):
+* `lat=LAT` Latitude in degrees. Defaults are my location. :)
+* `long=LONG` Longitude in degrees (-ve is West).
+* `lto=0` Local time offset in hours (-ve is West).
+
+Methods:
+* `set_day(day: int = 0)` The arg is the offset from the current system date.
+Calling this with a changed arg causes the rise and set times to be updated.
+* `sunrise(variant: int = 0)` See below for details and the `variant` arg.
+* `sunset(variant: int = 0)`
+* `moonrise(variant: int = 0)`
+* `moonset(variant: int = 0)`
+* `moonphase()` Return current phase as a float: 0.0 <= result < 1.0. 0.0 is new
+moon, 0.5 is full.
+
+The return value of the rise and set method is determined by the `variant` arg.
+In all cases rise and set events are identified which occur in the current 24
+hour period. Note that a given event may be absent in the period: this can occur
+with the moon at most locations, and with the sun in polar regions.
+
+Variants:
+* 0 Return integer seconds since midnight local time (or `None` if no event).
+* 1 Return integer seconds since since epoch of the MicroPython platform
+ (or `None`).
+* 2 Return text of form hh:mm:ss (or --:--:--) being local time.
+
+# The moonphase function
+
+This is a simple function whose provenance is uncertain. I have a lunar clock
+which uses a C version of this which has run for 14 years without issue, but I
+can't vouch for its absolute accuracy over long time intervals. The Montenbruck
+and Pfleger version is very much more involved but they claim accuracy over
+centuries.
+
+Args (all integers):
+* `year` 4-digit year
+* `month` 1..12
+* `day` Day of month 1..31
+* `hour` 0..23
+
+Return value:  
+A float in range 0.0 <= result < 1.0, 0 being new moon, 0.5 being full moon.

astronomy/sun_moon.py

@@ -6,17 +6,18 @@
 
 # Source "Astronomy on the Personal Computer" by Montenbruck and Pfleger
 # ISBN 978-3-540-67221-0
+# Also contributions from Raul Kompaß and Marcus Mendenhall: see
+# https://github.com/orgs/micropython/discussions/13075
 
 import time
 from math import sin, cos, sqrt, fabs, atan, radians, floor
 
 LAT = 53.29756504536339  # Local defaults
 LONG = -2.102811634540558
-MOON_PHASE_LENGTH = 29.530588853
 
 
 def quad(ym, yz, yp):
-    # See Astronomy on the PC P48-49
+    # See Astronomy on the PC P48-49, plus contribution from Marcus Mendenhall
     # finds the parabola throuh the three points (-1,ym), (0,yz), (1, yp)
     # and returns the values of x where the parabola crosses zero
     # (roots of the quadratic)
@@ -29,9 +30,11 @@ def quad(ym, yz, yp):
     ye = (a * xe + b) * xe + c
     dis = b * b - 4.0 * a * c  # discriminant of y=a*x^2 +bx +c
     if dis > 0:  # parabola has roots
-        dx = 0.5 * sqrt(dis) / fabs(a)
+        if b < 0:
-        z1 = xe - dx
+            z2 = (-b + sqrt(dis)) / (2 * a)  # z2 is larger root in magnitude
-        z2 = xe + dx
+        else:
+            z2 = (-b - sqrt(dis)) / (2 * a)
+        z1 = (c / a) / z2  # z1 is always closer to zero
         if fabs(z1) <= 1.0:
             nz += 1
         if fabs(z2) <= 1.0:
@@ -53,8 +56,8 @@ def quad(ym, yz, yp):
 # (date(2000, 1, 1) - date(1970, 1, 1)).days * 24*60*60 = 946684800
 # (date(2000, 1, 1) - date(1970, 1, 1)).days = 10957
 
-# Re platform comparisons get_mjd does integer arithmetic and returns the same
+# Re platform comparisons get_mjd returns the same value regardless of
-# value regardless of the platform's epoch
+# the platform's epoch: integer days since 00:00 on 17 November 1858.
 def get_mjd(ndays: int = 0) -> int:  # Days offset from today
     secs_per_day = 86400  # 24 * 3600
     tsecs = time.time()  # Time now in secs since epoch
@@ -78,6 +81,7 @@ def to_int(x):
 
 
 # Approximate moon phase in range 0.0..1.0 0.0 is new moon, 0.5 full moon
+# Provenance of this cde is uncertain.
 def moonphase(year, month, day, hour):
     fty = year - floor((12.0 - month) / 10.0)
     itm = month + 9
@@ -89,7 +93,7 @@ def moonphase(year, month, day, hour):
     tmp = term1 + term2 + day + 59 + hour / 24.0
     if tmp > 2299160.0:
         tmp = tmp - term3
-    phi = (tmp - 2451550.1) / MOON_PHASE_LENGTH
+    phi = (tmp - 2451550.1) / 29.530588853  # 29.530588853 is length of lunar cycle (days)
     return phi % 1
 
 
@@ -181,53 +185,74 @@ def minimoon(t):
 
 
 class RiSet:
-    def __init__(self, lat=LAT, long=LONG):  # Local defaults
+    def __init__(self, lat=LAT, long=LONG, lto=0):  # Local defaults
         self.sglat = sin(radians(lat))
         self.cglat = cos(radians(lat))
         self.long = long
+        self.lto = round(lto * 3600)  # Localtime offset in secs
         self.mjd = None  # Current integer MJD
-        # Times in integer secs from midnight on current day
+        # Times in integer secs from midnight on current day (in local time)
-        self._sr = None  # Sunrise
+        # [sunrise, sunset, moonrise, moonset]
-        self._ss = None  # Sunset
+        self._times = [None] * 4
-        self._mr = None  # Moonrise
-        self._ms = None  # Moon set
         self.set_day()  # Initialise to today's date
 
     # ***** API start *****
     # 109μs on PBD-SF2W 166μs on ESP32-S3 394μs on RP2 (standard clocks)
-    def set_day(self, day=0):
+    def set_day(self, day: int = 0):
         mjd = get_mjd(day)
         if self.mjd is None or self.mjd != mjd:
             spd = 86400  # Secs per day
-            self._t0 = ((round(time.time()) + day * spd) // spd) * spd  # Midnight on target day
+            # ._t0 is time of midnight (local time) in secs since MicroPython epoch
+            # time.time() assumes MicroPython clock is set to local time
+            self._t0 = ((round(time.time()) + day * spd) // spd) * spd
+            self._times = [None] * 4
+            self._ms = None  # Moon set
+            for day in range(3):
+                self.mjd = mjd + day - 1
+                sr, ss = self.rise_set(True)  # Sun
+                mr, ms = self.rise_set(False)  # Moon
+                # Adjust for local time. Store in ._times if value is in 24-hour
+                # local time window
+                self.adjust(sr, day, 0)
+                self.adjust(ss, day, 1)
+                self.adjust(mr, day, 2)
+                self.adjust(ms, day, 3)
             self.mjd = mjd
-            self._sr, self._ss = self.rise_set(True)  # Sun
+            t = time.gmtime(time.time() + day * spd)
-            self._mr, self._ms = self.rise_set(False)  # Moon
-            t = time.gmtime(time.time() + day * 86400)
             self._phase = moonphase(*t[:4])
+        return self  # Allow r.set_day().sunrise()
 
-    def sunrise(self, to=0):
+    # variants: 0 secs since 00:00:00 localtime. 1 secs since MicroPython epoch
-        return self._format(self._sr, to)
+    # (relies on system being set to localtime). 2 human-readable text.
+    def sunrise(self, variant: int = 0):
+        return self._format(self._times[0], variant)
 
-    def sunset(self, to=0):
+    def sunset(self, variant: int = 0):
-        return self._format(self._ss, to)
+        return self._format(self._times[1], variant)
 
-    def moonrise(self, to=0):
+    def moonrise(self, variant: int = 0):
-        return self._format(self._mr, to)
+        return self._format(self._times[2], variant)
 
-    def moonset(self, to=0):
+    def moonset(self, variant: int = 0):
-        return self._format(self._ms, to)
+        return self._format(self._times[3], variant)
 
     def moonphase(self):
         return self._phase
 
     # ***** API end *****
-    def _format(self, n, to):
+    def adjust(self, n, day, idx):
-        if to == 0:  # Default: secs since Midnight
+        if n is not None:
+            n += self.lto + (day - 1) * 86400
+            h = n // 3600
+            if 0 <= h < 24:
+                self._times[idx] = n
+
+    def _format(self, n, variant):
+        if variant == 0:  # Default: secs since Midnight (local time)
             return n
-        elif to == 1:  # Secs since epoch
+        elif variant == 1:  # Secs since epoch of MicroPython platform
             return None if n is None else n + self._t0
-        # to == 3
+        # variant == 3
         if n is None:
             return "--:--:--"
         else:
@@ -235,26 +260,31 @@ class RiSet:
             mi, sec = divmod(tmp, 60)
             return f"{hr:02d}:{mi:02d}:{sec:02d}"
 
-    def lmst(self, mjd):
+    # See https://github.com/orgs/micropython/discussions/13075
+    def lmstt(self, t):
         # Takes the mjd and the longitude (west negative) and then returns
         # the local sidereal time in hours. Im using Meeus formula 11.4
         # instead of messing about with UTo and so on
-        d = mjd - 51544.5
+        # modified to use the pre-computed 't' value from sin_alt
-        t = d / 36525.0
+        d = t * 36525
-        lst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - t * t * t / 38710000
+        df = frac(d)
+        c1 = 360
+        c2 = 0.98564736629
+        dsum = c1 * df + c2 * d  # no large integer * 360 here
+        lst = 280.46061837 + dsum + 0.000387933 * t * t - t * t * t / 38710000
         return (lst % 360) / 15.0 + self.long / 15
 
     def sin_alt(self, hour, func):
         # Returns the sine of the altitude of the object (moon or sun)
         # at an hour relative to the current date (mjd)
-        mjd = self.mjd + hour / 24.0
+        mjd1 = (self.mjd - 51544.5) + hour / 24.0
-        t = (mjd - 51544.5) / 36525.0
+        t1 = mjd1 / 36525.0
-        dec, ra = func(t)
+        # print(f"sin_alt mjd0={mjd0:.7f} t0={t0:.9f} mjd1={mjd1:.7f} t1={t1:.9f}")
+        dec, ra = func(t1)
         # hour angle of object: one hour = 15 degrees. Note lmst() uses longitude
-        tau = 15.0 * (self.lmst(mjd) - ra)
+        tau = 15.0 * (self.lmstt(t1) - ra)
         # sin(alt) of object using the conversion formulas
-        salt = self.sglat * sin(radians(dec)) + self.cglat * cos(radians(dec)) * cos(radians(tau))
+        return self.sglat * sin(radians(dec)) + self.cglat * cos(radians(dec)) * cos(radians(tau))
-        return salt
 
     # Modified to find sunrise and sunset only, not twilight events.
     def rise_set(self, sun):
@@ -296,8 +326,17 @@ class RiSet:
 
 
 r = RiSet()
+# Seattle RiSet(lat=47.61, long=-122.35, lto=-8)
+# t = time.ticks_us()
+# r.set_day()
+# print("Elapsed us", time.ticks_diff(time.ticks_us(), t))
 for d in range(7):
     print(f"Day {d}")
     r.set_day(d)
     print(f"Sun rise {r.sunrise(3)} set {r.sunset(3)}")
     print(f"Moon rise {r.moonrise(3)} set {r.moonset(3)}")
+
+print(r.set_day().sunrise(0))
+# for d in range(30):
+#    r.set_day(d)
+#    print(round(r.moonphase() * 1000))