astronomy: Add moonphase.py.

master
Peter Hinch 2023-12-25 10:04:33 +00:00
rodzic ecd9b4d94e
commit 03c924d6b4
4 zmienionych plików z 364 dodań i 60 usunięć

Wyświetl plik

@ -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

Wyświetl plik

@ -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

Wyświetl plik

@ -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"
}

Wyświetl plik

@ -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