From d3790626c45dc8112ee758a04b9202efd846937e Mon Sep 17 00:00:00 2001 From: Peter Hinch Date: Sat, 2 Dec 2023 17:01:06 +0000 Subject: [PATCH] sun_moon_test.py: Add precomputed expected results. --- astronomy/sun_moon.py | 46 ++++++++++++-------------- astronomy/sun_moon_test.py | 66 ++++++++++++++++++++++++++++++++++++-- 2 files changed, 85 insertions(+), 27 deletions(-) diff --git a/astronomy/sun_moon.py b/astronomy/sun_moon.py index c1befc1..56f98a0 100644 --- a/astronomy/sun_moon.py +++ b/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: diff --git a/astronomy/sun_moon_test.py b/astronomy/sun_moon_test.py index 6376014..25a5ea4 100644 --- a/astronomy/sun_moon_test.py +++ b/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