//! Axisymmetric AQUAL field and quantum interface response. No dependencies.
//! Units: acceleration=the model scale, length=sqrt(G M_sun / scale).
use std::{env, fs, io::Write, time::Instant};

#[derive(Clone, Copy, Debug)]
enum Model {
    Newton,
    Standard,
    Linear,
    Tanh,
    Table,
}
impl Model {
    fn parse(s: &str) -> Self {
        match s {
            "newton" => Self::Newton,
            "standard" => Self::Standard,
            "linear" => Self::Linear,
            "tanh" => Self::Tanh,
            "table" => Self::Table,
            _ => panic!("Unknown model"),
        }
    }
    fn scale(self) -> f64 {
        match self {
            Self::Newton => 1.2e-10,
            Self::Standard => 1.622565958040207e-10,
            Self::Linear => 9.967065795357357e-11,
            Self::Tanh => 1.4217448781905785e-10,
            Self::Table => panic!("A tabulated response requires an explicit frozen scale"),
        }
    }
}

fn gauss_legendre(n: usize) -> (Vec<f64>, Vec<f64>) {
    let mut xs = vec![0.; n];
    let mut ws = vec![0.; n];
    for i in 0..n.div_ceil(2) {
        let mut z = (std::f64::consts::PI * (i as f64 + 0.75) / (n as f64 + 0.5)).cos();
        let mut dp = 0.;
        for _ in 0..30 {
            let (mut p0, mut p1) = (1., z);
            for j in 2..=n {
                let p = ((2 * j - 1) as f64 * z * p1 - (j - 1) as f64 * p0) / j as f64;
                p0 = p1;
                p1 = p;
            }
            dp = n as f64 * (z * p1 - p0) / (z * z - 1.);
            let next = z - p1 / dp;
            if (next - z).abs() < 2e-15 {
                z = next;
                break;
            }
            z = next;
        }
        let weight = 1. / ((1. - z * z) * dp * dp);
        xs[i] = (1. - z) / 2.;
        xs[n - 1 - i] = (1. + z) / 2.;
        ws[i] = weight;
        ws[n - 1 - i] = weight;
    }
    (xs, ws)
}

fn exact_mu(model: Model, x: f64, nodes: &[f64], weights: &[f64]) -> f64 {
    match model {
        Model::Table => panic!("Use the physical response table, not an analytic placeholder"),
        Model::Newton => 1.,
        Model::Standard => x / x.hypot(1.),
        Model::Linear => {
            if x < 1e-3 {
                return 0.75 * x - 0.3 * x * x + 3. * x.powi(4) / 56.;
            }
            1. - 3.
                * nodes
                    .iter()
                    .zip(weights)
                    .map(|(&u, &w)| w * u * u / ((x * u).hypot(1.) + x * u))
                    .sum::<f64>()
        }
        Model::Tanh => {
            if x < 1e-3 {
                return 0.75 * x + 0.5 * x * x - (5. / 6.) * x.powi(3);
            }
            let upper = x.min(24.);
            let value = nodes
                .iter()
                .zip(weights)
                .map(|(&u, &w)| {
                    let t = upper * u;
                    let f = t.tanh();
                    w * t * f / (t.cosh().powi(2) * (f.hypot(1.) + f))
                })
                .sum::<f64>()
                * upper;
            1. - 3. * value / x.powi(3)
        }
    }
}

