//! The embedded-sphere manufactured solution shared by the embedded3
//! wall gates (items 9–11): the fields, the source, the boundary data,
//! the exact surface integrals, and one steady march measured.

use rtx_cfd::solvers::incompressible::ConvectionScheme;
use rtx_cfd::solvers::incompressible::embedded3::{
    Body, FaceKind, Field, Fluid, Grid, Parameters, Solver, WallScheme,
};
use std::f64::consts::PI;

pub const RHO: f64 = 1.0;
pub const MU: f64 = 0.05;
/// The sphere's default centre (off-centre so the exact force is not zero by symmetry).
pub const C: (f64, f64, f64) = (0.6, 0.45, 0.5);
pub const R: f64 = 0.2;

pub fn u3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).cos() * (PI * z).cos()
}
pub fn v3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).cos() * (PI * y).sin() * (PI * z).cos()
}
pub fn w3(x: f64, y: f64, z: f64) -> f64 {
    -2.0 * (PI * x).cos() * (PI * y).cos() * (PI * z).sin()
}
pub fn p3(x: f64, y: f64, z: f64) -> f64 {
    (PI * x).sin() * (PI * y).sin() * (PI * z).sin()
}
/// The velocity gradient ∂u_i/∂x_j and the pressure gradient.
pub fn grads(x: f64, y: f64, z: f64) -> ([[f64; 3]; 3], [f64; 3]) {
    let (sx, cx) = (PI * x).sin_cos();
    let (sy, cy) = (PI * y).sin_cos();
    let (sz, cz) = (PI * z).sin_cos();
    (
        [
            [PI * cx * cy * cz, -PI * sx * sy * cz, -PI * sx * cy * sz],
            [-PI * sx * sy * cz, PI * cx * cy * cz, -PI * cx * sy * sz],
            [
                2.0 * PI * sx * cy * sz,
                2.0 * PI * cx * sy * sz,
                -2.0 * PI * cx * cy * cz,
            ],
        ],
        [PI * cx * sy * sz, PI * sx * cy * sz, PI * sx * sy * cz],
    )
}
pub fn source3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let (g, gp) = grads(x, y, z);
    let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
    let lap = -3.0 * PI * PI;
    let conv = |i: usize| u[0] * g[i][0] + u[1] * g[i][1] + u[2] * g[i][2];
    (
        RHO * conv(0) + gp[0] - MU * lap * u[0],
        RHO * conv(1) + gp[1] - MU * lap * u[1],
        RHO * conv(2) + gp[2] - MU * lap * u[2],
    )
}
pub fn boundary3(x: f64, y: f64, z: f64) -> (f64, f64, f64) {
    let u = if x <= 0.0 || x >= 1.0 {
        0.0
    } else {
        u3(x, y, z)
    };
    let v = if y <= 0.0 || y >= 1.0 {
        0.0
    } else {
        v3(x, y, z)
    };
    let w = if z <= 0.0 || z >= 1.0 {
        0.0
    } else {
        w3(x, y, z)
    };
    (u, v, w)
}

/// Exact force `∮ (−p I + μ(∇u + ∇uᵀ)) n dA` and momentum flux `∮ ρ u (u·n) dA`
/// over the sphere by a fine Fibonacci quadrature.
pub fn exact_force_and_flux(c: (f64, f64, f64)) -> ([f64; 3], [f64; 3]) {
    let n = 200_000;
    let golden = PI * (3.0 - 5.0_f64.sqrt());
    let (mut f, mut m) = ([0.0; 3], [0.0; 3]);
    let da = 4.0 * PI * R * R / n as f64;
    for k in 0..n {
        let zz = 1.0 - 2.0 * (k as f64 + 0.5) / n as f64;
        let rr = (1.0 - zz * zz).sqrt();
        let th = golden * k as f64;
        let nrm = [rr * th.cos(), rr * th.sin(), zz];
        let (x, y, z) = (c.0 + R * nrm[0], c.1 + R * nrm[1], c.2 + R * nrm[2]);
        let (g, _) = grads(x, y, z);
        let p = p3(x, y, z);
        let u = [u3(x, y, z), v3(x, y, z), w3(x, y, z)];
        let un = u[0] * nrm[0] + u[1] * nrm[1] + u[2] * nrm[2];
        for i in 0..3 {
            let mut t = -p * nrm[i];
            for j in 0..3 {
                t += MU * (g[i][j] + g[j][i]) * nrm[j];
            }
            f[i] += t * da;
            m[i] += RHO * u[i] * un * da;
        }
    }
    (f, m)
}

