kopia lustrzana https://github.com/peterhinch/micropython-samples
Update astronomy WIP.
Peter Hinch
rodzic
d29fb19eb4
commit
0bc573837a
|
|
@ -0,0 +1,75 @@
|
|||
Astronomy on the Personal Computer; 4th Edition
|
||||
|
||||
Springer-Verlag GmbH & Co.KG Berlin Heidelberg and the user of this CD-ROM conclude the following Agreement concerning:
|
||||
|
||||
Use and Registration
|
||||
|
||||
Art. 1 Copyright and License
|
||||
|
||||
1 The data stored on the CD-ROM, the software and the manual, hereinafter referred to as the Objects are protected by copyright; in relation to the user any rights in them lie exclusively with Springer-Verlag.
|
||||
Independently thereof, the parties hereto agree to apply the copyright principles to the Objects.
|
||||
|
||||
2 The user is, not exclusively on the basis of the law of obligations, entitled to use the Objects in the way described in the manuals. Any other ways or possibilities of using the Objects are inadmissible, in particular any translation, reproduction, decompilation, transformation in a machine-readable language and public communication; this applies to all Objects as a whole and to any of their parts.
|
||||
|
||||
3 The Objects shall at any one time be used only on one computer and at one terminal. The person using them must pertain to the user’s institution (for example, employees of his or her enterprise or library users). If the Objects are used on computers with more than one terminal or in networks, the user must obtain additional licenses from Springer-Verlag. Such a license is obtained by acquiring the corresponding multiuser or network
|
||||
version of the CD-ROM and paying the purchase price which is valid for that version at the time. In any other respect that version is subject to the same conditions of use.
|
||||
|
||||
4 The user shall store the Objects with due care in order to prevent third persons from accessing and abusing the Objects. The data stored on the Objects may not be reproduced.
|
||||
|
||||
Art. 2 Transmittal
|
||||
|
||||
1 Any transmittal (for example, sale) of the Objects to third persons, and therefore any passing on of the right and the possibility to use them, is only admissible with written authorization from Springer-Verlag.
|
||||
|
||||
2 Springer-Verlag shall provide such authorization if the previous user files a written application, supplying a representation given by the new user to the effect that he or she will comply with the provisions of the present Agreement. When the authorization is received by the previous user, the license expires, rendering the transmittal admissible.
|
||||
|
||||
|
||||
Art. 3 Consultancy
|
||||
|
||||
1 Springer-Verlag offers the possibility of asking the originator questions with respect to the use of the CD-ROM. There is, however, no legal title to the availability of this service.
|
||||
|
||||
2 The questions may refer to the installation of the program or to problems concerning its operation and use.
|
||||
|
||||
3 Questions shall be addressed in writing or via mail to Springer-Verlag (em-helpdesk@springer.de). Springer-Verlag forwards the answers given by the originator or manufacturer without verifying them. The answers will normally be given in the order in which they were received. It will not be possible to answer all questions.
|
||||
|
||||
Art. 4 Warranty
|
||||
|
||||
1 Springer-Verlag is not the originator of the data and programs but only makes them available. The user is conscious of the fact that it is not possible to create faultless databases and software; users will take appropriate steps to verify the correctness of the results of their queries.
|
||||
|
||||
2 In the case of faulty material, manufacturing defects, absence of warranted characteristics, or damage in transit, Springer-Verlag shall exchange the Object. Further claims shall only be admitted if the user has purchased the Objects from Springer-Verlag directly. The warranty requires the user to supply a detailed written description of
|
||||
any fault immediately.
|
||||
|
||||
Art. 5 Liabilities of Springer-Verlag
|
||||
|
||||
1 Springer-Verlag will only be liable for damages, whatever the legal ground, in the case of intent or gross negligence and with respect to warranted characteristics. A warranty of specific characteristics is given only in individual cases to a specific user and requires explicit written representation. Neither the originator nor the manufacturer is liable for information given under Article 4. Liability under the product liability act is not
|
||||
affected hereby. Springer-Verlag may always claim a contributory fault on the part of the user.
|
||||
|
||||
2 The originator or manufacturer named on the Objects will only be liable to the user, whatever the legal ground, in the case of intent or gross negligence.
|
||||
|
||||
Art. 6 Liabilities of the User
|
||||
|
||||
1 The user agrees to comply with the rules for use and transmittal (Articles 1 and 2). Violations of these rules may be criminal offences and may also give rise to claims for damages against the user from the licensers of Springer-Verlag.
|
||||
|
||||
2 In the case of serious violations committed by the user, Springer-Verlag may revoke the license and demand that Objects are returned.
|
||||
|
||||
Art. 7 Data Protection
|
||||
|
||||
The user consents to the computerized storage and processing of his or her personal data.
|
||||
|
||||
Art. 8 Conclusion of the Agreement
|
||||
|
||||
The user agrees that he or she shall not receive, together with this agreement, a declaration of consent given by Springer-Verlag.
|
||||
|
||||
Art. 9 Final Provisions
|
||||
|
||||
1 The present Agreement applies to Objects already delivered and still to be delivered.
|
||||
2 If any provision of the Agreement is or becomes invalid or if the Agreement is incomplete, the remainder of the Agreement is not affected. The invalid provision shall then be replaced by a legally valid provision which comes as close as possible to the invalid provision in its economic effect. The same applies to possible gaps in the Agreement.
|
||||
3 This Agreement falls under the jurisdiction of the courts at Heidelberg, if the user is a merchant who has been entered in the commercial register, a legal entity under public law, or a public special fund, or if the user has no residence or place of business in Germany.
|
||||
|
||||
4 This Agreement is subject to the laws of Germany to the exclusion of the uncitral trading rules.
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
||||
|
|
@ -8,33 +8,42 @@ 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.
|
||||
|
||||
## Licensing and acknowledgements
|
||||
|
||||
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.
|
||||
|
||||
# The RiSet class
|
||||
|
||||
## Constructor
|
||||
|
||||
Args (float):
|
||||
* `lat=LAT` Latitude in degrees. Defaults are my location. :)
|
||||
* `lat=LAT` Latitude in degrees (-ve is South). Defaults are my location. :)
|
||||
* `long=LONG` Longitude in degrees (-ve is West).
|
||||
* `lto=0` Local time offset in hours (-ve is West).
|
||||
* `lto=0` Local time offset in hours to UTC (-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.
|
||||
* `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.
|
||||
* `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.
|
||||
|
||||
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
|
||||
|
|
@ -47,6 +56,13 @@ Variants:
|
|||
(or `None`).
|
||||
* 2 Return text of form hh:mm:ss (or --:--:--) being local time.
|
||||
|
||||
Example constructor invocations:
|
||||
```python
|
||||
r = RiSet() # UK near Manchester
|
||||
r = RiSet(lat=47.609722, long=-122.3306, lto=-8) # Seattle 47°36′35″N 122°19′59″W
|
||||
r = RiSet(lat=-33.87667, long=151.21, lto=11) # Sydney 33°52′04″S 151°12′36″E
|
||||
```
|
||||
|
||||
# The moonphase function
|
||||
|
||||
This is a simple function whose provenance is uncertain. I have a lunar clock
|
||||
|
|
@ -55,11 +71,11 @@ 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
|
||||
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.
|
||||
|
|
|
|||
|
|
@ -1,8 +1,7 @@
|
|||
# sun_moon.py MicroPython Port of lunarmath.c
|
||||
# Calculate sun and moon rise and set times for any date and location
|
||||
|
||||
# Copyright (c) 2023 Peter Hinch
|
||||
# Released under the MIT License (MIT) - see LICENSE file
|
||||
# Licensing and copyright: see README.md
|
||||
|
||||
# Source "Astronomy on the Personal Computer" by Montenbruck and Pfleger
|
||||
# ISBN 978-3-540-67221-0
|
||||
|
|
@ -65,6 +64,7 @@ def get_mjd(ndays: int = 0) -> int: # Days offset from today
|
|||
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)
|
||||
mjepoch = 40587 # Modified Julian date of C epoch (1 Jan 70)
|
||||
if time.gmtime(0)[0] == 2000:
|
||||
mjepoch += 10957
|
||||
|
|
@ -98,7 +98,7 @@ def moonphase(year, month, day, hour):
|
|||
|
||||
|
||||
def minisun(t):
|
||||
# Output dec, ra
|
||||
# Output sin(dec), cos(dec), ra
|
||||
# returns the ra and dec of the Sun
|
||||
# in decimal hours, degs referred to the equinox of date and using
|
||||
# obliquity of the ecliptic at J2000.0 (small error for +- 100 yrs)
|
||||
|
|
@ -107,19 +107,17 @@ def minisun(t):
|
|||
coseps = 0.91748
|
||||
sineps = 0.39778
|
||||
|
||||
M = p2 * frac(0.993133 + 99.997361 * t)
|
||||
DL = 6893.0 * sin(M) + 72.0 * sin(2 * M)
|
||||
L = p2 * frac(0.7859453 + M / p2 + (6191.2 * t + DL) / 1296000)
|
||||
SL = sin(L)
|
||||
X = cos(L)
|
||||
Y = coseps * SL
|
||||
Z = sineps * SL
|
||||
RHO = sqrt(1 - Z * Z)
|
||||
dec = (360.0 / p2) * atan(Z / RHO)
|
||||
ra = (48.0 / p2) * atan(Y / (X + RHO))
|
||||
if ra < 0:
|
||||
ra += 24
|
||||
return dec, ra
|
||||
m = p2 * frac(0.993133 + 99.997361 * t)
|
||||
dl = 6893.0 * sin(m) + 72.0 * sin(2 * m)
|
||||
l = p2 * frac(0.7859453 + m / p2 + (6191.2 * t + dl) / 1296000)
|
||||
sl = sin(l)
|
||||
x = cos(l)
|
||||
y = coseps * sl
|
||||
z = sineps * sl
|
||||
rho = sqrt(1 - z * z)
|
||||
# dec = (360.0 / p2) * atan(z / rho)
|
||||
ra = (48.0 / p2) * atan(y / (x + rho)) % 24
|
||||
return z, rho, ra
|
||||
|
||||
|
||||
def minimoon(t):
|
||||
|
|
@ -132,56 +130,54 @@ def minimoon(t):
|
|||
coseps = 0.91748
|
||||
sineps = 0.39778
|
||||
|
||||
L0 = frac(0.606433 + 1336.855225 * t) # mean longitude of moon
|
||||
L = p2 * frac(0.374897 + 1325.552410 * t) # mean anomaly of Moon
|
||||
LS = p2 * frac(0.993133 + 99.997361 * t) # mean anomaly of Sun
|
||||
D = p2 * frac(0.827361 + 1236.853086 * t) # difference in longitude of moon and sun
|
||||
F = p2 * frac(0.259086 + 1342.227825 * t) # mean argument of latitude
|
||||
l0 = frac(0.606433 + 1336.855225 * t) # mean longitude of moon
|
||||
l = p2 * frac(0.374897 + 1325.552410 * t) # mean anomaly of Moon
|
||||
ls = p2 * frac(0.993133 + 99.997361 * t) # mean anomaly of Sun
|
||||
d = p2 * frac(0.827361 + 1236.853086 * t) # difference in longitude of moon and sun
|
||||
f = p2 * frac(0.259086 + 1342.227825 * t) # mean argument of latitude
|
||||
|
||||
# corrections to mean longitude in arcsec
|
||||
DL = 22640 * sin(L)
|
||||
DL += -4586 * sin(L - 2 * D)
|
||||
DL += +2370 * sin(2 * D)
|
||||
DL += +769 * sin(2 * L)
|
||||
DL += -668 * sin(LS)
|
||||
DL += -412 * sin(2 * F)
|
||||
DL += -212 * sin(2 * L - 2 * D)
|
||||
DL += -206 * sin(L + LS - 2 * D)
|
||||
DL += +192 * sin(L + 2 * D)
|
||||
DL += -165 * sin(LS - 2 * D)
|
||||
DL += -125 * sin(D)
|
||||
DL += -110 * sin(L + LS)
|
||||
DL += +148 * sin(L - LS)
|
||||
DL += -55 * sin(2 * F - 2 * D)
|
||||
dl = 22640 * sin(l)
|
||||
dl += -4586 * sin(l - 2 * d)
|
||||
dl += +2370 * sin(2 * d)
|
||||
dl += +769 * sin(2 * l)
|
||||
dl += -668 * sin(ls)
|
||||
dl += -412 * sin(2 * f)
|
||||
dl += -212 * sin(2 * l - 2 * d)
|
||||
dl += -206 * sin(l + ls - 2 * d)
|
||||
dl += +192 * sin(l + 2 * d)
|
||||
dl += -165 * sin(ls - 2 * d)
|
||||
dl += -125 * sin(d)
|
||||
dl += -110 * sin(l + ls)
|
||||
dl += +148 * sin(l - ls)
|
||||
dl += -55 * sin(2 * f - 2 * d)
|
||||
|
||||
# simplified form of the latitude terms
|
||||
S = F + (DL + 412 * sin(2 * F) + 541 * sin(LS)) / arc
|
||||
H = F - 2 * D
|
||||
N = -526 * sin(H)
|
||||
N += +44 * sin(L + H)
|
||||
N += -31 * sin(-L + H)
|
||||
N += -23 * sin(LS + H)
|
||||
N += +11 * sin(-LS + H)
|
||||
N += -25 * sin(-2 * L + F)
|
||||
N += +21 * sin(-L + F)
|
||||
s = f + (dl + 412 * sin(2 * f) + 541 * sin(ls)) / arc
|
||||
h = f - 2 * d
|
||||
n = -526 * sin(h)
|
||||
n += +44 * sin(l + h)
|
||||
n += -31 * sin(-l + h)
|
||||
n += -23 * sin(ls + h)
|
||||
n += +11 * sin(-ls + h)
|
||||
n += -25 * sin(-2 * l + f)
|
||||
n += +21 * sin(-l + f)
|
||||
|
||||
# ecliptic long and lat of Moon in rads
|
||||
L_moon = p2 * frac(L0 + DL / 1296000)
|
||||
B_moon = (18520.0 * sin(S) + N) / arc
|
||||
l_moon = p2 * frac(l0 + dl / 1296000)
|
||||
b_moon = (18520.0 * sin(s) + n) / arc
|
||||
|
||||
# equatorial coord conversion - note fixed obliquity
|
||||
CB = cos(B_moon)
|
||||
X = CB * cos(L_moon)
|
||||
V = CB * sin(L_moon)
|
||||
W = sin(B_moon)
|
||||
Y = coseps * V - sineps * W
|
||||
Z = sineps * V + coseps * W
|
||||
RHO = sqrt(1.0 - Z * Z)
|
||||
dec = (360.0 / p2) * atan(Z / RHO)
|
||||
ra = (48.0 / p2) * atan(Y / (X + RHO))
|
||||
if ra < 0:
|
||||
ra += 24
|
||||
return dec, ra
|
||||
cb = cos(b_moon)
|
||||
x = cb * cos(l_moon)
|
||||
v = cb * sin(l_moon)
|
||||
W = sin(b_moon)
|
||||
y = coseps * v - sineps * W
|
||||
z = sineps * v + coseps * W
|
||||
rho = sqrt(1.0 - z * z)
|
||||
# dec = (360.0 / p2) * atan(z / rho)
|
||||
ra = (48.0 / p2) * atan(y / (x + rho)) % 24
|
||||
return z, rho, ra
|
||||
|
||||
|
||||
class RiSet:
|
||||
|
|
@ -197,29 +193,19 @@ class RiSet:
|
|||
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: int = 0):
|
||||
# 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):
|
||||
mjd = get_mjd(day)
|
||||
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
|
||||
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
|
||||
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()
|
||||
|
||||
# variants: 0 secs since 00:00:00 localtime. 1 secs since MicroPython epoch
|
||||
|
|
@ -236,10 +222,29 @@ class RiSet:
|
|||
def moonset(self, variant: int = 0):
|
||||
return self._format(self._times[3], variant)
|
||||
|
||||
def moonphase(self):
|
||||
def moonphase(self) -> float:
|
||||
return self._phase
|
||||
|
||||
def set_lto(self, t): # Update the offset from UTC
|
||||
pass # TODO
|
||||
|
||||
# ***** API end *****
|
||||
# Re-calculate rise and set times
|
||||
def update(self, mjd):
|
||||
self._times = [None] * 4
|
||||
days = (1, 2) if self.lto < 0 else (1,) if self.lto == 0 else (0, 1)
|
||||
for day in days:
|
||||
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
|
||||
|
||||
def adjust(self, n, day, idx):
|
||||
if n is not None:
|
||||
n += self.lto + (day - 1) * 86400
|
||||
|
|
@ -261,7 +266,7 @@ class RiSet:
|
|||
return f"{hr:02d}:{mi:02d}:{sec:02d}"
|
||||
|
||||
# See https://github.com/orgs/micropython/discussions/13075
|
||||
def lmstt(self, t):
|
||||
def lmst(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
|
||||
|
|
@ -278,13 +283,13 @@ class RiSet:
|
|||
# Returns the sine of the altitude of the object (moon or sun)
|
||||
# at an hour relative to the current date (mjd)
|
||||
mjd1 = (self.mjd - 51544.5) + hour / 24.0
|
||||
# mjd1 = self.mjd + 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)
|
||||
sin_dec, cos_dec, ra = func(t1)
|
||||
# hour angle of object: one hour = 15 degrees. Note lmst() uses longitude
|
||||
tau = 15.0 * (self.lmstt(t1) - ra)
|
||||
tau = 15.0 * (self.lmst(t1) - ra)
|
||||
# sin(alt) of object using the conversion formulas
|
||||
return self.sglat * sin(radians(dec)) + self.cglat * cos(radians(dec)) * cos(radians(tau))
|
||||
return self.sglat * sin_dec + self.cglat * cos_dec * cos(radians(tau))
|
||||
|
||||
# Modified to find sunrise and sunset only, not twilight events.
|
||||
def rise_set(self, sun):
|
||||
|
|
@ -325,8 +330,9 @@ class RiSet:
|
|||
return to_int(t_rise), to_int(t_set) # Convert to int preserving None values
|
||||
|
||||
|
||||
r = RiSet()
|
||||
# Seattle RiSet(lat=47.61, long=-122.35, lto=-8)
|
||||
# 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))
|
||||
|
|
|
|||
Ładowanie…
Reference in New Issue