struct Response {
    values: Vec<f64>,
    low: f64,
    step: f64,
    model: Model,
    floor: f64,
    low_slope: f64,
}
impl Response {
    fn new(model: Model, floor: f64) -> Self {
        let (nodes, weights) = gauss_legendre(96);
        let low = -25.;
        let high = 25.;
        let n = 10001;
        let step = (high - low) / (n - 1) as f64;
        let values = (0..n)
            .map(|i| exact_mu(model, (low + i as f64 * step).exp(), &nodes, &weights))
            .collect();
        Self {
            values,
            low,
            step,
            model,
            floor,
            low_slope: if matches!(model, Model::Standard | Model::Newton) {
                1.
            } else {
                0.75
            },
        }
    }
    fn from_table(csv: &str, floor: f64) -> Result<Self, String> {
        let mut lines = csv.lines();
        if lines.next().map(str::trim) != Some("x,mu") {
            return Err("Response table must have header x,mu".into());
        }
        let rows: Vec<(f64, f64)> = lines
            .filter(|line| !line.trim().is_empty())
            .map(|line| {
                let columns: Vec<&str> = line.trim().split(',').collect();
                if columns.len() != 2 {
                    return Err("Response row needs exactly two columns".into());
                }
                let x: f64 = columns[0]
                    .parse()
                    .map_err(|_| "Invalid response acceleration")?;
                let mu: f64 = columns[1].parse().map_err(|_| "Invalid response value")?;
                if !x.is_finite() || !mu.is_finite() || x <= 0. || mu <= 0. {
                    return Err(
                        "Response coordinates and values must be finite and positive".into(),
                    );
                }
                Ok((x, mu))
            })
            .collect::<Result<_, String>>()?;
        if rows.len() < 3 || !floor.is_finite() || floor < 0. {
            return Err("Need at least three response rows and a nonnegative floor".into());
        }
        let low = rows[0].0.ln();
        let step = (rows.last().unwrap().0.ln() - low) / (rows.len() - 1) as f64;
        if step <= 0. {
            return Err("Response accelerations must be increasing".into());
        }
        for (i, &(x, mu)) in rows.iter().enumerate() {
            if (x.ln() - low - i as f64 * step).abs() > 1e-9 {
                return Err("Response table must have a regular logarithmic mesh".into());
            }
            if i > 0 && x * mu <= rows[i - 1].0 * rows[i - 1].1 {
                return Err("The tabulated static map x*mu must be increasing".into());
            }
            if i > 0 && mu.min(rows[i - 1].1) + (mu - rows[i - 1].1) / step <= 0. {
                return Err(
                    "Log-linear interpolation would violate positive radial stiffness".into(),
                );
            }
        }
        if (rows.last().unwrap().1 - 1.).abs() > 1e-8 {
            return Err("Response table must reach the Newtonian high-field limit".into());
        }
        Ok(Self {
            values: rows.iter().map(|row| row.1).collect(),
            low,
            step,
            model: Model::Table,
            floor,
            low_slope: rows[0].1 / rows[0].0,
        })
    }
    fn get(&self, x: f64) -> f64 {
        if matches!(self.model, Model::Newton) {
            return 1.;
        }
        if x <= self.low.exp() {
            return (self.low_slope * x).max(self.floor);
        }
        let f = (x.ln() - self.low) / self.step;
        let i = f as usize;
        if i >= self.values.len() - 1 {
            return 1.;
        }
        let t = f - i as f64;
        ((1. - t) * self.values[i] + t * self.values[i + 1]).max(self.floor)
    }
    fn log_slope(&self, x: f64) -> f64 {
        let h = 0.01;
        (self.get(x * f64::exp(h)).ln() - self.get(x * f64::exp(-h)).ln()) / (2. * h)
    }
}

struct Grid {
    nr: usize,
    nu: usize,
    rf: Vec<f64>,
    r: Vec<f64>,
    dr: Vec<f64>,
    u: Vec<f64>,
    du: f64,
    eta: f64,
    boundary: Vec<f64>,
}
impl Grid {
    fn new(nr: usize, nu: usize, rmin: f64, rmax: f64, eta: f64, response: &Response) -> Self {
        let rf: Vec<f64> = (0..=nr)
            .map(|i| rmin * (rmax / rmin).powf(i as f64 / nr as f64))
            .collect();
        let r: Vec<f64> = (0..nr).map(|i| (rf[i] + rf[i + 1]) / 2.).collect();
        let dr: Vec<f64> = (0..nr).map(|i| rf[i + 1] - rf[i]).collect();
        let du = 2. / nu as f64;
        let u: Vec<f64> = (0..nu).map(|j| -1. + (j as f64 + 0.5) * du).collect();
        let mue = response.get(eta);
        let le = response.log_slope(eta);
        let boundary = u
            .iter()
            .map(|&u| (1. - 1. / (mue * (1. + le * (1. - u * u)).sqrt())) / rmax)
            .collect();
        Self {
            nr,
            nu,
            rf,
            r,
            dr,
            u,
            du,
            eta,
            boundary,
        }
    }
    fn id(&self, i: usize, j: usize) -> usize {
        i * self.nu + j
    }
}

