diff --git a/astronomy/README.md b/astronomy/README.md index c001dbf..2922e3c 100644 --- a/astronomy/README.md +++ b/astronomy/README.md @@ -37,8 +37,6 @@ the dates and times of lunar quarters to be calculated. 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 `moonphase` module is currently under development: API changes are possible. - Moon phase options have been removed from `sun_moon` because accuracy was poor. ## 1.1 Applications @@ -73,9 +71,9 @@ licence. ## 1.3 Installation Installation copies files from the `astronomy` directory to a directory -`\lib\sched` on the target. This is for optional use with the +`\lib\sched` on the target. This directory eases optional use with the [schedule module](https://github.com/peterhinch/micropython-async/blob/master/v3/docs/SCHEDULE.md). -This may be done with the official +Installation may be done with the official [mpremote](https://docs.micropython.org/en/latest/reference/mpremote.html): ```bash $ mpremote mip install "github:peterhinch/micropython-samples/astronomy" @@ -87,7 +85,8 @@ On networked platforms it may alternatively be installed with ``` Currently these tools install to `/lib` on the built-in Flash memory. To install to a Pyboard's SD card [rshell](https://github.com/dhylands/rshell) may be used. -Move to `micropython-samples` on the PC, run `rshell` and issue: +Clone the repo and move to `micropython-samples` on the PC, run `rshell` and +issue: ``` > rsync astronomy /sd/sched ``` @@ -104,13 +103,14 @@ from sched.sun_moon import RiSet Time is a slippery concept when applied to global locations. This document uses the following conventions: * `UTC` The international time standard based on the Greenwich meridian. -* `LT (Local time)` Time as told on a clock at the device's location. May include +* `LT (Local time)` Time on a clock at the device's location. May include daylight saving time (`DST`). * `MT (Machine time)` Time defined by the platform's hardware clock. -* `LTO (Local time offset)` A `RiSet` instance contains a user supplied `LTO`. -The class computes rise and set times in UTC, using `LTO` to output results in -`LT` via `LT = UTC + LTO`. If an application maintains `LTO` to match `DST`, the -rise and set times will be in `LT`. +* `LTO (Local time offset)` A `RiSet` instance contains a user supplied `LTO` +intended for timezone support. The class computes rise and set times in UTC, +using `LTO` to compute results using `RESULT = UTC + LTO`. For output in `LT` +there are two options: periodically adjust `LTO` to handle DST or (better) +provide a `dst` function so that conversion is automatic. # 2. The RiSet class @@ -140,10 +140,16 @@ time (`MT`). (6 is Civil, 12 is Nautical, 18 is Astronomical). By default twilight times are not computed, saving some processor time. Offsets are positive numbers representing degrees below the horizon where twilight is deemed to start and end. +* `dst=lambda x: x` This is an optional user defined function for Daylight +Saving Time (DST). The assumption is that machine time is not changed, typically +permanently in winter time. A `dst` function handles seasonal changes. The +default assumes no DST is applicable. For how to write a DST function for a +given country see [section 6.4](./README.md#64-dst). By default when an application instantiates `RiSet` for the first time the constructor prints the system date and time. This can be inhibited by setting -the class variable `verbose` to `False`. +the class variable `verbose` to `False`. The purpose is to alert the user to a +common source of error where machine time is not set. ## 2.2 Methods @@ -162,8 +168,9 @@ horizon. * `has_risen(sun: bool)->bool` Returns `True` if the selected object has risen. * `has_set(sun: bool)->bool` Returns `True` if the selected object has set. * `set_lto(t)` Set local time offset `LTO` in hours relative to UTC. Primarily -intended for system longitude. The value is checked to ensure -`-15.0 < lto < 15.0`. See [section 2.3](./README.md#23-effect-of-local-time). +intended for timezone support, but this function can be used to support DST. The +value is checked to ensure `-15.0 < lto < 15.0`. See +[section 2.3](./README.md#23-effect-of-local-time). 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 @@ -173,7 +180,8 @@ with the moon at most locations, and with the sun in polar regions. Variants: * 0 Return integer seconds since midnight `LT` (or `None` if no event). * 1 Return integer seconds since since epoch of the MicroPython platform - (or `None`). This is machine time (`MT`) as per `time.time()`. + (or `None`). This allows comparisons with machine time (`MT`) as per + `time.time()`. * 2 Return text of form hh:mm:ss (or --:--:--) being local time (`LT`). Example constructor invocations: @@ -184,13 +192,16 @@ r = RiSet(lat=-33.87667, long=151.21, lto=11) # Sydney 33°52′04″S 151°12 ``` ## 2.3 Effect of local time -MicroPython has no concept of local time. The hardware platform has a clock +MicroPython has no concept of timezones. The hardware platform has a clock which reports machine time (`MT`): this might be set to local winter time or summer time. The `RiSet` instances' `LTO` should be set to represent the difference between `MT` and `UTC`. In continuously running applications it is best to avoid changing the hardware clock (`MT`) for reasons discussed below. -Daylight savings time should be implemented by changing the `RiSet` instances' -`LTO`. +Daylight savings time may be implemented in one of two ways: +* By changing the `RiSet` instances' `LTO` accordingly. +* Or by providing a `dst` function as discussed in +[section 6.4](./README.md#64-dst). This is the preferred solution as DST is then +handled automatically. Rise and set times are computed relative to UTC and then adjusted using the `RiSet` instance's `LTO` before being returned (see `.adjust()`). This means @@ -199,7 +210,7 @@ is used in determining rise and set times. The `.has_risen()`, `.has_set()` and `.is_up()` methods do use machine time (`MT`) and rely on `MT == UTC + LTO`: if `MT` has drifted, precision will be -reduced. +lost at times close to rise and set events. The constructor and the `set_day()` method set the instance's date relative to `MT`. They use only the date component of `MT`, hence they may be run at any @@ -231,16 +242,12 @@ synchronisation is required it is best done frequently to minimise the size of jumps. For this reason changing system time to accommodate daylight saving time is a -bad idea. It is usually best to run winter time all year round. Where a DST -change occurs, the `RiSet.set_lto()` method should be run to ensure that `RiSet` -operates in current local time. +bad idea. It is usually best to run winter time all year round and to use the +`dst` constructor arg to handle time changes. # 3. Utility functions `now_days() -> int` Returns the current time as days since the platform epoch. -`abs_to_rel_days(days: int) -> int` Takes a number of days since the Unix epoch -(1970,1,1) and returns a number of days relative to the current date. Platform -independent. This facilitates testing with pre-determined target dates. # 4. Demo script @@ -288,6 +295,10 @@ Maximum error 0. Expect 0 on 64-bit platform, 30s on 32-bit ``` Code comments show times retrieved from `timeanddate.com`. +The script includes some commented out code at the end. This tests `is_up`, +`has_risen` and `has_set` over 365 days. It is commented out to reduce printed +output. + # 5. Scheduling events A likely use case is to enable events to be timed relative to sunrise and set. @@ -436,7 +447,8 @@ to produce a time-precise value. The five quarters are calculated for the lunation including the midnight at the start of the specified day. * `set_lto(t:float)` Redefine the local time offset, `t` being in hours as per the constructor arg. -* `datum(text: bool = True)` Returns the current datum. +* `datum(text: bool = True)` Returns the current datum in secs since local epoch +or in human-readable text form. ## 6.3 Usage examples @@ -506,4 +518,6 @@ than 3s. ## 7.2 MoonPhase class -TODO +This uses Python's arbitrary precision integers to overcome the limitations of +32-bit floating point units. Results on 32 bit platforms match those on 64-bits +to within ~1 minute. Results match those on `timeanddate.com` within ~3 minutes. diff --git a/astronomy/sun_moon.py b/astronomy/sun_moon.py index de6fe87..dfdf421 100644 --- a/astronomy/sun_moon.py +++ b/astronomy/sun_moon.py @@ -33,22 +33,11 @@ LONG = -2.102811634540558 # right number of days for platform epoch at UTC. def now_days() -> int: secs_per_day = 86400 # 24 * 3600 - t = time.time() + t = RiSet.mtime() # Machine time as int. Can be overridden for test. t -= t % secs_per_day # Previous Midnight return round(t / secs_per_day) # Days since datum -# Convert number of days relative to the Unix epoch (1970,1,1) to a number of -# days relative to the current date. e.g. 19695 = 4th Dec 2023 -# Platform independent. -def abs_to_rel_days(days: int) -> int: - secs_per_day = 86400 # 24 * 3600 - now = now_days() # Days since platform epoch - if time.gmtime(0)[0] == 2000: # Machine epoch - now += 10957 - return days - now - - def quad(ym, yz, yp): # 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) @@ -186,16 +175,30 @@ def minimoon(t): class RiSet: verbose = True + # Riset.mtime() returns machine time as an int. The class variable tim is for + # test purposes only and allows the hardware clock to be overridden + tim = None - def __init__(self, lat=LAT, long=LONG, lto=0, tl=None): # Local defaults + @classmethod + def mtime(cls): + return round(time.time()) if cls.tim is None else cls.tim + + @classmethod + def set_time(cls, t): # Given time from Unix epoch set time + if time.gmtime(0)[0] == 2000: # Machine epoch + t -= 10957 * 86400 + cls.tim = t + + def __init__(self, lat=LAT, long=LONG, lto=0, tl=None, dst=lambda x: x): # Local defaults self.sglat = sin(radians(lat)) self.cglat = cos(radians(lat)) self.long = long self.check_lto(lto) # -15 < lto < 15 self.lto = round(lto * 3600) # Localtime offset in secs self.tlight = sin(radians(tl)) if tl is not None else tl + self.dst = dst self.mjd = None # Current integer MJD - # Times in integer secs from midnight on current day (in local time) + # Times in integer secs from midnight on current day (in machine time adjusted for DST) # [sunrise, sunset, moonrise, moonset, cvend, cvstart] self._times = [None] * 6 self.set_day() # Initialise to today's date @@ -212,9 +215,8 @@ class RiSet: if self.mjd is None or self.mjd != mjd: spd = 86400 # Secs per 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 - t = time.gmtime(time.time() + day * spd) + # time.time() assumes MicroPython clock is set to geographic local time + self._t0 = ((self.mtime() + day * spd) // spd) * spd self.update(mjd) # Recalculate rise and set times return self # Allow r.set_day().sunrise() @@ -243,30 +245,55 @@ class RiSet: self.lto = round(t * 3600) # Localtime offset in secs def has_risen(self, sun: bool): - now = round(time.time()) # Machine time - rt = self.sunrise(1) if sun else self.moonrise(1) # Machine time - if rt is None: - now += self.lto # UTC - t = (now % 86400) / 3600 # Time as UTC hour of day (float) - return self.sin_alt(t, sun) > 0 # Above horizon - return rt < now + return self.has_x(True, sun) def has_set(self, sun: bool): - now = round(time.time()) - st = self.sunset(1) if sun else self.moonset(1) - if st is None: - now += self.lto # UTC - t = (now % 86400) / 3600 # Time as UTC hour of day (float) - return self.sin_alt(t, sun) < 0 - return st < now + return self.has_x(False, sun) - def is_up(self, sun: bool): # Return current state of sun or moon - return self.has_risen(sun) and not self.has_set(sun) + # Return current state of sun or moon. The moon has a special case where it + # rises and sets in a 24 hour period. If its state is queried after both these + # events or before either has occurred, the current state depends on the order + # in which they occurred (the sun always sets afer it rises). + # The case is (.has_risen(False) and .has_set(False)) and if it occurs then + # .moonrise() and .moonset() must return valid times (not None). + def is_up(self, sun: bool): + hr = self.has_risen(sun) + hs = self.has_set(sun) + rt = self.sunrise() if sun else self.moonrise() + st = self.sunset() if sun else self.moonset() + if rt is None and st is None: # No event in 24hr period. + return self.above_horizon(sun) + # Handle special case: moon has already risen and set or moon has neither + # risen nor set, yet there is a rise and set event in the day + if not (hr ^ hs): + if not ((rt is None) or (st is None)): + return rt > st + if not (hr or hs): # No event has yet occurred + return rt is None + + return hr and not hs # Default case: up if it's risen but not set # ***** API end ***** + + # Generic has_risen/has_set function + def has_x(self, risen: bool, sun: bool): + if risen: + st = self.sunrise(1) if sun else self.moonrise(1) # DST- adjusted machine time + else: + st = self.sunset(1) if sun else self.moonset(1) + if st is not None: + return st < self.dst(self.mtime()) # Machine time + return False + + def above_horizon(self, sun: bool): + now = self.mtime() + self.lto # UTC + tutc = (now % 86400) / 3600 # Time as UTC hour of day (float) + return self.sin_alt(tutc, sun) > 0 # Object is above horizon + # Re-calculate rise and set times def update(self, mjd): - self._times = [None] * 6 + for x in range(len(self._times)): + self._times[x] = None # Assume failure days = (1, 2) if self.lto < 0 else (1,) if self.lto == 0 else (0, 1) tr = None # Assume no twilight calculations ts = None @@ -277,8 +304,8 @@ class RiSet: if self.tlight is not None: tr, ts = self.rise_set(True, True) mr, ms = self.rise_set(False, False) # Moon - # Adjust for local time. Store in ._times if value is in 24-hour - # local time window + # Adjust for local time and DST. Store in ._times if value is in + # 24-hour local time window self.adjust((sr, ss, mr, ms, tr, ts), day) self.mjd = mjd @@ -286,6 +313,7 @@ class RiSet: for idx, n in enumerate(times): if n is not None: n += self.lto + (day - 1) * 86400 + n = self.dst(n) # Adjust for DST on day of n h = n // 3600 if 0 <= h < 24: self._times[idx] = n @@ -332,7 +360,6 @@ class RiSet: tl = self.lstt(t, hour) + self.long # Local mean sidereal time adjusted for logitude return self.sglat * z + self.cglat * (x * cos(radians(tl)) + y * sin(radians(tl))) - # Modified to find sunrise and sunset only, not twilight events. # Calculate rise and set times of sun or moon for the current MJD. Times are # relative to that 24 hour period. def rise_set(self, sun, tl): diff --git a/astronomy/sun_moon_test.py b/astronomy/sun_moon_test.py index 25a5ea4..0b2b098 100644 --- a/astronomy/sun_moon_test.py +++ b/astronomy/sun_moon_test.py @@ -9,9 +9,19 @@ # import sun_moon_test try: - from .sun_moon import RiSet, abs_to_rel_days + from .sun_moon import RiSet except ImportError: # Running on PC in astronomy directory - from sun_moon import RiSet, abs_to_rel_days + from sun_moon import RiSet +import time + + +def mtime(h, m, t=None): + if t is None: + t = round(time.time()) + tm = (t // 86400) * 86400 + h * 3600 + m * 60 + print(time.gmtime(tm)) + return tm + nresults = [] # Times in seconds from local midnight @@ -20,27 +30,40 @@ 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() print("4th Dec 2023: Seattle UTC-8") rs = RiSet(lat=47.61, long=-122.35, lto=-8) # Seattle 47°36′35″N 122°19′59″W -rs.set_day(abs_to_rel_days(19695)) # 4th Dec 2023 +RiSet.set_time(19695 * 86400) +rs.set_day() show(rs) -print() print("4th Dec 2023: Sydney UTC+11") rs = RiSet(lat=-33.86, long=151.21, lto=11) # Sydney 33°52′04″S 151°12′36″E -rs.set_day(abs_to_rel_days(19695)) # 4th Dec 2023 +RiSet.set_time(19695 * 86400) +rs.set_day() show(rs) -print() print("From 4th Dec 2023: UK, UTC") rs = RiSet() for day in range(7): - rs.set_day(abs_to_rel_days(19695 + day)) # Start 4th Dec 2023 + RiSet.set_time(19695 * 86400) + rs.set_day(day) + # rs.set_day(abs_to_rel_days(19695 + day)) # Start 4th Dec 2023 print(f"Day: {day}") show(rs) + +print("4th Dec 2023: Sydney UTC+11 - test DST") +# Sydney 33°52′04″S 151°12′36″E +rs = RiSet(lat=-33.86, long=151.21, lto=11, dst=lambda x: x + 3600) +RiSet.set_time(19695 * 86400 + 86400 / 2) +rs.set_day() +# rs.set_day(abs_to_rel_days(19695)) # 4th Dec 2023 +show(rs) + + # Expected results as computed on Unix build (64-bit FPU) exp = [ 27628, @@ -79,6 +102,10 @@ exp = [ 57019, 19082, 50384, + 20212 + 3600, + 71598 + 3600, + 2747 + 3600, + 41257 + 3600, ] print() max_error = 0 @@ -88,6 +115,7 @@ for act, requirement in zip(nresults, exp): 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 @@ -97,3 +125,52 @@ print(f"Maximum error {max_error}. Expect 0 on 64-bit platform, 30s on 32-bit") # Sunrise 5:37 sunset 19:53 Moonrise 00:45 Moonset 11:28 # UK # Sunrise 8:04 sunset 15:52 Moonrise 23:02 Moonset 13:01 + + +def testup(t): # Time in secs since machine epoch + t = round((t // 86400) * 86400 + 60) # 1 minute past midnight + rs = RiSet() + RiSet.set_time(t) + rs.set_day() + tr = rs.moonrise() + ts = rs.moonset() + print(f"testup rise {rs.moonrise(2)} set {rs.moonset(2)}") + if tr is None and ts is None: + print(time.gmtime(t), "No moon events") + print( + f"Is up {rs.is_up(False)} Has risen {rs.has_risen(False)} has set {rs.has_set(False)}" + ) + return + # Initial state: not risen or set + assert not rs.has_set(False) + assert not rs.has_risen(False) + if tr is not None and (ts is None or ts > tr): + assert not rs.is_up(False) + rs.set_time(t + tr) + assert rs.has_risen(False) + assert rs.is_up(False) + assert not rs.has_set(False) + if ts is not None: + rs.set_time(t + ts) + assert rs.has_risen(False) + assert not rs.is_up(False) + assert rs.has_set(False) + return + if ts is not None: + assert rs.is_up(False) + rs.set_time(t + ts) + assert not rs.has_risen(False) + assert not rs.is_up(False) + assert rs.has_set(False) + if tr is not None: + rs.set_time(t + tr) + assert rs.has_risen(False) + assert rs.is_up(False) + assert rs.has_set(False) + return + print(f"Untested state tr {tr} ts {ts}") + + +# t = time.time() +# for d in range(365): +# testup(t + d * 86400)