137 lines
5.2 KiB
Python
137 lines
5.2 KiB
Python
"""Pure arithmetic of the Dütsch Re 100 / KC 5 protocol (numpy only, no CUDA).
|
|
|
|
docs/validation/duetsch_re100_kc5_protocol.md is the contract: elapsed time s
|
|
in native steps from the right turning point, paper phase
|
|
phi = (360 s / T - 90) mod 360, paper cycle k >= 1 spans
|
|
s/T in [0.25 + (k-1), 0.25 + k); cycle 0 is the startup quarter.
|
|
"""
|
|
from __future__ import annotations
|
|
|
|
import math
|
|
from functools import lru_cache
|
|
|
|
import numpy as np
|
|
|
|
KC, RE = 5.0, 100.0
|
|
STATIONS = (-0.6, 0.0, 0.6, 1.2) # lab-fixed x/D from the mean position
|
|
PHASES = (180.0, 210.0, 330.0) # degrees; bins [p, p + 1)
|
|
X2 = tuple(np.round(np.arange(-1.5, 1.5001, 0.05), 10)) # declared x2/D points
|
|
DISCARD, LAST_SCORED, N_PAIRS = 2, 11, 3 # cycles 0..2 discarded, cap 12 T
|
|
# A-priori bands. Settling is one fifth of the tightest comparison band of the
|
|
# same observable, and comparisons add both cases' residual window drift
|
|
# (conservative: |difference| + drift_a + drift_b <= band).
|
|
DOMAIN = dict(wave=.01, morison=.02, profile=.02) # vs the 60 D reference at D = 40
|
|
CROSS = dict(wave=.02, profile=.02) # D40 vs D80
|
|
SETTLE = dict(wave=.002, scalar=.002, profile=.004, morison=.004)
|
|
|
|
|
|
def amplitude(D):
|
|
return KC * D / (2.0 * math.pi)
|
|
|
|
|
|
def position(s, T, A):
|
|
"""Centre position relative to the mean: x = A cos(2 pi s / T), evaluated from the
|
|
native step modulo T so every cycle repeats the same geometry bit for bit."""
|
|
return A * math.cos(2.0 * math.pi * (int(s) % T) / T)
|
|
|
|
|
|
def phase(s, T):
|
|
"""Paper phase in degrees from the native step modulo T (the unwrapped cycle
|
|
index is kept separately by ``cycle``)."""
|
|
return (360.0 * (int(s) % T) / T - 90.0) % 360.0
|
|
|
|
|
|
def cycle(s, T):
|
|
return 0 if s < 0.25 * T else int(math.floor(s / T - 0.25)) + 1
|
|
|
|
|
|
def bin_phase(s, T):
|
|
"""Scored phase whose 1-degree bin [p, p + 1) contains time s, else None."""
|
|
phi = phase(s, T)
|
|
for p in PHASES:
|
|
if p <= phi < p + 1.0:
|
|
return p
|
|
return None
|
|
|
|
|
|
@lru_cache(maxsize=8)
|
|
def _quadrature_rule(points):
|
|
"""Reuse the nodes and weights, which depend only on quadrature order."""
|
|
return np.polynomial.legendre.leggauss(points)
|
|
|
|
|
|
def interval_bases(n, T, points=16):
|
|
"""Single-step interval averages over [n, n + 1] of cos phi |cos phi| and sin phi
|
|
(phi in radians = 2 pi s / T - pi / 2), Gauss-Legendre quadrature."""
|
|
x, w = _quadrature_rule(points)
|
|
s = int(n) % T + 0.5 * (x + 1.0) # periodic: identical bases in every cycle
|
|
phi = 2.0 * np.pi * s / T - 0.5 * np.pi
|
|
c = np.cos(phi)
|
|
return 0.5 * float(np.sum(w * c * np.abs(c))), 0.5 * float(np.sum(w * np.sin(phi)))
|
|
|
|
|
|
def morison_fit(f_norm, b1, b2):
|
|
"""Least squares of f_norm = c_d b1 - (pi^2 c_i / KC) b2 (paper phase, no shift)."""
|
|
a = np.column_stack([b1, -(math.pi ** 2 / KC) * np.asarray(b2)])
|
|
(c_d, c_i), *_ = np.linalg.lstsq(a, np.asarray(f_norm), rcond=None)
|
|
return float(c_d), float(c_i)
|
|
|
|
|
|
def inside_circle(x, y, centre, radius):
|
|
return (x - centre[0]) ** 2 + (y - centre[1]) ** 2 < radius * radius
|
|
|
|
|
|
def sample(field, fluid, x, y, centre, radius):
|
|
"""Profile sample: MISSING if the exact point is inside the body (checked before
|
|
any corner interpolation), else ``bilinear_fluid``."""
|
|
if inside_circle(x, y, centre, radius):
|
|
return None, (), 0.0
|
|
return bilinear_fluid(field, fluid, x, y)
|
|
|
|
|
|
def bilinear_fluid(field, fluid, x, y):
|
|
"""Bilinear value at (x, y) in cell coordinates from fluid cells only.
|
|
|
|
Returns (value, ((ix, iy, weight), ...), fluid_weight_sum) with weights
|
|
renormalized over the fluid corners and the original sum of their bilinear
|
|
weights, or (None, (), 0.0) when no corner is fluid (MISSING, never zero).
|
|
"""
|
|
x0, y0 = int(math.floor(x)), int(math.floor(y))
|
|
fx, fy = x - x0, y - y0
|
|
corners = ((x0, y0, (1 - fx) * (1 - fy)), (x0 + 1, y0, fx * (1 - fy)),
|
|
(x0, y0 + 1, (1 - fx) * fy), (x0 + 1, y0 + 1, fx * fy))
|
|
used = [(i, j, w) for i, j, w in corners if w > 0.0 and fluid[j, i]]
|
|
total = sum(w for _, _, w in used)
|
|
if total == 0.0:
|
|
return None, (), 0.0
|
|
used = tuple((i, j, w / total) for i, j, w in used)
|
|
return float(sum(w * field[j, i] for i, j, w in used)), used, float(total)
|
|
|
|
|
|
def settled_start(pairs_pass):
|
|
"""K_case from {k: pass of pair (k-1 -> k)}: first cycle after three consecutive
|
|
passing pairs among cycles > DISCARD, or None."""
|
|
run = 0
|
|
for k in sorted(pairs_pass):
|
|
if k - 1 <= DISCARD:
|
|
continue
|
|
run = run + 1 if pairs_pass[k] else 0
|
|
if run == N_PAIRS:
|
|
return k + 1
|
|
return None
|
|
|
|
|
|
def common_window(k_cases, pairs_pass_by_case, last_scored=LAST_SCORED):
|
|
"""Common scored cycles K..K+2 with K = max K_case, or None (FAIL-CLOSED) when a
|
|
case never settled, K + 2 exceeds the cap, or a case fails any pair p_j for
|
|
j = its own K_case .. K + 2 (continuous settling, including the wait for K)."""
|
|
if any(k is None for k in k_cases):
|
|
return None
|
|
K = max(k_cases)
|
|
if K + 2 > last_scored:
|
|
return None
|
|
for k_case, pairs in zip(k_cases, pairs_pass_by_case):
|
|
if not all(pairs.get(j, False) for j in range(k_case, K + 3)):
|
|
return None
|
|
return K, K + 1, K + 2
|