diff --git a/astronomy/README.md b/astronomy/README.md new file mode 100644 index 0000000..e69de29 diff --git a/astronomy/sun_moon.py b/astronomy/sun_moon.py new file mode 100644 index 0000000..85e3526 --- /dev/null +++ b/astronomy/sun_moon.py @@ -0,0 +1,311 @@ +# 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 + +# Source "Astronomy on the Personal Computer" by Montenbruck and Pfleger +# ISBN 978-3-540-67221-0 + +import time +from math import sin, cos, sqrt, fabs, atan, radians, floor + +LAT = 53.29756504536339 # Local defaults +LONG = -2.102811634540558 +MOON_PHASE_LENGTH = 29.530588853 + + +def quad(ym, yz, yp): + # See Astronomy on the PC P48-49 + # finds the parabola throuh the three points (-1,ym), (0,yz), (1, yp) + # and returns the values of x where the parabola crosses zero + # (roots of the quadratic) + # and the number of roots (0, 1 or 2) within the interval [-1, 1] + nz = 0 + a = 0.5 * (ym + yp) - yz + b = 0.5 * (yp - ym) + c = yz + xe = -b / (2 * a) + ye = (a * xe + b) * xe + c + dis = b * b - 4.0 * a * c # discriminant of y=a*x^2 +bx +c + if dis > 0: # parabola has roots + dx = 0.5 * sqrt(dis) / fabs(a) + z1 = xe - dx + z2 = xe + dx + if fabs(z1) <= 1.0: + nz += 1 + if fabs(z2) <= 1.0: + nz += 1 + if z1 < -1.0: + z1 = z2 + return nz, z1, z2, ye + + +# **** 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 +# 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 does integer arithmetic and returns the same +# value regardless of the platform's epoch +def get_mjd(ndays: int = 0) -> int: # Days offset from today + 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 = 40587 # Modified Julian date of C epoch (1 Jan 70) + if time.gmtime(0)[0] == 2000: + mjepoch += 10957 + return mjepoch + days_from_epoch # Convert days from 1 Jan 70 to MJD + + +def frac(x): + return x % 1 + + +# Convert rise or set time to int. These can be None (no event). +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 +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) / MOON_PHASE_LENGTH + return phi % 1 + + +def minisun(t): + # Output 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) + # takes t centuries since J2000.0. Claimed good to 1 arcmin + p2 = 6.283185307 + 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 + + +def minimoon(t): + # takes t and returns the geocentric ra and dec + # claimed good to 5' (angle) in ra and 1' in dec + # tallies with another approximate method and with ICE for a couple of dates + + p2 = 6.283185307 + arc = 206264.8062 + 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 + + # 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) + + # 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) + + # ecliptic long and lat of Moon in rads + 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 + + +class RiSet: + def __init__(self, lat=LAT, long=LONG): # Local defaults + self.sglat = sin(radians(lat)) + self.cglat = cos(radians(lat)) + self.long = long + self.mjd = None # Current integer MJD + # Times in integer secs from midnight on current day + self._sr = None # Sunrise + self._ss = None # Sunset + self._mr = None # Moonrise + self._ms = None # Moon set + 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=0): + mjd = get_mjd(day) + if self.mjd is None or self.mjd != mjd: + spd = 86400 # Secs per day + self._t0 = ((round(time.time()) + day * spd) // spd) * spd # Midnight on target day + self.mjd = mjd + self._sr, self._ss = self.rise_set(True) # Sun + self._mr, self._ms = self.rise_set(False) # Moon + t = time.gmtime(time.time() + day * 86400) + self._phase = moonphase(*t[:4]) + + def sunrise(self, to=0): + return self._format(self._sr, to) + + def sunset(self, to=0): + return self._format(self._ss, to) + + def moonrise(self, to=0): + return self._format(self._mr, to) + + def moonset(self, to=0): + return self._format(self._ms, to) + + def moonphase(self): + return self._phase + + # ***** API end ***** + def _format(self, n, to): + if to == 0: # Default: secs since Midnight + return n + elif to == 1: # Secs since epoch + return None if n is None else n + self._t0 + # to == 3 + if n is None: + return "--:--:--" + else: + hr, tmp = divmod(n, 3600) + mi, sec = divmod(tmp, 60) + return f"{hr:02d}:{mi:02d}:{sec:02d}" + + def lmst(self, mjd): + # 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 + d = mjd - 51544.5 + t = d / 36525.0 + lst = 280.46061837 + 360.98564736629 * d + 0.000387933 * t * t - t * t * t / 38710000 + return (lst % 360) / 15.0 + self.long / 15 + + def sin_alt(self, hour, func): + # Returns the sine of the altitude of the object (moon or sun) + # at an hour relative to the current date (mjd) + mjd = self.mjd + hour / 24.0 + t = (mjd - 51544.5) / 36525.0 + dec, ra = func(t) + # hour angle of object: one hour = 15 degrees. Note lmst() uses longitude + tau = 15.0 * (self.lmst(mjd) - ra) + # sin(alt) of object using the conversion formulas + salt = self.sglat * sin(radians(dec)) + self.cglat * cos(radians(dec)) * cos(radians(tau)) + return salt + + # Modified to find sunrise and sunset only, not twilight events. + def rise_set(self, sun): + # this is my attempt to encapsulate most of the program in a function + # then this function can be generalised to find all the Sun events. + t_rise = None # Rise and set times in secs from midnight + t_set = None + sinho = sin(radians(-0.833)) if sun else sin(radians(8 / 60)) + # moonrise taken as centre of moon at +8 arcmin + # sunset upper limb simple refraction + # The loop finds the sin(alt) for sets of three consecutive + # hours, and then tests for a single zero crossing in the interval + # or for two zero crossings in an interval for for a grazing event + func = minisun if sun else minimoon + yp = self.sin_alt(0, func) - sinho + for hour in range(1, 24, 2): + ym = yp + yz = self.sin_alt(hour, func) - sinho + yp = self.sin_alt(hour + 1, func) - sinho + nz, z1, z2, ye = quad(ym, yz, yp) # Find horizon crossings + if nz == 1: # One crossing found + if ym < 0.0: + t_rise = 3600 * (hour + z1) + else: + t_set = 3600 * (hour + z1) + # case where two events are found in this interval + # (rare but whole reason we are not using simple iteration) + elif nz == 2: + if ye < 0.0: + t_rise = 3600 * (hour + z2) + t_set = 3600 * (hour + z1) + else: + t_rise = 3600 * (hour + z1) + t_set = 3600 * (hour + z2) + + 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() +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.sunrise(1)) +for d in range(30): + r.set_day(d) + print(round(r.moonphase() * 1000))