struct Operator {
    diag: Vec<f64>,
    radial: Vec<f64>,
    angular: Vec<f64>,
    rhs: Vec<f64>,
    nu: usize,
}
impl Operator {
    fn apply(&self, x: &[f64], out: &mut [f64]) {
        for k in 0..x.len() {
            let mut a = self.diag[k] * x[k];
            if k + self.nu < x.len() {
                a -= self.radial[k] * x[k + self.nu];
            }
            if k >= self.nu {
                a -= self.radial[k - self.nu] * x[k - self.nu];
            }
            if k % self.nu + 1 < self.nu {
                a -= self.angular[k] * x[k + 1];
            }
            if k % self.nu > 0 {
                a -= self.angular[k - 1] * x[k - 1];
            }
            out[k] = a;
        }
    }
}

fn assemble(g: &Grid, response: &Response, w: &[f64]) -> Operator {
    let n = w.len();
    let mut wr = vec![0.; n];
    let mut wu = vec![0.; n];
    for i in 0..g.nr {
        for j in 0..g.nu {
            let k = g.id(i, j);
            wr[k] = if i == 0 {
                (w[g.id(i + 1, j)] - w[k]) / (g.r[i + 1] - g.r[i])
            } else if i + 1 == g.nr {
                (g.boundary[j] - w[g.id(i - 1, j)]) / (g.rf[g.nr] - g.r[i - 1])
            } else {
                (w[g.id(i + 1, j)] - w[g.id(i - 1, j)]) / (g.r[i + 1] - g.r[i - 1])
            };
            wu[k] = if j == 0 {
                (w[k + 1] - w[k]) / g.du
            } else if j + 1 == g.nu {
                (w[k] - w[k - 1]) / g.du
            } else {
                (w[k + 1] - w[k - 1]) / (2. * g.du)
            };
        }
    }
    let mut mr = vec![1.; (g.nr + 1) * g.nu];
    for i in 1..=g.nr {
        for j in 0..g.nu {
            let r = g.rf[i];
            let u = g.u[j];
            let left = g.id(i - 1, j);
            let (normal, tangent) = if i == g.nr {
                let db = if j == 0 {
                    (g.boundary[1] - g.boundary[0]) / g.du
                } else if j + 1 == g.nu {
                    (g.boundary[j] - g.boundary[j - 1]) / g.du
                } else {
                    (g.boundary[j + 1] - g.boundary[j - 1]) / (2. * g.du)
                };
                (
                    (g.boundary[j] - w[left]) / (r - g.r[i - 1]),
                    0.5 * (wu[left] + db),
                )
            } else {
                let right = g.id(i, j);
                (
                    (w[right] - w[left]) / (g.r[i] - g.r[i - 1]),
                    0.5 * (wu[right] + wu[left]),
                )
            };
            let gr = 1. / (r * r) - g.eta * u + normal;
            let gu = (-g.eta * r + tangent) * (1. - u * u).sqrt() / r;
            mr[i * g.nu + j] = response.get(gr.hypot(gu));
        }
    }
    let mut mu = vec![0.; g.nr * (g.nu + 1)];
    for i in 0..g.nr {
        for j in 1..g.nu {
            let uf = -1. + j as f64 * g.du;
            let r = g.r[i];
            let l = g.id(i, j - 1);
            let h = g.id(i, j);
            let gr = 1. / (r * r) - g.eta * uf + 0.5 * (wr[l] + wr[h]);
            let gu = (-g.eta * r + (w[h] - w[l]) / g.du) * (1. - uf * uf).sqrt() / r;
            mu[i * (g.nu + 1) + j] = response.get(gr.hypot(gu));
        }
    }
    let mut op = Operator {
        diag: vec![0.; n],
        radial: vec![0.; n],
        angular: vec![0.; n],
        rhs: vec![0.; n],
        nu: g.nu,
    };
    for i in 0..g.nr {
        for j in 0..g.nu {
            let k = g.id(i, j);
            let r = g.r[i];
            let u = g.u[j];
            let rlo = g.rf[i];
            let rhi = g.rf[i + 1];
            let ulo = -1. + j as f64 * g.du;
            let uhi = ulo + g.du;
            let mrl = mr[i * g.nu + j];
            let mrh = mr[(i + 1) * g.nu + j];
            let mul = mu[i * (g.nu + 1) + j];
            let muh = mu[i * (g.nu + 1) + j + 1];
            // Analytic harmonic background -1/r-eta*r*u prevents monopole discretization pollution.
            let frh = mrh * (1. - g.eta * u * rhi * rhi);
            let frl = mrl * (1. - g.eta * u * rlo * rlo);
            let fuh = (1. - uhi * uhi) * muh * (-g.eta * r);
            let ful = (1. - ulo * ulo) * mul * (-g.eta * r);
            op.rhs[k] = (frh - frl) * g.du + (fuh - ful) * g.dr[i];
            if i > 0 {
                op.diag[k] += rlo * rlo * mrl * g.du / (g.r[i] - g.r[i - 1]);
            }
            let ch = rhi * rhi * mrh * g.du
                / if i + 1 == g.nr {
                    rhi - r
                } else {
                    g.r[i + 1] - r
                };
            op.diag[k] += ch;
            if i + 1 == g.nr {
                op.rhs[k] += ch * g.boundary[j];
            } else {
                op.radial[k] = ch;
            }
            if j > 0 {
                op.diag[k] += (1. - ulo * ulo) * mul * g.dr[i] / g.du;
            }
            if j + 1 < g.nu {
                let cu = (1. - uhi * uhi) * muh * g.dr[i] / g.du;
                op.diag[k] += cu;
                op.angular[k] = cu;
            }
        }
    }
    op
}

