mangrovedigital/tide-engine-api

1"""Harmonic tidal analysis -> the ASTRONOMICAL tide baseline.
2
3We fit a constant + linear trend + a set of standard tidal constituents
4(each as a cos/sin pair) to observed water level by least squares. The fitted
5model, evaluated at any time, IS the astronomical (uncorrected) prediction.
6
7Deriving the baseline from the gauge's own record makes it self-consistent at
8the exact target point, instead of importing a distant station's datum/phase.
9"""
10
11import math
12from . import linalg
13
14# Standard tidal constituents: speed in degrees per mean solar hour.
15CONSTITUENTS = {
16    "M2": 28.9841042,   # principal lunar semidiurnal
17    "S2": 30.0000000,   # principal solar semidiurnal
18    "N2": 28.4397295,   # larger lunar elliptic semidiurnal
19    "K2": 30.0821373,   # lunisolar semidiurnal
20    "K1": 15.0410686,   # lunar diurnal
21    "O1": 13.9430356,   # lunar diurnal
22    "P1": 14.9589314,   # solar diurnal
23    "Q1": 13.3986609,   # larger lunar elliptic diurnal
24    "M4": 57.9682084,   # shallow-water overtide of M2
25    "MS4": 58.9841042,  # shallow-water quarter-diurnal
26    "M6": 86.9523127,   # shallow-water overtide
27    "MF": 1.0980331,    # lunar fortnightly (long period)
28    "MM": 0.5443747,    # lunar monthly
29    "SA": 0.0410686,    # solar annual
30    "SSA": 0.0821373,   # solar semiannual
31}
32
33# rad/sec for each constituent
34_OMEGA = {k: deg_per_hr * math.pi / 180.0 / 3600.0
35          for k, deg_per_hr in CONSTITUENTS.items()}
36
37
38class HarmonicModel:
39    def __init__(self, t0, coeffs, names, tscale):
40        self.t0 = t0          # reference epoch (mean of training times)
41        self.coeffs = coeffs  # least-squares coefficients
42        self.names = names    # column names parallel to coeffs
43        self.tscale = tscale  # seconds, for the conditioning of the trend term
44
45    def design_row(self, epoch):
46        dt = epoch - self.t0
47        row = []
48        for name in self.names:
49            if name == "const":
50                row.append(1.0)
51            elif name == "trend":
52                row.append(dt / self.tscale)
53            else:
54                cons, kind = name.rsplit("_", 1)
55                w = _OMEGA[cons]
56                row.append(math.cos(w * dt) if kind == "cos" else math.sin(w * dt))
57        return row
58
59    def predict(self, epochs):
60        return [sum(r[i] * self.coeffs[i] for i in range(len(self.coeffs)))
61                for r in (self.design_row(e) for e in epochs)]
62
63    def amplitudes(self):
64        """Return {constituent: (amplitude, phase_deg)} for inspection."""
65        out = {}
66        idx = {n: i for i, n in enumerate(self.names)}
67        for cons in CONSTITUENTS:
68            ci, si = f"{cons}_cos", f"{cons}_sin"
69            if ci in idx and si in idx:
70                a = self.coeffs[idx[ci]]
71                b = self.coeffs[idx[si]]
72                out[cons] = (math.hypot(a, b),
73                             math.degrees(math.atan2(-b, a)) % 360.0)
74        return out
75
76
77def fit(t, y, constituents=None, ridge=1e-7, include_trend=False):
78    """Fit harmonic model to (t, y). constituents: iterable of names or None=all.
79
80    include_trend=False by default: a linear trend extrapolates unboundedly and
81    blows up the baseline out of sample. Seasonal variation is carried by the
82    periodic SA/SSA/MM/MF constituents instead, which stay bounded.
83    """
84    if constituents is None:
85        constituents = list(CONSTITUENTS.keys())
86    t0 = sum(t) / len(t)
87    tscale = max(1.0, (max(t) - min(t)) / 2.0)
88    names = ["const"]
89    if include_trend:
90        names.append("trend")
91    for c in constituents:
92        names += [f"{c}_cos", f"{c}_sin"]
93    X = []
94    for e in t:
95        dt = e - t0
96        row = [1.0]
97        if include_trend:
98            row.append(dt / tscale)
99        for c in constituents:
100            w = _OMEGA[c]
101            row.append(math.cos(w * dt))
102            row.append(math.sin(w * dt))
103        X.append(row)
104    coeffs = linalg.solve_lstsq(X, y, ridge=ridge)
105    return HarmonicModel(t0, coeffs, names, tscale)
106
107
108# --------------------------------------------------------------------------
109if __name__ == "__main__":
110    # SELFTEST: synthesise a tide from known constituents + noise and confirm
111    # the analysis recovers the astronomical signal.
112    import random
113    random.seed(1)
114    truth = {"M2": (0.6, 40.0), "S2": (0.2, 90.0), "K1": (0.15, 200.0),
115             "O1": (0.1, 150.0)}
116    t0 = 1.7e9
117    t = [t0 + i * 3600 for i in range(24 * 90)]  # 90 days hourly
118
119    def synth(e):
120        v = 0.05  # mean
121        for c, (amp, ph) in truth.items():
122            w = _OMEGA[c]
123            v += amp * math.cos(w * (e - t0) - math.radians(ph))
124        return v
125    y = [synth(e) + random.gauss(0, 0.02) for e in t]
126    model = fit(t, y)
127    pred = model.predict(t)
128    print("selftest harmonic RMSE vs noisy obs: %.4f m (noise sd=0.02)"
129          % linalg.rmse(pred, y))
130    clean = [synth(e) for e in t]
131    print("selftest harmonic RMSE vs clean tide: %.4f m" % linalg.rmse(pred, clean))
132    amps = model.amplitudes()
133    for c in truth:
134        print(f"  {c}: fit amp={amps[c][0]:.3f} (true {truth[c][0]:.3f})")
135    ok = linalg.rmse(pred, clean) < 0.01
136    print("HARMONIC SELFTEST:", "PASS" if ok else "FAIL")
137