From 03c924d6b4a22962f5a5e7f94d6084b5c0980980 Mon Sep 17 00:00:00 2001 From: Peter Hinch Date: Mon, 25 Dec 2023 10:04:33 +0000 Subject: [PATCH] astronomy: Add moonphase.py. --- astronomy/README.md | 171 ++++++++++++++++++++++++------- astronomy/moonphase.py | 228 +++++++++++++++++++++++++++++++++++++++++ astronomy/package.json | 3 +- astronomy/sun_moon.py | 22 +--- 4 files changed, 364 insertions(+), 60 deletions(-) create mode 100644 astronomy/moonphase.py diff --git a/astronomy/README.md b/astronomy/README.md index 32b8705..8eea3d9 100644 --- a/astronomy/README.md +++ b/astronomy/README.md @@ -10,29 +10,36 @@ 2.2 [Methods](./README.md#22-methods) 2.3 [Effect of local time](./README.md#23-effect-of-local-time) 2.4 [Continuously running applications](./README.md#24-continuously-running-applications) -3. [The moonphase function](./README.md#3-the-moonphase-function) -4. [Utility functions](./README.md#4-utility-functions) -5. [Demo script](./README.md#5-demo-script) -6. [Scheduling events](./README.md#6-scheduling-events) +3. [Utility functions](./README.md#3-utility-functions) +4. [Demo script](./README.md#4-demo-script) +5. [Scheduling events](./README.md#5-scheduling-events) +6. [The moonphase module](./README.md#6-the-moonphase-module) + 6.1 [Constructor](./README.md#61-constructor) + 6.2 [Methods](./README.md#62-methods) + 6.3 [Usage examples](./README.md#63-usage-examples) + 6.4 [DST](./README.md#64-dst) Daylight savings time. 7. [Performance and accuracy](./README.md#7-performance-and-accuracy) + 7.1 [RiSet class](./README.md#71-riset-class) + 7.2 [moonphase class](./README.md#72-moonphase-class) # 1. Overview -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 specified as the current day plus an offset in -days. +The `sun_moon` 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 specified as the current day plus an +offset in days. It can also compute Civil, Nautical or Astronomical twilight +times. -It can also compute Civil, Nautical or Astronomical twilight times. A -`moonphase` function is also provided enabling the moon phase to be determined -for any date. +The `moonphase` module enables the moon phase to be determined for any date, and +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 code is currently under development but I don't anticipate breaking changes -to the API. +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 @@ -43,6 +50,8 @@ lunar clocks such as this one - the "lunartick": ## 1.2 Licensing and acknowledgements +#### sun_moon.py + Some code was ported from C/C++ as presented in "Astronomy on the Personal Computer" by Montenbruck and Pfleger, with mathematical improvements contributed by Marcus Mendenhall and Raul Kompaß. I (Peter Hinch) performed the port and @@ -56,6 +65,11 @@ source. I have not spotted any restrictions on use in the book. I am not a lawyer; I have no idea of the legal status of code based on sourcecode in a published work. +#### moonphase.py + +This was derived from unrestricted public sources and is released under the MIT +licence. + ## 1.3 Installation Installation copies files from the `astronomy` directory to a directory @@ -146,11 +160,8 @@ instance. 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. -* `moonphase()->float` Return current phase: 0.0 <= result < 1.0. 0.0 is new -moon, 0.5 is full. See [section 3](./README.md#3-the-moonphase-function) for -observations about this. * `set_lto(t)` Set local time offset `LTO` in hours relative to UTC. Primarily -intended for daylight saving time. The value is checked to ensure +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). The return value of the rise and set method is determined by the `variant` arg. @@ -223,28 +234,14 @@ 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. -# 3. The moonphase function - -This is a simple function whose provenance is uncertain. It appears to produce -valid results but I plan to implement a better solution. - -Args: -* `year: int` 4-digit year -* `month: int` 1..12 -* `day: int` Day of month 1..31 -* `hour: int` 0..23 - -Return value: -A float in range 0.0 <= result < 1.0, 0 being new moon, 0.5 being full moon. - -# 4. Utility functions +# 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. -# 5. Demo script +# 4. Demo script This produces output for the fixed date 4th Dec 2023 at three geographical locations. It can therefore be run on platforms where the system time is wrong. @@ -290,7 +287,7 @@ Maximum error 0. Expect 0 on 64-bit platform, 30s on 32-bit ``` Code comments show times retrieved from `timeanddate.com`. -# 6. Scheduling events +# 5. Scheduling events A likely use case is to enable events to be timed relative to sunrise and set. In simple cases this can be done with `asyncio`. This will execute a payload at @@ -386,9 +383,108 @@ try: finally: _ = asyncio.new_event_loop() ``` +# 6. The moonphase module + +This contains a single class `MoonPhase`. The term "machine time" below refers +to the time reported by the MicroPython `time` module. The "local time offset" +(LTO) passed to the constructor specifies the difference between machine time +and UTC based on system longitude. "Daylight saving time" (DST) allows reported +times to be offset to compensate for DST. Internally phases are calculated in +UTC, but where times are output they are adjusted for LTO and DST. + +It is recommended that the machine clock is not adjusted for DST because large +changes can play havoc with program timing as described above. To accommodate +DST, a `dst` function can be provided to the constructor. The module uses this +to adjust reported times. + +A `MoonPhase` instance has a time `datum`, which defaults to the instantiation +time. Phases are calculated with respect to this datum. It may be changed using +`.set_day` to enable future and past phases to be determined or to enable long +running applications to track time. + +## 6.1 Constructor + +* `lto:float=0, dst = lambda x: x` Local time offset in hours to UTC (-ve is +West); the value is checked to ensure `-15 < lto < 15`. `dst` is an optional +user defined function for Daylight Saving Time (DST). See +[section 6.4](./README.md#64-dst) + +## 6.2 Methods + +* `quarter(q: int, text: bool = True)` Return the time of a given quarter. Five +quarters are calculated around the instance datum. By default the time +is last midnight machine time with an optional offset in days `doff` added. The +`quarter` arg specifies the quarter with 0 and 4 being new moons and quarter 2 +being full. The `text` arg determines how the value is returned: as text or as +`int` is secs from the machine epoch. Results are adjusted for DST if a `dst` +function is provided to the constructor. +* `phase() -> float)` Returns moon phase where 0.0 <= phase < 1.0 with 0.5 being +full moon. The phase is that pertaining to the datum. +* `nextphase(, text: bool = True)` This is a generator function. Each iteration +of the generator returns three values: the phase number, the lunation number and +the datetime of the phase. The `text` arg is as per `.quarter()`, defining the +format of the datetime. +* `set_day(doff: float = 0)` Set the `MoonPhase` datum time to machine time plus +an offset in days: this may include a fractional part if `.phase()` is required +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. + +## 6.3 Usage examples + +```python +from moonphase import MoonPhase +mp = MoonPhase() # datum is midnight last night +print(f"Full moon, current lunation {mp.quarter(2)}") +mp.set_day(0.5) # Adjust datum to noon today machine time +print(f"Phase at Noon {mp.phase()}") +mp.set_day(182) # Set datum ahead 6 months +print(f"Lunation 1st new moon: {mp.quarter(0)}, 2nd new moon: {mp.quarter(4)}") +mp.set_day(0) # Reset datum to today +n = mp.nextphase() # Instantiate generator +for _ in range(8): + print(next(n)) +``` + +## 6.4 DST + +Daylight saving time depends on country and geographic location, and there is no +built-in MicroPython support. The moonphase module supports DST via an optional +user supplied function. DST does not affect the calculation of quarters or phase +which is based on the machine clock. If the machine clock runs at a fixed offset +to UTC (which is recommended), a DST function can be used to enable reported +results to reflect local time. + +A DST function takes as input a time measured in seconds since the machine epoch +(as returned by `time.time()`) and returns that number adjusted for local time. +The following example is for UK time, which adds one hour at 2:00 on the last +Sunday in March, reverting to winter time at 2:00 on the last Sunday in October. + +```python +def uk_dst(secs_epoch: int): # Change in March (3) and Oct (10) + t = time.gmtime(secs_epoch) + month = t[1] + mday = t[2] + wday = t[6] + winter = secs_epoch + summer = secs_epoch + 3600 # +1hr + if month in (1, 2, 11, 12): # Simple cases: depend only on month + return winter + if not month in (3, 10): + return summer # +1 hr + # We are in March or October. Find the day in month of last Sunday. + ld = (wday + 31 - mday) % 7 # weekday of 31st. + lsun = 31 - (1 + ld) % 7 + thresh = time.mktime((t[0], month, lsun, 2, 0, 0, 6, 0)) + return summer if ((secs_epoch >= thresh) ^ (month == 10)) else winter +``` # 7. Performance and accuracy +## 7.1 RiSet class + A recalculation is triggered whenever the 24 hour local time window is changed, such as calling `.set_day()` where the stored date changes. Normally two days of data are calculated, except where the local time is UTC where only one day is @@ -401,7 +497,6 @@ checked values corresponded with data computed on a platform with 64 bit floating point unit. The loss of precision from using a 32 bit FPU was no more than 3s. -The lunar phase calculation is poor. It is adequate for displaying a phase icon -or adjusting a pointer, but not good enough for predicting lunar quarters. +## 7.2 MoonPhase class -I plan to improve this, but it may be via a separate module. +TODO diff --git a/astronomy/moonphase.py b/astronomy/moonphase.py new file mode 100644 index 0000000..1726f49 --- /dev/null +++ b/astronomy/moonphase.py @@ -0,0 +1,228 @@ +# moonphase.py Calculate lunar phases + +# Source Tech\ Notes/Astronomy/astro_references/moontool.c +# The information for this was drawn from public domain sources including C code +# written by John Walker and Ron Hitchens in 1987-88 and released with the "licence" +# Do what thou wilt shall be the whole of the law". + +# Uses Python arbitrary length integers to maintain accuracy on platforms with +# 32-bit floating point. +# Copyright (c) Peter Hinch 2023 Released under the MIT license. + + +# Exports calc_phases() + +from math import radians, sin, cos, floor +import time +import array + +SYNMONTH = 29.53058868 # Synodic month (new Moon to new Moon) + + +# MEANPHASE -- Calculates time of the mean new Moon for a given base date. +# This argument K to this function is the precomputed synodic month index, given by: +# K = (year - 1900) * 12.3685 +# where year is expressed as a year and fractional year. +# sdate is days from 1900 January 0.5. Returns days from 1900 January 0.5 + + +def meanphase(sdate: float, k: int) -> float: + # Time in Julian centuries from 1900 January 0.5 + t = sdate / 36525 + t2 = t * t # Square for frequent use + t3 = t2 * t # Cube for frequent use + nt1 = 0.75933 + SYNMONTH * k + 0.0001178 * t2 - 0.000000155 * t3 + return nt1 + 0.00033 * sin(radians(166.56 + 132.87 * t - 0.009173 * t2)) + + +# TRUEPHASE -- Given a K value used to determine the mean phase of the new moon, +# and a phase no. (0..3), return the true, corrected phase time +# as integer Julian seconds. + + +def truephase(k: int, phi: int) -> int: + k += (0, 0.25, 0.5, 0.75)[phi] # Add phase to new moon time + t = k / 1236.85 # Time in Julian centuries from 1900 January 0.5 + t2 = t * t # Square for frequent use + t3 = t2 * t # Cube for frequent use + # Sun's mean anomaly + m = 359.2242 + 29.10535608 * k - 0.0000333 * t2 - 0.00000347 * t3 + # Moon's mean anomaly + mprime = 306.0253 + 385.81691806 * k + 0.0107306 * t2 + 0.00001236 * t3 + # Moon's argument of latitude + f = 21.2964 + 390.67050646 * k - 0.0016528 * t2 - 0.00000239 * t3 + if phi in (0, 2): # Corrections for New and Full Moon + pt = (0.1734 - 0.000393 * t) * sin(radians(m)) + pt += 0.0021 * sin(radians(2 * m)) + pt -= 0.4068 * sin(radians(mprime)) + pt += 0.0161 * sin(radians(2 * mprime)) + pt -= 0.0004 * sin(radians(3 * mprime)) + pt += 0.0104 * sin(radians(2 * f)) + pt -= 0.0051 * sin(radians(m + mprime)) + pt -= 0.0074 * sin(radians(m - mprime)) + pt += 0.0004 * sin(radians(2 * f + m)) + pt -= 0.0004 * sin(radians(2 * f - m)) + pt -= 0.0006 * sin(radians(2 * f + mprime)) + pt += 0.0010 * sin(radians(2 * f - mprime)) + pt += 0.0005 * sin(radians(m + 2 * mprime)) + else: # First or last quarter + pt = (0.1721 - 0.0004 * t) * sin(radians(m)) + pt += 0.0021 * sin(radians(2 * m)) + pt -= 0.6280 * sin(radians(mprime)) + pt += 0.0089 * sin(radians(2 * mprime)) + pt -= 0.0004 * sin(radians(3 * mprime)) + pt += 0.0079 * sin(radians(2 * f)) + pt -= 0.0119 * sin(radians(m + mprime)) + pt -= 0.0047 * sin(radians(m - mprime)) + pt += 0.0003 * sin(radians(2 * f + m)) + pt -= 0.0004 * sin(radians(2 * f - m)) + pt -= 0.0006 * sin(radians(2 * f + mprime)) + pt += 0.0021 * sin(radians(2 * f - mprime)) + pt += 0.0003 * sin(radians(m + 2 * mprime)) + pt += 0.0004 * sin(radians(m - 2 * mprime)) + pt -= 0.0003 * sin(radians(2 * m + mprime)) + if phi < 2: # First quarter correction + pt += 0.0028 - 0.0004 * cos(radians(m)) + 0.0003 * cos(radians(mprime)) + else: # Last quarter correction + pt += -0.0028 + 0.0004 * cos(radians(m)) - 0.0003 * cos(radians(mprime)) + pt = round(pt * 86400) # Integer seconds from here + pt += round(2_953_058_868 * 864 * k) // 1000_000 # round(SYNMONTH * k * 86400) + qq = 0.0001178 * t2 - 0.000000155 * t3 + qq += 0.00033 * sin(radians(166.56 + 132.87 * t - 0.009173 * t2)) + pt += round(qq * 86400) # qq amounts to 2s + return pt + 208_657_793_606 + + +def dt_to_text(tim): # Convert a time to text + t = time.localtime(tim) + return f"{t[2]:02}/{t[1]:02}/{t[0]:4} {t[3]:02}:{t[4]:02}:{t[5]:02}" + + +class MoonPhase: + verbose = True + + def __init__(self, lto: float = 0, dst=lambda x: x): + self.lto_s = self._check_lto(lto) # -15 < lto < 15 + # local time = UTC + lto .lto_s = offset in secs + self.dst = dst + # Datetimes in secs since hardware epoch based on UTC + # With epoch 1970 this could need long ints. + self.phases = array.array("q", (0,) * 5) + # Calculate Julian date of machine epoch + # Multiply by 100 to avoid fraction + jepoch = 244058750 # Julian date of Unix epoch (1st Jan 1970) * 100 + if time.gmtime(0)[0] == 2000: # Machine epoch + jepoch += 1095700 + jepoch *= 864 # Seconds from epoch + self.jepoch = jepoch + self.secs = 0 # Time of calling .set_day in secs UTC + self.set_day() # Populate array and .secs + if MoonPhase.verbose: + print(f"Machine time: {dt_to_text(time.time())}") + MoonPhase.verbose = False + + # Take offset in days from today, return time of last midnight in secs from machine epoch + # Take time of last midnight machine time in secs since machine epoch. Add a + # passed offset in days. Convert to UTC using LTO. The returned value is as + # if the hardware clock were running UTC. + def _midnight(self, doff: float = 0): # Midnight last night + days offset (UTC) + tl = round((time.time() // 86400 + doff) * 86400) # Target in local time + return tl - self.lto_s + + def set_lto(self, t: float): # Update the offset from UTC + self.lto_s = self._check_lto(t) # Localtime offset in secs + + def set_day(self, doff: float = 0): + self.secs = round(time.time() + doff * 86400 - self.lto_s) + start = self._midnight(doff) # Phases are calculated around this time (UTC) + self._populate(start) # Immediate return if .phases already OK + + def datum(self, text: bool = True): + t = self.secs + self.lto_s + return dt_to_text(t) if text else t + + def quarter(self, q: int, text: bool = True): + if not 0 <= q <= 4: + raise ValueError("Quarter nos must be from 0 to 4.") + tutc = self.phases[q] # Time of phase in secs UTC + # Adjust time: t is machine time in secs since machine epoch + t = self.dst(tutc + self.lto_s) # UTC secs from hardware epoch -> local time + return dt_to_text(t) if text else t # Secs since machine epoch + + # Return moon phase as 0.0 <= n < 1.0 by defaut for current datetime. + def phase(self) -> float: # doff: days offset with optional fraction + t = self.secs # As set by .set_day() + if not (self.phases[0] <= t <= self.phases[4]): # set_day was not called + self.set_day() # Assume today + prev = self.phases[0] + for n, phi in enumerate(self.phases): + if phi > t: + break # phi is upcoming phase time + prev = phi # Last phase before now + if prev == phi: # Day is day of new moon: use synodic month/4 + r = (t - prev) * 0.25 / 637860.715488 + if r < 0: + r = 1 - r + else: + r = (n - 1) * 0.25 + (t - prev) * 0.25 / (phi - prev) + return min(r, 0.999999) # Rare pathological results where r slightly > 1.0 + + def _next_lunation(self): # Use approx time of next full moon to advance + self._populate(round(self.phases[2] + SYNMONTH * 86400)) + + # toff: days offset with optional fraction + def nextphase(self, text: bool = True): + n = 0 + lun = 0 # Skip historic quarters + while True: + yield n, lun, self.quarter(n, text) + n += 1 + n %= 4 + if n == 0: + self._next_lunation() + lun += 1 + + def _check_lto(self, lto: float) -> int: + if not -15 < lto < 15: + raise ValueError("Invalid local time offset.") + return round(lto * 3600) + + # Populate the phase array. Fast return if phases are alrady correct. + # Find time of phases of the moon which surround the passed datetime. + # Five phases are found, starting and ending with the new moons which bound + # the specified lunation. + # Passed time, and the result in .phases, are seconds since hardware epoch + # adjusted for UTC: i.e. as if the RTC were running UTC rather than local time. + def _populate(self, t: int): + if self.phases[0] < t < self.phases[4]: + return # Nothing to do + # Return days since Jan 0.5 1900 as a float. Returns same value on 32 and 64 bits + def jd1900(t: int) -> float: + y, m, mday = time.localtime(t)[:3] + if m <= 2: + m += 12 + y -= 1 + b = round(y / 400 - y / 100 + y / 4) + mjm = 365 * y - 679004 + b + int(30.6001 * (m + 1)) + mday + return mjm - 15019.5 + + sdate: float = jd1900(t) # Days since 1900 January 0.5 + adate: float = sdate - 45 + yy, mm, dd = time.localtime(t)[:3] + k1: int = floor((yy + ((mm - 1) * (1.0 / 12.0)) - 1900) * 12.3685) # 365.25/SYNMONTH + adate = meanphase(adate, k1) # Find new moon well before current date + nt1: float = adate + while True: + adate += SYNMONTH # For each lunar month + k2: int = k1 + 1 + nt2: float = meanphase(adate, k2) + if nt1 <= sdate and nt2 > sdate: + break + nt1 = nt2 + k1 = k2 + # k is integer days since start of 1900, being the lunation number + # 1533, 1534 on both platforms. + for n, k in enumerate((k1, k1, k1, k1, k2)): + phi: int = truephase(k, n % 4) # Args lunation no., phase no. 0..3 + self.phases[n] = phi - self.jepoch # Julian datetime to secs since hardware epoch + # Datetimes in secs since hardware epoch based on UTC diff --git a/astronomy/package.json b/astronomy/package.json index df7ec76..889de37 100644 --- a/astronomy/package.json +++ b/astronomy/package.json @@ -1,7 +1,8 @@ { "urls": [ ["sched/sun_moon.py", "github:peterhinch/micropython-samples/astronomy/sun_moon.py"], - ["sched/sun_moon_test.py", "github:peterhinch/micropython-samples/astronomy/sun_moon_test.py"] + ["sched/sun_moon_test.py", "github:peterhinch/micropython-samples/astronomy/sun_moon_test.py"], + ["sched/moonphase.py", "github:peterhinch/micropython-samples/astronomy/moonphase.py"] ], "version": "0.1" } diff --git a/astronomy/sun_moon.py b/astronomy/sun_moon.py index f4e5491..a4f8eca 100644 --- a/astronomy/sun_moon.py +++ b/astronomy/sun_moon.py @@ -12,6 +12,7 @@ # Raul Kompaß perfomed major simplification of the maths for deriving rise and # set_times with improvements in precision with 32-bit floats. +# Moon phase now in separate module import time from math import sin, cos, sqrt, fabs, atan, radians, floor, pi @@ -103,23 +104,6 @@ def to_int(x): return None if x is None else round(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 - if itm >= 12: - itm -= 12 - term1 = floor(365.25 * (fty + 4712)) - term2 = floor(30.6 * itm + 0.5) - term3 = floor(floor((fty / 100) + 49) * 0.75) - 38 - tmp = term1 + term2 + day + 59 + hour / 24.0 - if tmp > 2299160.0: - tmp = tmp - term3 - phi = (tmp - 2451550.1) / 29.530588853 # 29.530588853 is length of lunar cycle (days) - return phi % 1 - - def minisun(t): # Output sin(dec), cos(dec), ra # returns the ra and dec of the Sun @@ -231,7 +215,6 @@ class RiSet: # 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) - self._phase = moonphase(*t[:4]) self.update(mjd) # Recalculate rise and set times return self # Allow r.set_day().sunrise() @@ -255,9 +238,6 @@ class RiSet: def tend(self, variant: int = 0): return self._format(self._times[5], variant) - def moonphase(self) -> float: - return self._phase - def set_lto(self, t): # Update the offset from UTC self.check_lto(t) # No need to recalc beause date is unchanged self.lto = round(t * 3600) # Localtime offset in secs