fn dot(x: &[f64], y: &[f64]) -> f64 {
    x.iter().zip(y).map(|(a, b)| a * b).sum()
}
fn pcg(op: &Operator, x: &mut [f64], tol: f64) -> (usize, f64) {
    let n = x.len();
    let mut ax = vec![0.; n];
    op.apply(x, &mut ax);
    let mut r: Vec<f64> = op.rhs.iter().zip(&ax).map(|(b, a)| b - a).collect();
    let norm = dot(&op.rhs, &op.rhs).sqrt().max(1e-20);
    let mut z: Vec<f64> = r.iter().zip(&op.diag).map(|(r, d)| r / d).collect();
    let mut p = z.clone();
    let mut rz = dot(&r, &z);
    if dot(&r, &r).sqrt() / norm < tol {
        return (0, dot(&r, &r).sqrt() / norm);
    }
    for it in 0..10000 {
        op.apply(&p, &mut ax);
        let pap = dot(&p, &ax);
        if pap <= 0. || !pap.is_finite() {
            panic!("PCG operator not positive: {pap}");
        }
        let alpha = rz / pap;
        for k in 0..n {
            x[k] += alpha * p[k];
            r[k] -= alpha * ax[k];
        }
        let error = dot(&r, &r).sqrt() / norm;
        if error < tol {
            return (it + 1, error);
        }
        for k in 0..n {
            z[k] = r[k] / op.diag[k];
        }
        let next = dot(&r, &z);
        let beta = next / rz;
        for k in 0..n {
            p[k] = z[k] + beta * p[k];
        }
        rz = next;
    }
    (10000, dot(&r, &r).sqrt() / norm)
}