pub struct Measurement {
    pub l2_velocity: f64,
    pub max_div: f64,
    pub ghost_correction: f64,
    /// The scheme's wall route (the traction sampler on the binary wall,
    /// the operator route on the cut wall).
    pub force_surface: [f64; 3],
    pub skipped: usize,
    pub force_cv: [f64; 3],
    /// The traction sampler on the cut wall (S2-1); the wall route again
    /// on the binary wall.
    pub force_sampler: [f64; 3],
}

/// March the manufactured solution with the sphere at `c` to steady state on grid `n`.
pub fn measure(n: usize, scheme: WallScheme, c: (f64, f64, f64)) -> Measurement {
    let h = 1.0 / n as f64;
    let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h);
    let mut solver = Solver::new(
        Fluid {
            density: RHO,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        Parameters {
            corrector_steps: 2,
            tolerance: 1e-8,
            // S2-7b: `RTX_E3_MMS_SCHEME=tvd` for a second-order interior.
            convection_scheme: if std::env::var("RTX_E3_MMS_SCHEME").is_ok_and(|v| v == "tvd") {
                ConvectionScheme::TvdVanAlbada
            } else {
                ConvectionScheme::Upwind
            },
            wall_scheme: scheme,
            momentum_volume_cell_mean: std::env::var("RTX_E3_CELL_MEAN").is_ok_and(|v| v == "1"),
            ..Parameters::default()
        },
    );
    solver.set_momentum_source(|x, y, z, _t| source3(x, y, z));
    solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z));
    solver.set_body(
        Body::sphere(move |_t| c, R)
            .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))),
    );
    let g = Grid::cubic(n, n, n, h);
    let mut f = Field::new(g);
    solver.initialize(&mut f);
    let mut last = solver.advance(&mut f, dt);
    for _ in 0..200_000 {
        let (bu, bv, bw) = (f.u.clone(), f.v.clone(), f.w.clone());
        last = solver.advance(&mut f, dt);
        let mut change = 0.0_f64;
        for (a, b) in
            f.u.iter()
                .zip(&bu)
                .chain(f.v.iter().zip(&bv))
                .chain(f.w.iter().zip(&bw))
        {
            change = change.max((a - b).abs());
        }
        if change / dt < 1e-6 {
            break;
        }
    }
    let mask = solver.mask().expect("mask");
    let (mut sq, mut vol) = (0.0, 0.0);
    let dv = h * h * h;
    for k in 0..n {
        for j in 0..n {
            for i in 1..n {
                if mask.u_kind(g.uface(k, j, i)) == FaceKind::Fluid {
                    let e = f.u[g.uface(k, j, i)]
                        - u3(i as f64 * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
        for j in 1..n {
            for i in 0..n {
                if mask.v_kind(g.vface(k, j, i)) == FaceKind::Fluid {
                    let e = f.v[g.vface(k, j, i)]
                        - v3((i as f64 + 0.5) * h, j as f64 * h, (k as f64 + 0.5) * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
    }
    for k in 1..n {
        for j in 0..n {
            for i in 0..n {
                if mask.w_kind(g.wface(k, j, i)) == FaceKind::Fluid {
                    let e = f.w[g.wface(k, j, i)]
                        - w3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, k as f64 * h);
                    sq += e * e * dv;
                    vol += dv;
                }
            }
        }
    }
    let body = solver.body().expect("body");
    let t = solver.time();
    // The apertured divergence per unit volume, the porous surface's flux
    // through the wall included (the plain divergence on the binary wall);
    // a virtually merged small cell's flux counts with its master's (only
    // the pair's continuity holds).
    let mut max_div = 0.0_f64;
    let mut at_vol = 1.0;
    let mut sum_flux = 0.0;
    let (wall_fluxes, _) = mask.wall_flux_table(body, t);
    let mut cell_flux = vec![0.0; g.cells()];
    for k in 0..n {
        for j in 0..n {
            for i in 0..n {
                let idx = g.cell(k, j, i);
                if mask.is_fluid_cell(idx) {
                    let flux = (mask.a_u(g.uface(k, j, i + 1)) * f.u[g.uface(k, j, i + 1)]
                        - mask.a_u(g.uface(k, j, i)) * f.u[g.uface(k, j, i)])
                        * h
                        * h
                        + (mask.a_v(g.vface(k, j + 1, i)) * f.v[g.vface(k, j + 1, i)]
                            - mask.a_v(g.vface(k, j, i)) * f.v[g.vface(k, j, i)])
                            * h
                            * h
                        + (mask.a_w(g.wface(k + 1, j, i)) * f.w[g.wface(k + 1, j, i)]
                            - mask.a_w(g.wface(k, j, i)) * f.w[g.wface(k, j, i)])
                            * h
                            * h
                        + wall_fluxes[idx];
                    cell_flux[mask.master(idx).unwrap_or(idx)] += flux;
                }
            }
        }
    }
    for (idx, &flux) in cell_flux.iter().enumerate() {
        sum_flux += flux.abs();
        if (flux / (h * h * h)).abs() > max_div {
            max_div = (flux / (h * h * h)).abs();
            at_vol = mask.vol(idx);
        }
    }
    println!(
        "    [{scheme:?} n {n}] max div {max_div:.2e} in a cell of fluid fraction {at_vol:.3e}; Σ|flux| {sum_flux:.2e}; last step residual {:.2e}",
        last.final_residual
    );
    // S2-7b: the pressure's error against the manufactured p (mean-free
    // over the full interior cells) in the full cells and in the cut cells,
    // and the six axis wall points' reads by the linear box probe and the
    // quadratic fit (r 2.5 h) — in units of the pressure's scale (1).
    {
        let (mut sum_full, mut n_full) = (0.0, 0usize);
        let exact = |idx: usize| {
            let (k, j, i) = g.kji(idx);
            p3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h)
        };
        let is_cut = |idx: usize| mask.vol(idx) < 1.0 - 1e-9;
        for idx in 0..g.cells() {
            if mask.is_fluid_cell(idx) && !is_cut(idx) && mask.master(idx).is_none() {
                sum_full += f.p[idx] - exact(idx);
                n_full += 1;
            }
        }
        let level = sum_full / n_full.max(1) as f64;
        let (mut sq_full, mut sq_cut, mut n_cut) = (0.0, 0.0, 0usize);
        // By fluid fraction (small < 0.5 ≤ large) and the cut cells' mean
        // error (a level) against their scatter about it.
        let (mut sq_small, mut n_small, mut sq_large, mut n_large, mut sum_cut) = (0.0, 0usize, 0.0, 0usize, 0.0);
        let (mut sum_small, mut sum_large) = (0.0, 0.0);
        for idx in 0..g.cells() {
            if !mask.is_fluid_cell(idx) || mask.master(idx).is_some() {
                continue;
            }
            let e = f.p[idx] - level - exact(idx);
            if is_cut(idx) {
                sq_cut += e * e;
                n_cut += 1;
                sum_cut += e;
                if mask.vol(idx) < 0.5 {
                    sq_small += e * e;
                    sum_small += e;
                    n_small += 1;
                } else {
                    sq_large += e * e;
                    sum_large += e;
                    n_large += 1;
                }
            } else {
                sq_full += e * e;
            }
        }
        let mean_cut = sum_cut / n_cut.max(1) as f64;
        println!(
            "    [{scheme:?} n {n}] cut cells' pressure error: mean {mean_cut:+.3e}, scatter about it {:.3e}; fraction < 0.5: rms {:.3e} mean {:+.3e} ({n_small}), ≥ 0.5: rms {:.3e} mean {:+.3e} ({n_large})",
            ((sq_cut / n_cut.max(1) as f64) - mean_cut * mean_cut).max(0.0).sqrt(),
            (sq_small / n_small.max(1) as f64).sqrt(),
            sum_small / n_small.max(1) as f64,
            (sq_large / n_large.max(1) as f64).sqrt(),
            sum_large / n_large.max(1) as f64
        );
        // The six axis wall points' signed errors: linear box / quadratic
        // r 2.5 h / quadratic through full cells only (r 2.5 h and 3.5 h).
        let mut wl = String::new();
        for (dx, dy, dz) in [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, 1.0), (0.0, 0.0, -1.0)] {
            let (x, y, z) = (c.0 + R * dx, c.1 + R * dy, c.2 + R * dz);
            let pe = p3(x, y, z) + level;
            let f1 = |v: Option<f64>| v.map_or("n/a".to_string(), |v| format!("{:+.4}", v - pe));
            wl += &format!(
                " [{}]",
                [
                    f1(mask.pressure_at(&f.p, x, y, z)),
                    f1(mask.pressure_at_quadratic(&f.p, x, y, z, 2.5 * h)),
                    f1(mask.pressure_at_quadratic_from(&f.p, x, y, z, 2.5 * h, true)),
                    f1(mask.pressure_at_quadratic_from(&f.p, x, y, z, 3.5 * h, true)),
                ]
                .join(" ")
            );
        }
        println!("    [{scheme:?} n {n}] wall-point signed errors (linear, quad, full-only quad r2.5, r3.5):{wl}");
        let (mut sq_lin, mut sq_quad, mut n_lin, mut n_quad) = (0.0, 0.0, 0usize, 0usize);
        for (dx, dy, dz) in [(1.0, 0.0, 0.0), (-1.0, 0.0, 0.0), (0.0, 1.0, 0.0), (0.0, -1.0, 0.0), (0.0, 0.0, 1.0), (0.0, 0.0, -1.0)] {
            let (x, y, z) = (c.0 + R * dx, c.1 + R * dy, c.2 + R * dz);
            let pe = p3(x, y, z);
            if let Some(pl) = mask.pressure_at(&f.p, x, y, z) {
                sq_lin += (pl - level - pe).powi(2);
                n_lin += 1;
            }
            if let Some(pq) = mask.pressure_at_quadratic(&f.p, x, y, z, 2.5 * h) {
                sq_quad += (pq - level - pe).powi(2);
                n_quad += 1;
            }
        }
        println!(
            "    [{scheme:?} n {n}] pressure error rms: full cells {:.3e} ({n_full}), cut cells {:.3e} ({n_cut}); wall-point reads rms: linear box {:.3e} ({n_lin}/6), quadratic r2.5h {:.3e} ({n_quad}/6)",
            (sq_full / n_full.max(1) as f64).sqrt(),
            (sq_cut / n_cut.max(1) as f64).sqrt(),
            (sq_lin / n_lin.max(1) as f64).sqrt(),
            (sq_quad / n_quad.max(1) as f64).sqrt()
        );
    }
    let surface = match scheme {
        WallScheme::GhostBinary => mask.surface_force(body, &f, MU, t, 0.5 * h),
        WallScheme::CutCell => rtx_cfd::solvers::incompressible::embedded3::SurfaceForce {
            f: mask.cut_wall_force(body, &f, MU, t).expect("cut wall"),
            samples: 0,
            skipped: 0,
        },
    };
    let (i0, i1) = (n / 8, n - n / 8);
    let src = |x: f64, y: f64, z: f64| source3(x, y, z);
    let force_cv = mask.control_volume_force(&f, dt, RHO, MU, Some(&src), (i0, i1, i0, i1, i0, i1));
    let force_sampler = match scheme {
        WallScheme::GhostBinary => surface.f,
        WallScheme::CutCell => mask
            .cut_wall_force_reconstructed(body, &f, MU, t, None)
            .expect("reconstructed"),
    };
    Measurement {
        l2_velocity: (sq / vol).sqrt(),
        max_div,
        ghost_correction: solver.ghost_correction().abs(),
        force_surface: surface.f,
        skipped: surface.skipped,
        force_cv,
        force_sampler,
    }
}

/// S2-7b: the discrete operator applied to the EXACT manufactured field
/// valued at the faces' open-part centroids — one predictor, one
/// corrector; the predictor's acceleration per fluid face (recovered from
/// the one correction) in units of the source's largest acceleration, by
/// aperture band, plus the full faces next to cut cells and the interior;
/// the corrected field's divergence per cut cell over its largest face
/// flux; the correction's pressure at cut cells over the pressure scale.
pub fn operator_probe(n: usize, c: (f64, f64, f64)) {
    let h = 1.0 / n as f64;
    let dt = 0.4 * (h * h / (4.0 * MU / RHO)).min(h);
    let mut params = Parameters {
        corrector_steps: 1,
        tolerance: 1e-10,
        convection_scheme: if std::env::var("RTX_E3_MMS_SCHEME").is_ok_and(|v| v == "tvd") {
            ConvectionScheme::TvdVanAlbada
        } else {
            ConvectionScheme::Upwind
        },
        wall_scheme: WallScheme::CutCell,
        ..Parameters::default()
    };
    params.momentum_volume_cell_mean = false;
    let mut solver = Solver::new(
        Fluid {
            density: RHO,
            viscosity: MU,
            reference_velocity: 1.0,
            reference_length: 1.0,
        },
        params,
    );
    solver.set_momentum_source(|x, y, z, _t| source3(x, y, z));
    solver.set_boundary_velocity(|x, y, z, _t| boundary3(x, y, z));
    solver.set_body(
        Body::sphere(move |_t| c, R)
            .with_surface_velocity(|x, y, z, _t| (u3(x, y, z), v3(x, y, z), w3(x, y, z))),
    );
    let g = Grid::cubic(n, n, n, h);
    let mut f = Field::new(g);
    solver.initialize(&mut f);
    let tables = solver
        .mask()
        .expect("mask")
        .face_shift_tables()
        .expect("shift tables")
        .clone();
    for k in 0..n {
        for j in 0..n {
            for i in 0..=n {
                let idx = g.uface(k, j, i);
                let t = &tables[0][3 * idx..3 * idx + 3];
                f.u[idx] = u3(i as f64 * h + t[0], (j as f64 + 0.5) * h + t[1], (k as f64 + 0.5) * h + t[2]);
            }
        }
        for j in 0..=n {
            for i in 0..n {
                let idx = g.vface(k, j, i);
                let t = &tables[1][3 * idx..3 * idx + 3];
                f.v[idx] = v3((i as f64 + 0.5) * h + t[0], j as f64 * h + t[1], (k as f64 + 0.5) * h + t[2]);
            }
        }
    }
    for k in 0..=n {
        for j in 0..n {
            for i in 0..n {
                let idx = g.wface(k, j, i);
                let t = &tables[2][3 * idx..3 * idx + 3];
                f.w[idx] = w3((i as f64 + 0.5) * h + t[0], (j as f64 + 0.5) * h + t[1], k as f64 * h + t[2]);
            }
        }
    }
    for idx in 0..g.cells() {
        let (k, j, i) = g.kji(idx);
        f.p[idx] = p3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
    }
    {
        let (body, mask) = (solver.body().expect("body"), solver.mask().expect("mask"));
        mask.impose(body, &mut f.u, &mut f.v, &mut f.w, solver.time());
    }
    solver.advance(&mut f, dt);
    let mask = solver.mask().expect("mask");
    let pp = &f.p_prime;
    // The source's largest acceleration (the unit of the read).
    let mut a_max: f64 = 0.0;
    for k in 0..n {
        for j in 0..n {
            for i in 0..n {
                let s = source3((i as f64 + 0.5) * h, (j as f64 + 0.5) * h, (k as f64 + 0.5) * h);
                a_max = a_max.max((s.0 * s.0 + s.1 * s.1 + s.2 * s.2).sqrt() / RHO);
            }
        }
    }
    let is_cut = |idx: usize| mask.vol(idx) < 1.0 - 1e-9;
    // Bands: α < 1/16, 1/16–1/8, 1/8–1/4, ¼–½, ½–¾, ¾–1, full next to a cut
    // cell, interior. Each entry: the acceleration and the FORCE per unit
    // `ρ h A` (the acceleration times the floored fraction max(α, 0.1)).
    let bin_of = |a: f64, near: bool| -> usize {
        if a < 1.0 / 16.0 {
            0
        } else if a < 1.0 / 8.0 {
            1
        } else if a < 0.25 {
            2
        } else if a < 1.0 {
            2 + ((a * 4.0).floor() as usize).min(3)
        } else if near {
            6
        } else {
            7
        }
    };
    let mut acc: [Vec<(f64, f64)>; 8] = Default::default();
    // u faces
    for k in 0..n {
        for j in 0..n {
            for i in 1..n {
                let idx = g.uface(k, j, i);
                if mask.u_kind(idx) != FaceKind::Fluid {
                    continue;
                }
                let (cm, cp) = (g.cell(k, j, i - 1), g.cell(k, j, i));
                let star = f.u[idx] + (dt / RHO) * mask.grad_weight(0, idx) * (pp[cp] - pp[cm]) / h;
                let a = (star - f.u_old[idx]) / dt / a_max;
                acc[bin_of(mask.a_u(idx), is_cut(cm) || is_cut(cp))].push((a, a * mask.a_u(idx).max(0.1)));
            }
        }
    }
    for k in 0..n {
        for j in 1..n {
            for i in 0..n {
                let idx = g.vface(k, j, i);
                if mask.v_kind(idx) != FaceKind::Fluid {
                    continue;
                }
                let (cm, cp) = (g.cell(k, j - 1, i), g.cell(k, j, i));
                let star = f.v[idx] + (dt / RHO) * mask.grad_weight(1, idx) * (pp[cp] - pp[cm]) / h;
                let a = (star - f.v_old[idx]) / dt / a_max;
                acc[bin_of(mask.a_v(idx), is_cut(cm) || is_cut(cp))].push((a, a * mask.a_v(idx).max(0.1)));
            }
        }
    }
    for k in 1..n {
        for j in 0..n {
            for i in 0..n {
                let idx = g.wface(k, j, i);
                if mask.w_kind(idx) != FaceKind::Fluid {
                    continue;
                }
                let (cm, cp) = (g.cell(k - 1, j, i), g.cell(k, j, i));
                let star = f.w[idx] + (dt / RHO) * mask.grad_weight(2, idx) * (pp[cp] - pp[cm]) / h;
                let a = (star - f.w_old[idx]) / dt / a_max;
                acc[bin_of(mask.a_w(idx), is_cut(cm) || is_cut(cp))].push((a, a * mask.a_w(idx).max(0.1)));
            }
        }
    }
    // The correction's pressure at cut cells and interior cells, over the
    // pressure scale (1), mean-free over the interior.
    let (mut sum_int, mut n_int) = (0.0, 0usize);
    for idx in 0..g.cells() {
        if mask.is_fluid_cell(idx) && !is_cut(idx) && mask.master(idx).is_none() {
            sum_int += pp[idx];
            n_int += 1;
        }
    }
    let lvl = sum_int / n_int.max(1) as f64;
    let (mut sq_int, mut sq_cut, mut n_cut) = (0.0, 0.0, 0usize);
    for idx in 0..g.cells() {
        if !mask.is_fluid_cell(idx) || mask.master(idx).is_some() {
            continue;
        }
        let e = pp[idx] - lvl;
        if is_cut(idx) {
            sq_cut += e * e;
            n_cut += 1;
        } else {
            sq_int += e * e;
        }
    }
    let rms = |v: &[f64]| (v.iter().map(|x| x * x).sum::<f64>() / v.len().max(1) as f64).sqrt();
    let names = ["α<1/16", "1/16–1/8", "1/8–1/4", "¼–½", "½–¾", "¾–1", "full next to cut", "interior"];
    let mut line = format!("    probe n {n} (source accel {a_max:.3e}):");
    for (b, name) in names.iter().enumerate() {
        let a: Vec<f64> = acc[b].iter().map(|x| x.0).collect();
        let fo: Vec<f64> = acc[b].iter().map(|x| x.1).collect();
        line += &format!(" {name}: {} acc {:.3e} force {:.3e};", acc[b].len(), rms(&a), rms(&fo));
    }
    line += &format!(
        " p' rms interior {:.3e} cut {:.3e} ({n_cut})",
        (sq_int / n_int.max(1) as f64).sqrt(),
        (sq_cut / n_cut.max(1) as f64).sqrt()
    );
    println!("{line}");
}