9756538f16
Organizes 11 projects for Cerbo GX/Venus OS into a single repository: - axiom-nmea: Raymarine LightHouse protocol decoder - dbus-generator-ramp: Generator current ramp controller - dbus-lightning: Blitzortung lightning monitor - dbus-meteoblue-forecast: Meteoblue weather forecast - dbus-no-foreign-land: noforeignland.com tracking - dbus-tides: Tide prediction from depth + harmonics - dbus-vrm-history: VRM cloud history proxy - dbus-windy-station: Windy.com weather upload - mfd-custom-app: MFD app deployment package - venus-html5-app: Custom Victron HTML5 app fork - watermaker: Watermaker PLC control UI Adds root README, .gitignore, project template, and per-project .gitignore files. Sensitive config files excluded via .gitignore with .example templates provided. Made-with: Cursor
385 lines
12 KiB
Python
385 lines
12 KiB
Python
"""
|
|
Pure-Python tidal harmonic prediction engine.
|
|
|
|
Implements the IHO standard harmonic method for predicting ocean tides.
|
|
Zero external dependencies -- uses only the Python standard library.
|
|
|
|
Reference: Schureman, "Manual of Harmonic Analysis and Prediction of Tides"
|
|
(US Coast and Geodetic Survey, Special Publication No. 98)
|
|
|
|
The prediction equation:
|
|
h(t) = Z0 + sum_i( f_i * A_i * cos(V_i(t) + u_i - G_i) )
|
|
|
|
Where:
|
|
Z0 = mean water level offset
|
|
A_i = amplitude of constituent i (from model, location-specific)
|
|
G_i = phase (Greenwich epoch) of constituent i (from model)
|
|
f_i = nodal amplitude factor (varies over 18.6-year lunar node cycle)
|
|
u_i = nodal phase correction (same cycle)
|
|
V_i = astronomical argument at time t
|
|
"""
|
|
|
|
import math
|
|
|
|
DEG2RAD = math.pi / 180.0
|
|
RAD2DEG = 180.0 / math.pi
|
|
|
|
# =====================================================================
|
|
# TIDAL CONSTITUENT DEFINITIONS
|
|
# =====================================================================
|
|
#
|
|
# 37 standard tidal constituents.
|
|
# Each: (name, Doodson_numbers, speed_deg_per_hour)
|
|
# Doodson numbers: (T, s, h, p, N_prime, ps)
|
|
# Speeds from Schureman Table 2.
|
|
|
|
CONSTITUENTS = [
|
|
# Semi-diurnal
|
|
("M2", (2, -2, 2, 0, 0, 0), 28.9841042),
|
|
("S2", (2, 0, 0, 0, 0, 0), 30.0000000),
|
|
("N2", (2, -3, 2, 1, 0, 0), 28.4397295),
|
|
("K2", (2, 0, 2, 0, 0, 0), 30.0821373),
|
|
("2N2", (2, -4, 2, 2, 0, 0), 27.8953548),
|
|
("MU2", (2, -4, 4, 0, 0, 0), 27.9682084),
|
|
("NU2", (2, -3, 4, -1, 0, 0), 28.5125831),
|
|
("L2", (2, -1, 2, -1, 0, 0), 29.5284789),
|
|
("T2", (2, 0, -1, 0, 0, 1), 29.9589333),
|
|
("R2", (2, 0, 1, 0, 0, -1), 30.0410667),
|
|
("LAM2", (2, -1, 0, 1, 0, 0), 29.4556253),
|
|
# Diurnal
|
|
("K1", (1, 0, 1, 0, 0, 0), 15.0410686),
|
|
("O1", (1, -2, 1, 0, 0, 0), 13.9430356),
|
|
("P1", (1, 0, -1, 0, 0, 0), 14.9589314),
|
|
("Q1", (1, -3, 1, 1, 0, 0), 13.3986609),
|
|
("J1", (1, 1, 1, -1, 0, 0), 15.5854433),
|
|
("OO1", (1, 2, 1, 0, 0, 0), 16.1391017),
|
|
("M1", (1, -1, 1, 0, 0, 0), 14.4920521),
|
|
("2Q1", (1, -4, 1, 2, 0, 0), 12.8542862),
|
|
("RHO1", (1, -3, 3, -1, 0, 0), 13.4715145),
|
|
# Long-period
|
|
("MF", (0, 2, 0, 0, 0, 0), 1.0980331),
|
|
("MM", (0, 1, 0, -1, 0, 0), 0.5443747),
|
|
("SSA", (0, 0, 2, 0, 0, 0), 0.0821373),
|
|
("MSM", (0, 1, -2, 1, 0, 0), 0.4715211),
|
|
("MS", (0, 2, -2, 0, 0, 0), 1.0158958),
|
|
# Shallow-water / compound
|
|
("M3", (3, -3, 3, 0, 0, 0), 43.4761563),
|
|
("M4", (4, -4, 4, 0, 0, 0), 57.9682084),
|
|
("MS4", (4, -2, 2, 0, 0, 0), 58.9841042),
|
|
("M6", (6, -6, 6, 0, 0, 0), 86.9523127),
|
|
("M8", (8, -8, 8, 0, 0, 0), 115.9364168),
|
|
("MN4", (4, -5, 4, 1, 0, 0), 57.4238337),
|
|
("S4", (4, 0, 0, 0, 0, 0), 60.0000000),
|
|
("S6", (6, 0, 0, 0, 0, 0), 90.0000000),
|
|
("2SM2", (2, 2, -2, 0, 0, 0), 31.0158958),
|
|
("MKS2", (2, -2, 4, 0, 0, 0), 29.0662415),
|
|
# Additional
|
|
("S1", (1, 0, 0, 0, 0, 0), 15.0000000),
|
|
("SA", (0, 0, 1, 0, 0, 0), 0.0410686),
|
|
# Terdiurnal
|
|
("MK3", (3, -2, 3, 0, 0, 0), 44.0251729),
|
|
("2MK3", (3, -4, 3, 0, 0, 0), 42.9271398),
|
|
]
|
|
|
|
_CONST_BY_NAME = {c[0]: (c[1], c[2]) for c in CONSTITUENTS}
|
|
CONSTITUENT_NAMES = [c[0] for c in CONSTITUENTS]
|
|
|
|
|
|
# =====================================================================
|
|
# ASTRONOMICAL ARGUMENTS
|
|
# =====================================================================
|
|
|
|
def _julian_centuries(unix_ts):
|
|
"""Convert Unix timestamp to Julian centuries from J2000.0."""
|
|
jd = unix_ts / 86400.0 + 2440587.5
|
|
return (jd - 2451545.0) / 36525.0
|
|
|
|
|
|
def _astro_angles(T):
|
|
"""Compute fundamental astronomical angles in degrees.
|
|
|
|
Args:
|
|
T: Julian centuries from J2000.0
|
|
|
|
Returns:
|
|
(s, h, p, N, ps) all in degrees
|
|
|
|
s = mean longitude of the Moon
|
|
h = mean longitude of the Sun
|
|
p = mean longitude of lunar perigee
|
|
N = longitude of Moon's ascending node
|
|
ps = mean longitude of solar perigee
|
|
|
|
Formulae from Meeus (1998), consistent with Schureman/IHO.
|
|
"""
|
|
s = (218.3164477
|
|
+ 481267.88123421 * T
|
|
- 0.0015786 * T * T
|
|
+ T * T * T / 538841.0
|
|
- T * T * T * T / 65194000.0)
|
|
|
|
h = 280.46646 + 36000.76983 * T + 0.0003032 * T * T
|
|
|
|
p = (83.3532465
|
|
+ 4069.0137287 * T
|
|
- 0.0103200 * T * T
|
|
- T * T * T / 80053.0
|
|
+ T * T * T * T / 18999000.0)
|
|
|
|
N = (125.04452
|
|
- 1934.13626 * T
|
|
+ 0.0020708 * T * T
|
|
+ T * T * T / 450000.0)
|
|
|
|
ps = 282.93735 + 1.71946 * T + 0.00046 * T * T
|
|
|
|
return s, h, p, N, ps
|
|
|
|
|
|
def _gmst_degrees(T):
|
|
"""Greenwich Mean Sidereal Time in degrees.
|
|
|
|
Args:
|
|
T: Julian centuries from J2000.0
|
|
|
|
Meeus (1998) eq. 12.4, valid for the 21st century.
|
|
"""
|
|
D = T * 36525.0
|
|
return (280.46061837 + 360.98564736629 * D
|
|
+ 0.000387933 * T * T
|
|
- T * T * T / 38710000.0) % 360.0
|
|
|
|
|
|
def _solar_hour_angle(theta_g, h):
|
|
"""Mean solar hour angle at Greenwich (degrees).
|
|
|
|
T = GMST - h, where h is mean longitude of the Sun.
|
|
The first Doodson number multiplies T (not lunar time tau).
|
|
Verified: speed(M2) = 2*(GMST_rate - h_rate) - 2*s_rate + 2*h_rate
|
|
= 2*GMST_rate - 2*s_rate = 28.984 deg/hr ✓
|
|
"""
|
|
return theta_g - h
|
|
|
|
|
|
def _V0(doodson, T_solar, s, h, p, N, ps):
|
|
"""Astronomical argument V0 for a constituent (degrees)."""
|
|
dT, ds, dh, dp, dN, dps = doodson
|
|
return dT * T_solar + ds * s + dh * h + dp * p + dN * N + dps * ps
|
|
|
|
|
|
# =====================================================================
|
|
# NODAL CORRECTIONS (Schureman 1958)
|
|
# =====================================================================
|
|
|
|
def _nodal_corrections(N_deg):
|
|
"""Compute nodal f (amplitude factor) and u (phase correction, degrees).
|
|
|
|
Args:
|
|
N_deg: longitude of Moon's ascending node in degrees
|
|
|
|
Returns:
|
|
dict {constituent_name: (f, u_degrees)}
|
|
"""
|
|
N = N_deg * DEG2RAD
|
|
cosN = math.cos(N)
|
|
cos2N = math.cos(2.0 * N)
|
|
sinN = math.sin(N)
|
|
sin2N = math.sin(2.0 * N)
|
|
|
|
corrections = {}
|
|
|
|
# M2 family (semi-diurnal lunar)
|
|
f_m2 = 1.0 - 0.03690 * cosN + 0.00256 * cos2N
|
|
u_m2 = -2.14 * sinN + 0.29 * sin2N
|
|
for name in ("M2", "N2", "2N2", "MU2", "NU2", "LAM2", "2SM2", "MKS2"):
|
|
corrections[name] = (f_m2, u_m2)
|
|
|
|
# S2 family (solar, no nodal dependence)
|
|
for name in ("S2", "T2", "R2"):
|
|
corrections[name] = (1.0, 0.0)
|
|
|
|
# K2
|
|
f_k2 = 1.0 + 0.2852 * cosN + 0.0324 * cos2N
|
|
u_k2 = -17.74 * sinN + 0.68 * sin2N
|
|
corrections["K2"] = (f_k2, u_k2)
|
|
|
|
# L2
|
|
lr = 1.0 - 0.25 * cosN
|
|
li = 0.25 * sinN
|
|
f_l2 = math.sqrt(lr * lr + li * li)
|
|
u_l2 = math.atan2(li, lr) * RAD2DEG
|
|
corrections["L2"] = (f_l2, u_l2)
|
|
|
|
# K1
|
|
kr = 1.0 + 0.1158 * cosN - 0.0029 * cos2N
|
|
ki = 0.1554 * sinN - 0.0029 * sin2N
|
|
f_k1 = math.sqrt(kr * kr + ki * ki)
|
|
u_k1 = -math.atan2(ki, kr) * RAD2DEG
|
|
corrections["K1"] = (f_k1, u_k1)
|
|
corrections["J1"] = (f_k1, u_k1)
|
|
corrections["OO1"] = (f_k1, u_k1)
|
|
|
|
# O1 family
|
|
o_r = 1.0 + 0.1886 * cosN - 0.0058 * cos2N
|
|
o_i = 0.1886 * sinN - 0.0058 * sin2N
|
|
f_o1 = math.sqrt(o_r * o_r + o_i * o_i)
|
|
u_o1 = math.atan2(o_i, o_r) * RAD2DEG
|
|
for name in ("O1", "Q1", "2Q1", "RHO1", "M1"):
|
|
corrections[name] = (f_o1, u_o1)
|
|
|
|
# P1, S1 -- no nodal
|
|
corrections["P1"] = (1.0, 0.0)
|
|
corrections["S1"] = (1.0, 0.0)
|
|
|
|
# Long-period
|
|
f_mf = 1.043 + 0.414 * cosN
|
|
u_mf = -23.74 * sinN + 2.68 * sin2N
|
|
corrections["MF"] = (f_mf, u_mf)
|
|
|
|
f_mm = 1.0 - 0.130 * cosN
|
|
corrections["MM"] = (f_mm, 0.0)
|
|
corrections["MSM"] = (f_mm, 0.0)
|
|
corrections["MS"] = (f_m2, u_m2)
|
|
|
|
for name in ("SSA", "SA"):
|
|
corrections[name] = (1.0, 0.0)
|
|
|
|
# Shallow water
|
|
corrections["M3"] = (f_m2 ** 1.5, 1.5 * u_m2)
|
|
corrections["M4"] = (f_m2 * f_m2, 2.0 * u_m2)
|
|
corrections["MN4"] = (f_m2 * f_m2, 2.0 * u_m2)
|
|
corrections["MS4"] = (f_m2, u_m2)
|
|
corrections["M6"] = (f_m2 ** 3, 3.0 * u_m2)
|
|
corrections["M8"] = (f_m2 ** 4, 4.0 * u_m2)
|
|
corrections["S4"] = (1.0, 0.0)
|
|
corrections["S6"] = (1.0, 0.0)
|
|
|
|
# Terdiurnal
|
|
corrections["MK3"] = (f_m2 * f_k1, u_m2 + u_k1)
|
|
corrections["2MK3"] = (f_m2 * f_m2 * f_k1, 2.0 * u_m2 + u_k1)
|
|
|
|
return corrections
|
|
|
|
|
|
# =====================================================================
|
|
# PREDICTION
|
|
# =====================================================================
|
|
|
|
def predict(constituents, timestamps, z0=0.0):
|
|
"""Predict tide heights at the given timestamps.
|
|
|
|
Args:
|
|
constituents: list of dicts with keys:
|
|
"name" - constituent name (e.g. "M2")
|
|
"amplitude" - amplitude in meters
|
|
"phase" - Greenwich phase in degrees
|
|
timestamps: list of Unix timestamps (float, seconds since epoch)
|
|
z0: mean water level offset (meters)
|
|
|
|
Returns:
|
|
list of predicted tide heights (meters)
|
|
"""
|
|
if not constituents or not timestamps:
|
|
return []
|
|
|
|
t0 = timestamps[0]
|
|
T0 = _julian_centuries(t0)
|
|
s, h, p, N, ps = _astro_angles(T0)
|
|
theta_g = _gmst_degrees(T0)
|
|
T_solar = _solar_hour_angle(theta_g, h)
|
|
nodal = _nodal_corrections(N)
|
|
|
|
precomp = []
|
|
for c in constituents:
|
|
name = c["name"]
|
|
amp = c["amplitude"]
|
|
phase_g = c["phase"]
|
|
|
|
if name not in _CONST_BY_NAME:
|
|
continue
|
|
if amp <= 0.0:
|
|
continue
|
|
|
|
doodson, speed = _CONST_BY_NAME[name]
|
|
f, u = nodal.get(name, (1.0, 0.0))
|
|
|
|
v0 = _V0(doodson, T_solar, s, h, p, N, ps)
|
|
speed_rad = speed * DEG2RAD
|
|
phase_offset = (v0 + u - phase_g) * DEG2RAD
|
|
|
|
precomp.append((amp * f, speed_rad, phase_offset))
|
|
|
|
heights = []
|
|
for ts in timestamps:
|
|
dt_hr = (ts - t0) / 3600.0
|
|
total = z0
|
|
for amp_f, speed_rad, phase_offset in precomp:
|
|
total += amp_f * math.cos(speed_rad * dt_hr + phase_offset)
|
|
heights.append(total)
|
|
|
|
return heights
|
|
|
|
|
|
def find_extrema(timestamps, heights, min_amplitude=0.3):
|
|
"""Find high and low tide events from a predicted curve.
|
|
|
|
Uses a multi-pass algorithm to avoid the cascading-rejection problem
|
|
that occurs with single-pass filtering when diurnal inequality creates
|
|
small-range tidal cycles (common in the Caribbean, Gulf of Mexico, etc.).
|
|
|
|
Pipeline:
|
|
1. Identify all local maxima and minima
|
|
2. Merge consecutive same-type extrema (enforce strict alternation)
|
|
3. Iteratively remove the adjacent pair with the smallest height
|
|
range until all remaining pairs exceed min_amplitude
|
|
|
|
Args:
|
|
timestamps: list of Unix timestamps
|
|
heights: corresponding predicted heights
|
|
min_amplitude: minimum height change between adjacent extrema
|
|
|
|
Returns:
|
|
list of (timestamp, height, "high"|"low") tuples
|
|
"""
|
|
if len(timestamps) < 3:
|
|
return []
|
|
|
|
candidates = []
|
|
for i in range(1, len(heights) - 1):
|
|
hp = heights[i - 1]
|
|
hc = heights[i]
|
|
hn = heights[i + 1]
|
|
|
|
if hc > hp and hc > hn:
|
|
candidates.append((timestamps[i], hc, "high"))
|
|
elif hc < hp and hc < hn:
|
|
candidates.append((timestamps[i], hc, "low"))
|
|
|
|
if not candidates:
|
|
return []
|
|
|
|
merged = [candidates[0]]
|
|
for ts, h, t in candidates[1:]:
|
|
if merged[-1][2] == t:
|
|
_, prev_h, _ = merged[-1]
|
|
if (t == "high" and h > prev_h) or (t == "low" and h < prev_h):
|
|
merged[-1] = (ts, h, t)
|
|
else:
|
|
merged.append((ts, h, t))
|
|
|
|
changed = True
|
|
while changed and len(merged) > 2:
|
|
changed = False
|
|
min_range = float('inf')
|
|
min_idx = -1
|
|
for i in range(len(merged) - 1):
|
|
r = abs(merged[i][1] - merged[i + 1][1])
|
|
if r < min_range:
|
|
min_range = r
|
|
min_idx = i
|
|
if min_range < min_amplitude and min_idx >= 0:
|
|
merged = merged[:min_idx] + merged[min_idx + 2:]
|
|
changed = True
|
|
|
|
return merged
|