fn quadrupole(g: &Grid, w: &[f64], low: f64, high: f64) -> f64 {
    let mut p: Vec<f64> = g.u.iter().map(|u| 0.5 * (3. * u * u - 1.)).collect();
    let mean = p.iter().sum::<f64>() / g.nu as f64;
    for x in &mut p {
        *x -= mean;
    }
    let norm = dot(&p, &p);
    let (mut numer, mut denom) = (0., 0.);
    for i in 0..g.nr {
        let r = g.r[i];
        if r < low || r > high {
            continue;
        }
        let c = dot(&w[i * g.nu..(i + 1) * g.nu], &p) / norm;
        numer += r * r * c;
        denom += r.powi(4);
    }
    numer / denom
}

fn response_csv(path: &str) {
    let (nodes, weights) = gauss_legendre(96);
    let mut f = fs::File::create(path).unwrap();
    writeln!(f, "x,standard,linear,tanh").unwrap();
    for i in 0..401 {
        let x = 10f64.powf(-5. + 10. * i as f64 / 400.);
        writeln!(
            f,
            "{x:.16e},{:.16e},{:.16e},{:.16e}",
            exact_mu(Model::Standard, x, &nodes, &weights),
            exact_mu(Model::Linear, x, &nodes, &weights),
            exact_mu(Model::Tanh, x, &nodes, &weights)
        )
        .unwrap();
    }
}

fn table_response_csv(input: &str, output: &str) {
    let csv = fs::read_to_string(input).expect("Cannot read response table");
    let response = Response::from_table(&csv, 0.).expect("Invalid response table");
    let mut file = fs::File::create(output).unwrap();
    writeln!(file, "x,mu,log_slope").unwrap();
    for i in 0..1201 {
        let x = 10f64.powf(-12. + 24. * i as f64 / 1200.);
        writeln!(
            file,
            "{x:.16e},{:.16e},{:.16e}",
            response.get(x),
            response.log_slope(x)
        )
        .unwrap();
    }
}

