sun_moon_test.py: Add precomputed expected results.

astronomy/sun_moon.py

@@ -9,7 +9,7 @@
 # https://github.com/orgs/micropython/discussions/13075
 
 import time
-from math import sin, cos, sqrt, fabs, atan, radians, floor
+from math import sin, cos, sqrt, fabs, atan, radians, floor, pi
 
 LAT = 53.29756504536339  # Local defaults
 LONG = -2.102811634540558
@@ -64,7 +64,8 @@ def quad(ym, yz, yp):
             nz += 1
         if z1 < -1.0:
             z1 = z2
-    return nz, z1, z2, ye
+        return nz, z1, z2, ye
+    return 0, 0, 0, 0  # No roots
 
 
 # **** GET MODIFIED JULIAN DATE FOR DAY RELATIVE TO TODAY ****
@@ -116,20 +117,19 @@ def minisun(t):
     # in decimal hours, degs referred to the equinox of date and using
     # obliquity of the ecliptic at J2000.0 (small error for +- 100 yrs)
     # takes t centuries since J2000.0. Claimed good to 1 arcmin
-    p2 = 6.283185307
-    coseps = 0.91748
-    sineps = 0.39778
+    coseps = 0.9174805004
+    sineps = 0.397780757938
 
-    m = p2 * frac(0.993133 + 99.997361 * t)
+    m = 2 * pi * frac(0.993133 + 99.997361 * t)
     dl = 6893.0 * sin(m) + 72.0 * sin(2 * m)
-    l = p2 * frac(0.7859453 + m / p2 + (6191.2 * t + dl) / 1296000)
+    l = 2 * pi * frac(0.7859453 + m / (2 * pi) + (6191.2 * t + dl) / 1296000)
     sl = sin(l)
     x = cos(l)
     y = coseps * sl
     z = sineps * sl
     rho = sqrt(1 - z * z)
-    # dec = (360.0 / p2) * atan(z / rho)
-    ra = (48.0 / p2) * atan(y / (x + rho)) % 24
+    # dec = (360.0 / 2 * pi) * atan(z / rho)
+    ra = ((48.0 / (2 * pi)) * atan(y / (x + rho))) % 24
     return z, rho, ra
 
 
@@ -137,17 +137,15 @@ def minimoon(t):
     # takes t and returns the geocentric ra and dec
     # claimed good to 5' (angle) in ra and 1' in dec
     # tallies with another approximate method and with ICE for a couple of dates
-
-    p2 = 6.283185307
     arc = 206264.8062
-    coseps = 0.91748
-    sineps = 0.39778
+    coseps = 0.9174805004
+    sineps = 0.397780757938
 
     l0 = frac(0.606433 + 1336.855225 * t)  # mean longitude of moon
-    l = p2 * frac(0.374897 + 1325.552410 * t)  # mean anomaly of Moon
-    ls = p2 * frac(0.993133 + 99.997361 * t)  # mean anomaly of Sun
-    d = p2 * frac(0.827361 + 1236.853086 * t)  # difference in longitude of moon and sun
-    f = p2 * frac(0.259086 + 1342.227825 * t)  # mean argument of latitude
+    l = 2 * pi * frac(0.374897 + 1325.552410 * t)  # mean anomaly of Moon
+    ls = 2 * pi * frac(0.993133 + 99.997361 * t)  # mean anomaly of Sun
+    d = 2 * pi * frac(0.827361 + 1236.853086 * t)  # difference in longitude of moon and sun
+    f = 2 * pi * frac(0.259086 + 1342.227825 * t)  # mean argument of latitude
 
     # corrections to mean longitude in arcsec
     dl = 22640 * sin(l)
@@ -177,19 +175,19 @@ def minimoon(t):
     n += +21 * sin(-l + f)
 
     # ecliptic long and lat of Moon in rads
-    l_moon = p2 * frac(l0 + dl / 1296000)
+    l_moon = 2 * pi * frac(l0 + dl / 1296000)
     b_moon = (18520.0 * sin(s) + n) / arc
 
     # equatorial coord conversion - note fixed obliquity
     cb = cos(b_moon)
     x = cb * cos(l_moon)
     v = cb * sin(l_moon)
-    W = sin(b_moon)
-    y = coseps * v - sineps * W
-    z = sineps * v + coseps * W
+    w = sin(b_moon)
+    y = coseps * v - sineps * w
+    z = sineps * v + coseps * w
     rho = sqrt(1.0 - z * z)
-    # dec = (360.0 / p2) * atan(z / rho)
-    ra = (48.0 / p2) * atan(y / (x + rho)) % 24
+    # dec = (360.0 / 2 * pi) * atan(z / rho)
+    ra = ((48.0 / (2 * pi)) * atan(y / (x + rho))) % 24
     return z, rho, ra
 
 
@@ -207,8 +205,6 @@ class RiSet:
 
     # ***** API start *****
     # Examine Julian dates either side of current one to cope with localtime.
-    # TODO: Allow localtime offset to be varied at runtime for DST.
-    # TODO: relative=True arg for set_day. Allows entering absolute dates e.g. for testing.
     def set_day(self, day: int = 0):
         mjd = get_mjd(day)
         if self.mjd is None or self.mjd != mjd:

astronomy/sun_moon_test.py

@@ -1,12 +1,25 @@
-# sun_moon_test.py
+# sun_moon_test.py Test script for sun_moon.py
 
+# Copyright (c) Peter Hinch 2023
+# Released under the MIT license (see LICENSE)
 
-from .sun_moon import RiSet, abs_to_rel_days
+# On mip-installed host:
+# import sched.sun_moon_test
+# On PC in astronomy directory:
+# import sun_moon_test
+
+try:
+    from .sun_moon import RiSet, abs_to_rel_days
+except ImportError:  # Running on PC in astronomy directory
+    from sun_moon import RiSet, abs_to_rel_days
+
+nresults = []  # Times in seconds from local midnight
 
 
 def show(rs):
     print(f"Sun rise {rs.sunrise(3)} set {rs.sunset(3)}")
     print(f"Moon rise {rs.moonrise(3)} set {rs.moonset(3)}")
+    nresults.extend([rs.sunrise(), rs.sunset(), rs.moonrise(), rs.moonset()])
 
 
 print("4th Dec 2023: Seattle UTC-8")
@@ -28,6 +41,55 @@ for day in range(7):
     print(f"Day: {day}")
     show(rs)
 
+# Expected results as computed on Unix build (64-bit FPU)
+exp = [
+    27628,
+    58714,
+    85091,
+    46417,
+    20212,
+    71598,
+    2747,
+    41257,
+    29049,
+    57158,
+    82965,
+    46892,
+    29130,
+    57126,
+    None,
+    47460,
+    29209,
+    57097,
+    892,
+    47958,
+    29285,
+    57072,
+    5244,
+    48441,
+    29359,
+    57051,
+    9625,
+    48960,
+    29430,
+    57033,
+    14228,
+    49577,
+    29499,
+    57019,
+    19082,
+    50384,
+]
+print()
+max_error = 0
+for act, requirement in zip(nresults, exp):
+    if act is not None and requirement is not None:
+        err = abs(requirement - act)
+        max_error = max(max_error, err)
+        if err > 30:
+            print(f"Error {requirement - act}")
+print(f"Maximum error {max_error}. Expect 0 on 64-bit platform, 30s on 32-bit")
+
 # Times from timeanddate.com
 # Seattle
 # Sunrise 7:40 sunset 16:18 Moonrise 23:37 Moonset 12:53