diff --git a/astronomy/README.md b/astronomy/README.md index 98dc38b..056b41c 100644 --- a/astronomy/README.md +++ b/astronomy/README.md @@ -20,11 +20,51 @@ 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. The sourcecode exists in the book and also on an accompanying CD-R. The file `CDR_license.txt` contains a copy of the -license file on the CD-R. I am not a lawyer; I have no idea of the legal status -of code translated from that in a published work. +license file on the disk, which contains source, executable code, and databases. +This module (obviously) only references the source. I am not a lawyer; I have no +idea of the legal status of code translated from sourcecode in a published work. + +## Installation + +Installation copies files from the `astronomy` directory to a directory +`\lib\sched` on the target. This is for 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 +[mpremote](https://docs.micropython.org/en/latest/reference/mpremote.html): +```bash +$ mpremote mip install "github:peterhinch/micropython-samples/astronomy" +``` +On networked platforms it may alternatively be installed with +[mip](https://docs.micropython.org/en/latest/reference/packages.html). +```py +>>> mip.install("github:peterhinch/micropython-samples/astronomy") +``` +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: +``` +> rsync astronomy /sd/sched +``` +`mip` installs the following files in the `sched` directory. +* `sun_moon.py` +* `sun_moon_test.py` A test/demo script. +After installation the `RiSet` class may be accessed with +```python +from sched.sun_moon import RiSet +``` # The RiSet class +This holds the local geographic coordinates and the localtime offset relative to +UTC. It is initialised to the current date and can provide the times of rise and +set events occurring within a 24 hour window starting at 00:00:00 local time. +The `RiSet` instance's date may be changed allowing rise and set times to be +retrieved for other 24 hour windows. + +Rise and set times may be retrieved in various formats including seconds from +local midnight: this may be used to enable the timing of actions relative to a +rise or set event. + ## Constructor Args (float): @@ -33,17 +73,17 @@ Args (float): * `lto=0` Local time offset in hours to UTC (-ve is West). Methods: -* `set_day(day: int = 0, relative=True)` `day` is the offset from the current -system date if `relative` is `True` otherwise it is the offset from the platform -epoch. If `day` is changed the rise and set times are updated. +* `set_day(day: int = 0)` `day` is the offset in days from the current system +date. If `day` is changed compared to the object's currently stored value its +rise and set times are updated. Returns the `RiSet` instance. * `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. -* `set_lto(t)` Update localtime offset to UTC (for daylight saving time). Rise -and set times are updated if the lto is changed. +* `set_lto(t)` Set localtime offset in hours relative to UTC. Primarily intended +for daylight saving time. Rise and set times are updated if the lto is changed. 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 @@ -66,7 +106,7 @@ r = RiSet(lat=-33.87667, long=151.21, lto=11) # Sydney 33°52′04″S 151°12 # 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 +which uses the original C code. This 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. @@ -79,3 +119,126 @@ Args: Return value: A float in range 0.0 <= result < 1.0, 0 being new moon, 0.5 being full moon. + +# 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. + +# 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. +To run issue: +```python +import sched.sun_moon_test +``` +Expected output: +```python +>>> import sched.sun_moon_test +4th Dec 2023: Seattle UTC-8 +Sun rise 07:40:09 set 16:18:15 +Moon rise 23:38:11 set 12:53:40 + +4th Dec 2023: Sydney UTC+11 +Sun rise 05:36:24 set 19:53:21 +Moon rise 00:45:55 set 11:27:14 + +From 4th Dec 2023: UK, UTC +Day: 0 +Sun rise 08:04:34 set 15:52:13 +Moon rise 23:03:15 set 13:01:04 +Day: 1 +Sun rise 08:05:54 set 15:51:42 +Moon rise --:--:-- set 13:10:35 +Day: 2 +Sun rise 08:07:13 set 15:51:13 +Moon rise 00:14:40 set 13:18:59 +Day: 3 +Sun rise 08:08:28 set 15:50:49 +Moon rise 01:27:12 set 13:27:08 +Day: 4 +Sun rise 08:09:42 set 15:50:28 +Moon rise 02:40:34 set 13:35:56 +Day: 5 +Sun rise 08:10:53 set 15:50:10 +Moon rise 03:56:44 set 13:46:27 +Day: 6 +Sun rise 08:12:01 set 15:49:56 +Moon rise 05:18:32 set 14:00:11 +>>> +``` +Code comments show times retrieved from `timeanddate.com`. + +# 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`. +```python +from sched.sun_moon import RiSet +import time +rs = RiSet() +tsecs = time.time() # Time now in secs since epoch +tsecs -= tsecs % 86400 # Last midnight in secs since epoch +tmidnight = tsecs +async def do_sunrise(): + while True: + toff = time.time() - tmidnight # Seconds relative to midnight + if toff > 0: # Midnight has passed, wait for sunrise + twait = rs.sunrise() - toff # Assumes a latitude where sun must rise + else: # Wait for tomorrow + twait = tmidnight + 86400 + toff + await asyncio.sleep(twait) + if toff > 0: + # Turn the lights off, or whatever +``` +An alternative, particularly suited to more complex cases, is to use the +[schedule module](https://github.com/peterhinch/micropython-async/blob/master/v3/docs/SCHEDULE.md). +This allows more intuitive coding without the epoch calculations. The following +is a minimal example: +```python +import uasyncio as asyncio +from sched.sched import schedule +from sched.sun_moon import RiSet + +async def turn_off_lights(rs): # Runs at 00:01:00 + rs.set_day() # Re-calculate for new daylight + asyncio.sleep(rs.sunrise() - 60) + # Actually turn them off + +async def main(): + rs = RiSet() # May need args for your location + await schedule(turn_off_lights, rs, hrs=0, mins=1) # Never terminates + +try: + asyncio.run(main()) +finally: + _ = asyncio.new_event_loop() +``` +This approach lends itself to additional triggers and events: +```python +import uasyncio as asyncio +from sched.sched import schedule, Sequence +from sched.sun_moon import RiSet + +async def turn_off_lights(t): + asyncio.sleep(t) + # Actually turn them off + +async def main(): + rs = RiSet() # May need args for your location + seq = Sequence() # A Sequence comprises one or more schedule instances + asyncio.create_task(schedule(seq, "off", hrs=0, mins=1)) + # Can schedule other events here + async for args in seq: + if args[0] == "off": # Triggered at 00:01 hrs (there might be other triggers) + rs.set_day() # Re-calculate for new day + asyncio.create_task(turn_off_lights(rs.sunrise() - 60)) + +try: + asyncio.run(main()) +finally: + _ = asyncio.new_event_loop() +``` diff --git a/astronomy/package.json b/astronomy/package.json new file mode 100644 index 0000000..df7ec76 --- /dev/null +++ b/astronomy/package.json @@ -0,0 +1,7 @@ +{ + "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"] + ], + "version": "0.1" +} diff --git a/astronomy/sun_moon.py b/astronomy/sun_moon.py index cc18274..c1befc1 100644 --- a/astronomy/sun_moon.py +++ b/astronomy/sun_moon.py @@ -14,6 +14,30 @@ from math import sin, cos, sqrt, fabs, atan, radians, floor LAT = 53.29756504536339 # Local defaults LONG = -2.102811634540558 +# MicroPython wanton epochs: +# time.gmtime(0)[0] = 1970 or 2000 depending on platform. +# On CPython: +# (date(2000, 1, 1) - date(1970, 1, 1)).days * 24*60*60 = 946684800 +# (date(2000, 1, 1) - date(1970, 1, 1)).days = 10957 + +# Return time now in days relative to platform epoch. +def now_days() -> int: + secs_per_day = 86400 # 24 * 3600 + t = time.time() + 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 @@ -45,26 +69,15 @@ def quad(ym, yz, yp): # **** GET MODIFIED JULIAN DATE FOR DAY RELATIVE TO TODAY **** -# Takes the system time in seconds from 1 Jan 70 & returns -# modified julian day number defined as mjd = jd - 2400000.5 +# Returns modified julian day number defined as mjd = jd - 2400000.5 # Deals only in integer MJD's: the JD of just after midnight will always end in 0.5 # hence the MJD of an integer day number will always be an integer -# MicroPython wanton epochs: -# time.gmtime(0)[0] = 1970 or 2000 depending on platform. -# (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 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 +def get_mjd(ndays: int = 0) -> int: secs_per_day = 86400 # 24 * 3600 - tsecs = time.time() # Time now in secs since epoch - tsecs -= tsecs % secs_per_day # Time last midnight - tsecs += secs_per_day * ndays # Specified day - days_from_epoch = round(tsecs / secs_per_day) # Days from 1 Jan 70 - # tsecs += secs_per_day # 2 # Noon - # mjepoch = -10957.5 # 40587 - 51544.5 # Modified Julian date of C epoch (1 Jan 70) + days_from_epoch = now_days() + ndays # Days since platform epoch mjepoch = 40587 # Modified Julian date of C epoch (1 Jan 70) if time.gmtime(0)[0] == 2000: mjepoch += 10957 @@ -196,7 +209,7 @@ class RiSet: # 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, relative: bool = True): + 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 @@ -226,7 +239,10 @@ class RiSet: return self._phase def set_lto(self, t): # Update the offset from UTC - pass # TODO + lto = round(t * 3600) # Localtime offset in secs + if self.lto != lto: # changed + self.lto = lto + self.update(self.mjd) # ***** API end ***** # Re-calculate rise and set times @@ -328,21 +344,3 @@ class RiSet: if t_rise is not None and t_set is not None: break # All done return to_int(t_rise), to_int(t_set) # Convert to int preserving None values - - -# r = RiSet() -# r = RiSet(lat=47.61, long=-122.35, lto=-8) # Seattle 47°36′35″N 122°19′59″W -r = RiSet(lat=-33.86, long=151.21, lto=11) # Sydney 33°52′04″S 151°12′36″E -# 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)) diff --git a/astronomy/sun_moon_test.py b/astronomy/sun_moon_test.py new file mode 100644 index 0000000..6376014 --- /dev/null +++ b/astronomy/sun_moon_test.py @@ -0,0 +1,37 @@ +# sun_moon_test.py + + +from .sun_moon import RiSet, abs_to_rel_days + + +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)}") + + +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 +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 +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 + print(f"Day: {day}") + show(rs) + +# Times from timeanddate.com +# Seattle +# Sunrise 7:40 sunset 16:18 Moonrise 23:37 Moonset 12:53 +# Sydney +# 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