"""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

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


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 = np.polynomial.legendre.leggauss(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