//! S2-9: the recorded FSI2 kinematics for the flag tests (`RTX_E3_FLAG_KINEMATICS=<csv>`
//! from `fsi2_overset_dump_centreline`): per knot t and per station (x, y, vx, vy);
//! between knots the cubic Hermite interpolant through the knots with their
//! velocities (the overset replay's form). The run's t = 0 maps to the recorded
//! `RTX_E3_FLAG_KIN_T0` (13.0 s, inside the limit cycle) and the motion ramps in
//! over `RTX_E3_FLAG_KIN_RAMP` (0.5 s) from the undeflected line at `cy`.
#![allow(dead_code)]

fn env_f(name: &str, default: f64) -> f64 {
    std::env::var(name).ok().and_then(|v| v.parse().ok()).unwrap_or(default)
}

pub struct Recorded {
    pub times: Vec<f64>,
    /// Per knot: per station (x, y, vx, vy).
    pub knots: Vec<Vec<(f64, f64, f64, f64)>>,
    pub t0: f64,
    pub ramp: f64,
}

pub fn recorded() -> Option<&'static Recorded> {
    static REC: std::sync::OnceLock<Option<Recorded>> = std::sync::OnceLock::new();
    REC.get_or_init(|| {
        let path = std::env::var("RTX_E3_FLAG_KINEMATICS").ok()?;
        let text = std::fs::read_to_string(&path).expect("kinematics csv");
        let mut times = Vec::new();
        let mut knots = Vec::new();
        for line in text.lines().skip(1) {
            let v: Vec<f64> = line.split(',').map(|x| x.trim().parse().expect("number")).collect();
            times.push(v[0]);
            knots.push(v[1..].chunks_exact(4).map(|c| (c[0], c[1], c[2], c[3])).collect());
        }
        assert!(times.len() >= 2, "kinematics: at least two knots");
        Some(Recorded {
            times,
            knots,
            t0: env_f("RTX_E3_FLAG_KIN_T0", 13.0),
            ramp: env_f("RTX_E3_FLAG_KIN_RAMP", 0.5),
        })
    })
    .as_ref()
}

impl Recorded {
    /// The largest recorded station speed (the CFL and band bound).
    pub fn max_speed(&self) -> f64 {
        self.knots
            .iter()
            .flat_map(|k| k.iter().map(|p| (p.2 * p.2 + p.3 * p.3).sqrt()))
            .fold(0.0f64, f64::max)
    }

    /// The stations at the run's time `t`: (x, y, vx, vy), ramped from the
    /// undeflected line at `cy` (each station at its first-knot x).
    pub fn at(&self, t: f64, cy: f64) -> Vec<(f64, f64, f64, f64)> {
        let tr = (self.t0 + t).min(*self.times.last().unwrap());
        let k = match self.times.binary_search_by(|x| x.partial_cmp(&tr).unwrap()) {
            Ok(i) => i.min(self.times.len() - 2),
            Err(i) => i.saturating_sub(1).min(self.times.len() - 2),
        };
        let h = self.times[k + 1] - self.times[k];
        let s = ((tr - self.times[k]) / h).clamp(0.0, 1.0);
        let (s2, s3) = (s * s, s * s * s);
        let (h00, h10, h01, h11) = (2.0 * s3 - 3.0 * s2 + 1.0, s3 - 2.0 * s2 + s, -2.0 * s3 + 3.0 * s2, s3 - s2);
        let (d00, d10, d01, d11) = (
            (6.0 * s2 - 6.0 * s) / h,
            (3.0 * s2 - 4.0 * s + 1.0) / h,
            (-6.0 * s2 + 6.0 * s) / h,
            (3.0 * s2 - 2.0 * s) / h,
        );
        let (r, rd) = if t >= self.ramp {
            (1.0, 0.0)
        } else {
            let u = t / self.ramp;
            (u * u * (3.0 - 2.0 * u), 6.0 * u * (1.0 - u) / self.ramp)
        };
        let (a, b) = (&self.knots[k], &self.knots[k + 1]);
        let first = &self.knots[0];
        a.iter()
            .zip(b)
            .zip(first)
            .map(|((p, q), f0)| {
                let x = h00 * p.0 + h10 * h * p.2 + h01 * q.0 + h11 * h * q.2;
                let y = h00 * p.1 + h10 * h * p.3 + h01 * q.1 + h11 * h * q.3;
                let vx = d00 * p.0 + d10 * h * p.2 + d01 * q.0 + d11 * h * q.2;
                let vy = d00 * p.1 + d10 * h * p.3 + d01 * q.1 + d11 * h * q.3;
                let (xr, yr) = (f0.0, cy);
                (
                    (1.0 - r) * xr + r * x,
                    (1.0 - r) * yr + r * y,
                    r * vx + rd * (x - xr),
                    r * vy + rd * (y - yr),
                )
            })
            .collect()
    }

    /// The point and velocity at arc fraction `s ∈ [0, 1]` along the
    /// stations' polyline at `t` (linear between stations by arc length).
    pub fn along(&self, pts: &[(f64, f64, f64, f64)], s: f64) -> (f64, f64, f64, f64) {
        let mut cum = vec![0.0; pts.len()];
        for m in 1..pts.len() {
            let (dx, dy) = (pts[m].0 - pts[m - 1].0, pts[m].1 - pts[m - 1].1);
            cum[m] = cum[m - 1] + (dx * dx + dy * dy).sqrt();
        }
        let target = s.clamp(0.0, 1.0) * cum[pts.len() - 1];
        let m = cum.partition_point(|&c| c < target).clamp(1, pts.len() - 1);
        let seg = cum[m] - cum[m - 1];
        let u = if seg > 0.0 { ((target - cum[m - 1]) / seg).clamp(0.0, 1.0) } else { 0.0 };
        let (a, b) = (pts[m - 1], pts[m]);
        (a.0 + u * (b.0 - a.0), a.1 + u * (b.1 - a.1), a.2 + u * (b.2 - a.2), a.3 + u * (b.3 - a.3))
    }
}