fn main() {
    let args: Vec<String> = env::args().collect();
    if args.len() == 4 && args[1] == "table-response" {
        table_response_csv(&args[2], &args[3]);
        return;
    }
    if args.len() > 1 && args[1] == "response" {
        response_csv(&args[2]);
        return;
    }
    if args.len() < 8 {
        eprintln!(
            "Usage: aqual MODEL Nr Nu rmin rmax g_ext_SI output_prefix [mu_floor] [table_csv frozen_scale_SI]"
        );
        std::process::exit(2);
    }
    let model = Model::parse(&args[1]);
    let nr = args[2].parse().unwrap();
    let nu = args[3].parse().unwrap();
    let rmin = args[4].parse().unwrap();
    let rmax = args[5].parse().unwrap();
    let ext: f64 = args[6].parse().unwrap();
    let prefix = &args[7];
    let floor = args.get(8).map(|x| x.parse().unwrap()).unwrap_or(1e-9);
    let (response, scale) = if matches!(model, Model::Table) {
        if args.len() != 11 {
            panic!(
                "Table mode requires mu_floor, table_csv and frozen_scale_SI after output_prefix"
            );
        }
        let csv = fs::read_to_string(&args[9]).expect("Cannot read response table");
        let response = Response::from_table(&csv, floor).expect("Invalid response table");
        let scale: f64 = args[10].parse().expect("Invalid frozen acceleration scale");
        assert!(scale.is_finite() && scale > 0.);
        (response, scale)
    } else {
        (Response::new(model, floor), model.scale())
    };
    let eta = ext / scale;
    let grid = Grid::new(nr, nu, rmin, rmax, eta, &response);
    let n = nr * nu;
    let mut w = vec![0.; n];
    let mut iterations = 0;
    let mut cg_total = 0;
    let start = Instant::now();
    let mut nonlinear_error = f64::INFINITY;
    let mut converged = false;
    for it in 0..180 {
        let op = assemble(&grid, &response, &w);
        let mut solved = w.clone();
        let (cg, linear_error) = pcg(&op, &mut solved, 2e-11);
        cg_total += cg;
        if cg >= 10000 {
            eprintln!("Linear solve did not converge: {linear_error}");
            break;
        }
        let mut delta: f64 = 0.;
        for k in 0..n {
            delta = delta.max((solved[k] - w[k]).abs());
            w[k] = 0.6 * solved[k] + 0.4 * w[k];
        }
        let current = assemble(&grid, &response, &w);
        let mut ax = vec![0.; n];
        current.apply(&w, &mut ax);
        for k in 0..n {
            ax[k] -= current.rhs[k];
        }
        nonlinear_error = dot(&ax, &ax).sqrt() / dot(&current.rhs, &current.rhs).sqrt().max(1e-20);
        iterations = it + 1;
        if it % 10 == 0 {
            eprintln!("iter={it} delta={delta:.3e} residual={nonlinear_error:.3e} cg={cg}");
        }
        if delta < 2e-10 && nonlinear_error < 2e-9 {
            converged = true;
            break;
        }
    }
    let radius = (6.67430e-11 * 1.98847e30 / scale).sqrt();
    let windows = [(0.02, 0.08), (0.04, 0.16), (0.08, 0.32)];
    let coeffs: Vec<f64> = windows
        .iter()
        .map(|&(l, h)| quadrupole(&grid, &w, l, h))
        .collect();
    let q2: Vec<f64> = coeffs.iter().map(|a| -3. * a * scale / radius).collect();
    let report = format!(
        concat!(
            "{{\n  \"model\": \"{:?}\",\n  \"nr\": {}, \"nu\": {},\n",
            "  \"rmin\": {}, \"rmax\": {}, \"external_acceleration_m_s2\": {},\n",
            "  \"model_scale_m_s2\": {}, \"eta\": {}, \"mu_floor\": {},\n",
            "  \"converged\": {}, \"nonlinear_residual\": {:.10e}, \"picard_iterations\": {}, \"cg_iterations\": {},\n",
            "  \"quadrupole_coefficients\": {:?},\n  \"Q2_s_minus2\": {:?},\n",
            "  \"fit_windows\": [[0.02,0.08],[0.04,0.16],[0.08,0.32]],\n",
            "  \"seconds\": {:.3}\n}}\n"
        ),
        model,
        nr,
        nu,
        rmin,
        rmax,
        ext,
        scale,
        eta,
        floor,
        converged,
        nonlinear_error,
        iterations,
        cg_total,
        coeffs,
        q2,
        start.elapsed().as_secs_f64()
    );
    fs::write(format!("{prefix}.json"), &report).unwrap();
    println!("{report}");
    let mut file = fs::File::create(format!("{prefix}.csv")).unwrap();
    writeln!(file, "r,u,correction").unwrap();
    for i in 0..nr {
        for j in 0..nu {
            writeln!(
                file,
                "{:.12e},{:.12e},{:.16e}",
                grid.r[i],
                grid.u[j],
                w[grid.id(i, j)]
            )
            .unwrap();
        }
    }
    if !converged {
        std::process::exit(3);
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    fn example_table(values: &[f64]) -> String {
        let mut text = String::from("x,mu\n");
        for (i, value) in values.iter().enumerate() {
            text.push_str(&format!("{:.16e},{value:.16e}\n", 10f64.powi(i as i32 - 2)));
        }
        text
    }
    #[test]
    fn response_table_preserves_super_newtonian_values_and_low_slope() {
        let table = Response::from_table(&example_table(&[0.02, 0.2, 1.2, 1.]), 1e-9).unwrap();
        assert!((table.get(1.) - 1.2).abs() < 1e-14);
        assert!((table.get(1e-3) - 0.002).abs() < 1e-14);
        assert_eq!(table.get(0.), 1e-9);
        assert_eq!(table.get(100.), 1.);
        assert!(table.get(3.) > 1.);
    }
    #[test]
    fn response_table_log_interpolation_and_slope_are_consistent() {
        let table = Response::from_table(&example_table(&[0.02, 0.2, 1.2, 1.]), 0.).unwrap();
        assert!((table.get(0.1f64.sqrt()) - 0.7).abs() < 1e-14);
        let analytic = (1.0 / 10f64.ln()) / 0.7;
        assert!((table.log_slope(0.1f64.sqrt()) - analytic).abs() < 2e-5);
    }
    #[test]
    fn response_table_rejects_invalid_values_and_meshes() {
        for csv in [
            "x,mu\n0,1\n1,1\n10,1\n",
            "x,mu\n1,1\n2,1\n10,1\n",
            "x,mu\n1,NaN\n10,1\n100,1\n",
            "x,mu\n1,1\n10,-1\n100,1\n",
            "wrong\n1,1\n10,1\n100,1\n",
        ] {
            assert!(Response::from_table(csv, 0.).is_err());
        }
    }
    #[test]
    fn response_table_rejects_noninvertible_or_nonnewtonian_maps() {
        assert!(Response::from_table(&example_table(&[0.2, 0.001, 1.]), 0.).is_err());
        assert!(Response::from_table(&example_table(&[0.02, 0.2, 0.9]), 0.).is_err());
        // Node-wise invertibility alone does not protect an interpolated interval.
        assert!(Response::from_table(&example_table(&[0.02, 0.005, 1.]), 0.).is_err());
    }
    #[test]
    fn quadrature_moments() {
        let (u, w) = gauss_legendre(96);
        for power in 0..12 {
            let m = u
                .iter()
                .zip(&w)
                .map(|(u, w)| w * u.powi(power))
                .sum::<f64>();
            assert!((m - 1. / (power + 1) as f64).abs() < 1e-13);
        }
    }
    #[test]
    fn harmonic_background_cancels() {
        let response = Response::new(Model::Newton, 1e-9);
        let grid = Grid::new(24, 20, 0.003, 20., 1.7, &response);
        let op = assemble(&grid, &response, &vec![0.; 24 * 20]);
        assert!(op.rhs.iter().all(|x| x.abs() < 1e-12));
    }
    #[test]
    fn finite_volume_matrix_symmetric() {
        let response = Response::new(Model::Tanh, 1e-9);
        let grid = Grid::new(18, 16, 0.003, 20., 1.5, &response);
        let n = 18 * 16;
        let op = assemble(&grid, &response, &vec![0.; n]);
        let a: Vec<f64> = (0..n).map(|i| (i as f64 * 0.37).sin()).collect();
        let b: Vec<f64> = (0..n).map(|i| (i as f64 * 0.17).cos()).collect();
        let mut aa = vec![0.; n];
        let mut ab = vec![0.; n];
        op.apply(&a, &mut aa);
        op.apply(&b, &mut ab);
        assert!((dot(&a, &ab) - dot(&b, &aa)).abs() < 1e-10);
        assert!(dot(&a, &aa) > 0.);
    }
    #[test]
    fn manufactured_regular_quadrupole() {
        let response = Response::new(Model::Newton, 1e-9);
        let mut grid = Grid::new(96, 72, 0.003, 2., 1.5, &response);
        let amplitude = 0.1;
        for (j, u) in grid.u.iter().enumerate() {
            grid.boundary[j] = amplitude * 4. * 0.5 * (3. * u * u - 1.);
        }
        let mut w = vec![0.; 96 * 72];
        let op = assemble(&grid, &response, &w);
        let (_, residual) = pcg(&op, &mut w, 1e-12);
        let recovered = quadrupole(&grid, &w, 0.04, 0.16);
        assert!(residual < 1e-11);
        assert!((recovered / amplitude - 1.).abs() < 0.01, "{recovered}");
    }
}
