From d29fb19eb47c3f3ccbeb335a52fa09957f6326ab Mon Sep 17 00:00:00 2001 From: Peter Hinch Date: Wed, 29 Nov 2023 12:13:28 +0000 Subject: [PATCH] astronomy/sun_moon.py: Beta version and docs. --- astronomy/README.md | 65 +++++++++++++++++++++++ astronomy/sun_moon.py | 121 ++++++++++++++++++++++++++++-------------- 2 files changed, 145 insertions(+), 41 deletions(-) diff --git a/astronomy/README.md b/astronomy/README.md index e69de29..d584515 100644 --- a/astronomy/README.md +++ b/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. diff --git a/astronomy/sun_moon.py b/astronomy/sun_moon.py index c14f11f..cdf0132 100644 --- a/astronomy/sun_moon.py +++ b/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) - z1 = xe - dx - z2 = xe + dx + if b < 0: + z2 = (-b + sqrt(dis)) / (2 * a) # z2 is larger root in magnitude + 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 -# value regardless of the platform's epoch +# Re platform comparisons get_mjd returns the same value regardless of +# 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 - self._sr = None # Sunrise - self._ss = None # Sunset - self._mr = None # Moonrise - self._ms = None # Moon set + # Times in integer secs from midnight on current day (in local time) + # [sunrise, sunset, moonrise, moonset] + self._times = [None] * 4 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 - self._mr, self._ms = self.rise_set(False) # Moon - t = time.gmtime(time.time() + day * 86400) + t = time.gmtime(time.time() + day * spd) self._phase = moonphase(*t[:4]) + return self # Allow r.set_day().sunrise() - def sunrise(self, to=0): - return self._format(self._sr, to) + # variants: 0 secs since 00:00:00 localtime. 1 secs since MicroPython epoch + # (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): - return self._format(self._ss, to) + def sunset(self, variant: int = 0): + return self._format(self._times[1], variant) - def moonrise(self, to=0): - return self._format(self._mr, to) + def moonrise(self, variant: int = 0): + return self._format(self._times[2], variant) - def moonset(self, to=0): - return self._format(self._ms, to) + def moonset(self, variant: int = 0): + return self._format(self._times[3], variant) def moonphase(self): return self._phase # ***** API end ***** - def _format(self, n, to): - if to == 0: # Default: secs since Midnight + def adjust(self, n, day, idx): + 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 - t = d / 36525.0 - lst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - t * t * t / 38710000 + # modified to use the pre-computed 't' value from sin_alt + d = t * 36525 + 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 - t = (mjd - 51544.5) / 36525.0 - dec, ra = func(t) + mjd1 = (self.mjd - 51544.5) + hour / 24.0 + t1 = mjd1 / 36525.0 + # 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 salt + return self.sglat * sin(radians(dec)) + self.cglat * cos(radians(dec)) * cos(radians(tau)) # 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))