physica 0.3.13

//! US 12,417,870 ("driving movement of a target object by a gravitational field generated by a
//! varying electromagnetic field", Zhang, Xiao & Li, 2025) rebuilt from its embodiment and run as
//! a simulation: a 0.35 g polyethylene ball on a cotton thread in a 10 cm vacuum tank, between two
//! open high-voltage coils, with every force and torque on the ball computed from the crates and
//! integrated in time. Then the same authors' weighing rig.
//!
//! The rig is a calibration case, not part of any crate's API: the physics it uses is general and
//! lives in the libraries (`gravitics::electrodes` for the field of any rings and wires,
//! `gravitics::electrogravity` for the forces and torques and Zhang's constants,
//! `plasmoid::discharge` for the gas, `magnetics` for Earth's field and the coil's), and this file
//! only puts the patent's geometry around them. `web/gravity-patent.html` runs the same model in
//! the browser and checks itself against [`reference_values_for_the_web_model`].
//!
//! The geometry is the patent's embodiment, with Zhang's book for the coils: two 12-turn coils of
//! 1.5 mm copper, 3.5 cm in radius, 6–18 cm above and below the tank's centre; the upper one fed by
//! the supply, the lower one grounded; their inner ends cut, bare and staggered beside the tank
//! (the patent's FIG. 3); their outer ends led to the supply and to ground.
//!
//! On the ball, every 1 ms step:
//! - gravity and buoyancy, and the pull of the field's own energy (general relativity,
//!   `GU/(c²d²)`, about 10⁻²⁷ m/s², below what a double can add to 9.8);
//! - the polarisation pull toward strong field, the Coulomb force on the ball's charge, the force
//!   on its charge dipole in the field's gradient, and its image in the tank wall;
//! - gas drag at any Knudsen number; breakdown, which clamps the voltage across the tank and parks
//!   the rest on its walls; the floating potential a glow drains the ball to; and the ion drag of
//!   the glow's current;
//! - the coil's charging and leakage current plus Earth's field (WMM2025 at Zhang's lab), acting
//!   on the moving charge and on the ball's diamagnetism;
//! - about the thread: its torsion and damping, the gas's spin drag, the field turning a ball that
//!   is not quite round and its charge dipole, and, when chosen, Zhang's vortex field with one of
//!   his constants.
//!
//! What the patent leaves out (the ball's charge, its out-of-roundness, the thread's stiffness and
//! the pressure) is in `Settings`, with defaults typical of a handled polymer ball on cotton at a
//! rotary pump's 1 Pa. Left out: the dielectrics (silicone insulation, glass), the thread's own
//! charge, the forces of the 30 kHz ripple itself (a 1% ripple changes the field-squared forces by
//! about 2%), radiometric forces, and vibration.
//!
//! ```text
//! cargo run -p physica --example gravity_patent
//! ```

use core::f64::consts::PI;

use gravitics::electrodes::{Element, Ring, charges, field, field_at, polyline, solve, stack};
use gravitics::electrogravity::{
    ZHANG_KMQ, alignment_torque, image_force, sphere_force, spheroid_polarisability, torsion_constant, vortex_constant, vortex_spin,
};
use gravitics::propulsion::ion_wind_thrust;
use magnetics::magnets::{solenoid_bz, susceptibility_force};
use metrology::{C, EPS_0, G, G0};
use plasmoid::discharge::{
    AIR, AIR_DIAMETER, AIR_ION_MOBILITY, AIR_MOLAR_MASS, AIR_VISCOSITY, N2_ION_MASS, floating_potential, gas_density,
    ion_mobility, mean_free_path, mean_speed, paschen, peek_onset, sphere_drag, sphere_spin_drag,
};

/// The copper's radius (m).
const WIRE: f64 = 0.75e-3;
/// The tank's radius (m); it lies along x.
const TANK_RADIUS: f64 = 0.05;
/// Where the thread is tied, at the top of the tank (m); the ball hangs at the origin.
const SUSPENSION: [f64; 3] = [0.0, 0.0, 0.05];
/// The coils' mean radius (m), 12 turns each, 1 cm apart.
const COIL_RADIUS: f64 = 0.035;
/// The current per volt the patent reports, 10 µA at 50 kV (A/V).
const LEAKAGE: f64 = 10e-6 / 50e3;
/// The frequency of the supply's ripple (Hz), from Zhang's ZVS driver.
const RIPPLE_FREQUENCY: f64 = 30e3;
/// Zhang's lab in Lujiang County, Anhui (degrees north and east).
const LUJIANG: (f64, f64) = (31.25, 117.29);
/// How long a glow takes to drain the ball to its floating charge (s): the ion flux of a weak glow
/// (10¹⁴ m⁻³ at the Bohm speed) onto a 9 mm ball.
const GLOW_CHARGING: f64 = 1e-3;
/// The fixed time step (s); [`Rig::run`] takes whole steps, so every run of the same settings
/// gives the same trajectory.
const DT: f64 = 1e-3;

// The leads as polylines (m): the coils' inner ends cut and staggered beside the tank, with 1.1 cm
// of bare copper at each tip; the outer ends to the supply and to ground.
const LEAD_TOP_INNER: [[f64; 3]; 5] = [[0.0, -0.035, 0.06], [-0.005, -0.054, 0.056], [-0.016, -0.061, 0.04], [-0.016, -0.061, -0.016], [-0.016, -0.061, -0.027]];
const LEAD_BOTTOM_INNER: [[f64; 3]; 5] = [[0.0, -0.035, -0.06], [0.005, -0.054, -0.056], [0.016, -0.061, -0.04], [0.016, -0.061, 0.016], [0.016, -0.061, 0.027]];
const LEAD_TOP_OUTER: [[f64; 3]; 6] = [[0.0, -0.035, 0.18], [0.0, -0.05, 0.205], [-0.09, -0.03, 0.215], [-0.21, 0.0, 0.15], [-0.25, 0.0, -0.06], [-0.26, 0.0, -0.176]];
const LEAD_BOTTOM_OUTER: [[f64; 3]; 4] = [[0.0, -0.035, -0.18], [0.0, -0.05, -0.205], [-0.06, -0.06, -0.225], [-0.12, -0.06, -0.234]];

/// The rig's conductors, and whether each is on the supply side (the upper coil and its leads).
fn geometry() -> (Vec<Element>, Vec<bool>) {
    let mut parts: Vec<(Element, bool)> = Vec::new();
    for (z0, step, live) in [(0.065, 0.01, true), (-0.065, -0.01, false)] {
        parts.extend(stack(COIL_RADIUS, z0, step, 12).into_iter().map(|r| (Element::Ring(r), live)));
    }
    for (points, step, live) in [
        (&LEAD_TOP_INNER[..], 0.01, true),
        (&LEAD_BOTTOM_INNER[..], 0.01, false),
        (&LEAD_TOP_OUTER[..], 0.02, true),
        (&LEAD_BOTTOM_OUTER[..], 0.02, false),
    ] {
        parts.extend(polyline(points, step).into_iter().map(|s| (Element::Segment(s), live)));
    }
    parts.into_iter().unzip()
}

/// How the supply drives the coils.
#[derive(Clone, Copy, Debug, PartialEq)]
enum Drive {
    /// The patent's: the upper coil at the supply's voltage, the lower one grounded.
    Single,
    /// The upper coil at +V/2 and the lower at −V/2, the test that removes the live side.
    Symmetric,
}

/// The suspended ball.
#[derive(Clone, Copy, Debug, PartialEq)]
struct Ball {
    /// Radius (m).
    pub radius: f64,
    /// Density (kg/m³).
    pub density: f64,
    /// Relative permittivity.
    pub eps_r: f64,
    /// Volume magnetic susceptibility (SI).
    pub susceptibility: f64,
    /// Long axis over short axes: 1 for a perfect sphere.
    pub aspect: f64,
    /// Net charge (C).
    pub charge: f64,
    /// The dipole moment of its charge patches (C m), lying along its long axis.
    pub dipole: f64,
}

impl Ball {
    /// The patent's ball: 0.35 g of polyethylene (950 kg/m³, ε_r = 2.3, χ ≈ −1×10⁻⁵), 1% out of
    /// round, carrying 0.2 nC and a patch dipole of 0.1 nC across its radius, as a handled polymer
    /// ball typically does.
    fn polyethylene() -> Ball {
        let radius = (3.0 * 0.35e-3 / 950.0 / (4.0 * PI)).cbrt();
        Ball { radius, density: 950.0, eps_r: 2.3, susceptibility: -9.7e-6, aspect: 1.01, charge: 0.2e-9, dipole: 0.1e-9 * radius }
    }

    /// The same ball in PTFE (2200 kg/m³, ε_r = 2.1): the mass-independence test.
    fn ptfe() -> Ball {
        Ball { density: 2200.0, eps_r: 2.1, susceptibility: -1.05e-5, ..Ball::polyethylene() }
    }

    fn volume(&self) -> f64 {
        4.0 / 3.0 * PI * self.radius.powi(3)
    }

    fn mass(&self) -> f64 {
        self.density * self.volume()
    }

    /// Moment of inertia about the thread (kg m²).
    fn inertia(&self) -> f64 {
        0.4 * self.mass() * self.radius * self.radius
    }
}

/// Everything the simulation needs that the patent does not fix.
#[derive(Clone, Debug, PartialEq)]
struct Settings {
    /// The supply's set voltage (V).
    pub volts: f64,
    pub drive: Drive,
    /// Gas pressure in the tank (Pa) and temperature (K).
    pub pressure: f64,
    pub temperature: f64,
    pub ball: Ball,
    /// The thread's torsion constant (N m/rad) and quality factor.
    pub torsion: f64,
    pub thread_q: f64,
    /// The supply's rise and decay time constants (s).
    pub rise: f64,
    pub fall: f64,
    /// The supply's ripple, as a fraction of its voltage.
    pub ripple: f64,
    /// The electron temperature of any glow (eV).
    pub electron_temperature: f64,
    /// Earth's field in the rig's frame (T): x along the tank (taken as north), y west, z up.
    pub earth_field: [f64; 3],
    /// Zhang's vortex constant `f` (kg s/C), when his model is switched on.
    pub vortex: Option<f64>,
}

impl Default for Settings {
    /// The patent at 50 kV, at 1 Pa, with a thread that lets the ball twist back and forth every
    /// 30 s (Q = 20), a supply that rises in 50 ms and decays in 1 s with 1% ripple.
    fn default() -> Settings {
        let ball = Ball::polyethylene();
        Settings {
            volts: 50e3,
            drive: Drive::Single,
            pressure: 1.0,
            temperature: 293.15,
            ball,
            torsion: torsion_constant(ball.inertia(), 30.0),
            thread_q: 20.0,
            rise: 0.05,
            fall: 1.0,
            ripple: 0.01,
            electron_temperature: 3.0,
            earth_field: earth_field(),
            vortex: None,
        }
    }
}

/// Earth's field at Zhang's lab (T), from WMM2025 at 2025.0, in the rig's frame (x north, y west,
/// z up).
fn earth_field() -> [f64; 3] {
    let e = magnetics::main_field::wmm2025().elements(2025.0, LUJIANG.0.to_radians(), LUJIANG.1.to_radians(), 0.0);
    [e.x * 1e-9, -e.y * 1e-9, -e.z * 1e-9]
}

/// The rig, solved for its drive.
#[derive(Clone, Debug)]
struct Rig {
    pub settings: Settings,
    elements: Vec<Element>,
    /// Charge on each element per volt of the supply (C/V).
    unit: Vec<f64>,
    /// The field at the tank's centre per volt (V/m per V).
    centre: f64,
    /// Charge on the supply side per volt: its capacitance (F).
    capacitance: f64,
    /// Field energy per volt squared (J/V²).
    energy: f64,
}

impl Rig {
    fn new(settings: Settings) -> Rig {
        let (elements, live) = geometry();
        let (top, bottom) = match settings.drive {
            Drive::Single => (1.0, 0.0),
            Drive::Symmetric => (0.5, -0.5),
        };
        let volts: Vec<f64> = live.iter().map(|&l| if l { top } else { bottom }).collect();
        let unit = solve(&elements, WIRE, &volts);
        let capacitance = unit.iter().zip(&live).filter(|(_, l)| **l).map(|(q, _)| q).sum();
        let energy = 0.5 * unit.iter().zip(&volts).map(|(q, v)| q * v).sum::<f64>();
        let [x, y, z] = field_at(&elements, &unit, [0.0; 3]);
        Rig { settings, elements, unit, centre: (x * x + y * y + z * z).sqrt(), capacitance, energy }
    }

    /// The rig at rest with the supply off, the ball carrying its charge and patches.
    fn at_rest(&self) -> State {
        State { charge: self.settings.ball.charge, dipole: self.settings.ball.dipole, ..State::default() }
    }

    /// The supply side's capacitance (F): its charge per volt.
    fn capacitance(&self) -> f64 {
        self.capacitance
    }

    /// The surface field (V/m) per volt of the supply at the bare copper tips, the stronger of the
    /// two: the solved line charge of each tip's end segment over `2πε₀a`.
    fn tip_field(&self) -> f64 {
        let tips = [LEAD_TOP_INNER[4], LEAD_BOTTOM_INNER[4]];
        let near = |p: [f64; 3], q: [f64; 3]| (0..3).map(|k| (p[k] - q[k]).powi(2)).sum::<f64>().sqrt() < 1e-9;
        self.elements
            .iter()
            .zip(&self.unit)
            .filter_map(|(e, q)| match e {
                Element::Segment(s) if tips.iter().any(|&t| near(s.b, t)) => {
                    let l = (0..3).map(|k| (s.b[k] - s.a[k]).powi(2)).sum::<f64>().sqrt();
                    Some((q / l).abs() / (2.0 * PI * EPS_0 * WIRE))
                }
                _ => None,
            })
            .fold(0.0, f64::max)
    }

    /// The field's energy per volt squared of the supply (J/V²).
    fn energy(&self) -> f64 {
        self.energy
    }

    /// The field (V/m) at `p` (m) per volt across the tank's gas.
    fn unit_field(&self, p: [f64; 3]) -> [f64; 3] {
        field_at(&self.elements, &self.unit, p)
    }

    /// The voltage across the tank's gas at which it breaks down, in volts of the supply: the
    /// Paschen voltage over the tank's 10 cm, over the drop the coils put across it. `None` when
    /// it cannot break down.
    fn breakdown(&self) -> Option<f64> {
        let span = 2.0 * TANK_RADIUS;
        paschen(AIR, 0.01, self.settings.pressure * span).map(|v| v / (self.centre * span))
    }

    /// Advance one step of [`DT`].
    fn step(&self, s: &mut State) -> Report {
        let set = &self.settings;
        let ball = &set.ball;
        let (m, a, vol) = (ball.mass(), ball.radius, ball.volume());

        // The supply.
        let target = if s.on { set.volts } else { 0.0 };
        let tau = if s.on { set.rise } else { set.fall };
        let before = s.supply;
        s.supply += (target - before) * (1.0 - (-DT / tau).exp());
        let ramp = (s.supply - before) / DT;

        // The gas: past breakdown it carries the excess to the tank's walls, clamping the voltage
        // across itself; the ripple alone may keep a glow lit.
        let wall_before = s.wall;
        let limit = self.breakdown();
        let clamping = match limit {
            // A part in 10⁹ keeps the supply's last round-off creep from counting as a discharge.
            Some(l) if (s.supply - s.wall).abs() > l * (1.0 + 1e-9) => {
                s.wall = s.supply - l * (s.supply - s.wall).signum();
                true
            }
            _ => false,
        };
        let sustained = limit.is_some_and(|l| set.ripple * s.supply.abs() > l);
        s.glow = clamping || sustained;
        let across = s.supply - s.wall;
        let mut glow_current = self.capacitance * (s.wall - wall_before).abs() / DT;
        if sustained {
            glow_current += 2.0 * PI * RIPPLE_FREQUENCY * self.capacitance * set.ripple * s.supply.abs();
        }
        let glow_current = glow_current.min(LEAKAGE * set.volts);

        // The field at the ball and the gradient of its square.
        let field = |p: [f64; 3]| self.unit_field(p).map(|e| e * across);
        let r = s.position;
        let e = field(r);
        let h = 1e-4;
        let shifted = |k: usize, d: f64| {
            let mut p = r;
            p[k] += d;
            field(p)
        };
        let grad_e2: [f64; 3] = core::array::from_fn(|k| (norm2(shifted(k, h)) - norm2(shifted(k, -h))) / (2.0 * h));

        // The ball's charge and its patches: kept in vacuum, drained by a glow to the floating
        // potential and to an even cover.
        if s.glow {
            let floating = 4.0 * PI * EPS_0 * a * floating_potential(set.electron_temperature, N2_ION_MASS);
            let decay = (-DT / GLOW_CHARGING).exp();
            s.charge = floating + (s.charge - floating) * decay;
            s.dipole *= decay;
        }
        let q = s.charge;

        // The gas.
        let rho_gas = gas_density(set.pressure, set.temperature, AIR_MOLAR_MASS);
        let knudsen = mean_free_path(set.pressure, set.temperature, AIR_DIAMETER) / a;

        // Forces.
        let weight = [0.0, 0.0, -m * G0];
        let buoyancy = [0.0, 0.0, rho_gas * vol * G0];
        let polarisation = grad_e2.map(|g| sphere_force(a, ball.eps_r, g));
        let coulomb = e.map(|x| q * x);
        let axis = [s.twist.cos(), s.twist.sin(), 0.0];
        let dipole_force = {
            let plus = field([r[0] + h * axis[0], r[1] + h * axis[1], r[2]]);
            let minus = field([r[0] - h * axis[0], r[1] - h * axis[1], r[2]]);
            core::array::from_fn(|k| s.dipole * (plus[k] - minus[k]) / (2.0 * h))
        };
        let image = {
            let off = r[1].hypot(r[2]);
            if off == 0.0 {
                [0.0; 3]
            } else {
                let f = image_force(q, TANK_RADIUS - off, 6.0) - image_force(q, TANK_RADIUS + off, 6.0);
                [0.0, f * r[1] / off, f * r[2] / off]
            }
        };
        let ion_drag = {
            let size = norm2(e).sqrt();
            let push = ion_wind_thrust(glow_current, 2.0 * TANK_RADIUS, ion_mobility(AIR_ION_MOBILITY, set.pressure, set.temperature)) * (a / TANK_RADIUS).powi(2);
            if size > 0.0 { e.map(|x| push * x / size) } else { [0.0; 3] }
        };
        let gas_drag = s.velocity.map(|v| -sphere_drag(a, AIR_VISCOSITY, knudsen) * v);
        let current = LEAKAGE * s.supply + self.capacitance * ramp;
        let coil = |z: f64| solenoid_bz(0.12 - z, COIL_RADIUS, 0.12, 100.0, current);
        let b = [set.earth_field[0], set.earth_field[1], set.earth_field[2] + coil(r[2])];
        let lorentz = cross(s.velocity, b).map(|x| q * x);
        let magnetisation = [0.0, 0.0, susceptibility_force(vol, ball.susceptibility, b[2] * (coil(r[2] + h) - coil(r[2] - h)) / (2.0 * h))];
        let field_gravity = match set.drive {
            Drive::Single => {
                let to_coil = [-r[0], -r[1], 0.12 - r[2]];
                let d = norm2(to_coil).sqrt();
                let pull = G * m * self.energy * s.supply * s.supply / (C * C * d * d);
                to_coil.map(|x| pull * x / d)
            }
            Drive::Symmetric => [0.0; 3],
        };
        let total: [f64; 3] = core::array::from_fn(|k| {
            weight[k] + buoyancy[k] + polarisation[k] + coulomb[k] + dipole_force[k] + image[k] + ion_drag[k] + gas_drag[k] + lorentz[k] + magnetisation[k] + field_gravity[k]
        });

        // The thread: integrate, then hold the ball on the sphere the thread allows.
        let mut v: [f64; 3] = core::array::from_fn(|k| s.velocity[k] + total[k] / m * DT);
        let moved: [f64; 3] = core::array::from_fn(|k| r[k] + v[k] * DT - SUSPENSION[k]);
        let length = SUSPENSION[2];
        let n = moved.map(|x| x / norm2(moved).sqrt());
        s.position = core::array::from_fn(|k| SUSPENSION[k] + length * n[k]);
        let radial = v[0] * n[0] + v[1] * n[1] + v[2] * n[2];
        v = core::array::from_fn(|k| v[k] - radial * n[k]);
        s.velocity = v;

        // Torques about the thread.
        let inertia = ball.inertia();
        let horizontal = e[0].hypot(e[1]);
        let azimuth = e[1].atan2(e[0]);
        let (long, short) = spheroid_polarisability(ball.eps_r, ball.aspect);
        let thread_torque = -set.torsion * s.twist;
        let damping_torque = -(inertia * (set.torsion / inertia).sqrt() / set.thread_q + sphere_spin_drag(a, AIR_VISCOSITY, rho_gas, mean_speed(set.temperature, AIR_MOLAR_MASS))) * s.spin;
        let alignment = -alignment_torque(vol, long - short, horizontal, s.twist - azimuth);
        let dipole_torque = s.dipole * horizontal * (azimuth - s.twist).sin();
        let vortex_torque = set.vortex.map_or(0.0, |f| inertia * vortex_spin(b[2], f));
        s.spin += (thread_torque + damping_torque + alignment + dipole_torque + vortex_torque) / inertia * DT;
        s.twist += s.spin * DT;
        s.t += DT;

        Report {
            field: e,
            across,
            glow: s.glow,
            glow_current,
            current,
            b,
            weight,
            buoyancy,
            polarisation,
            coulomb,
            dipole_force,
            image,
            ion_drag,
            gas_drag,
            lorentz,
            magnetisation,
            field_gravity,
            thread_torque,
            damping_torque,
            alignment,
            dipole_torque,
            vortex_torque,
        }
    }

    /// Advance `duration` (s) in whole steps of [`DT`], returning the last step's report.
    fn run(&self, s: &mut State, duration: f64) -> Report {
        let mut report = Report::default();
        for _ in 0..(duration / DT).round() as usize {
            report = self.step(s);
        }
        report
    }
}

/// The rig's state.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
struct State {
    /// Time (s), and whether the supply is switched on.
    pub t: f64,
    pub on: bool,
    /// The supply's voltage now (V), and the part of it the tank's walls hold off the gas (V).
    pub supply: f64,
    pub wall: f64,
    /// The ball's centre (m) and velocity (m/s); at rest it hangs at the origin.
    pub position: [f64; 3],
    pub velocity: [f64; 3],
    /// Its turn about the thread from the thread's rest (rad), and its rate (rad/s).
    pub twist: f64,
    pub spin: f64,
    /// Its charge (C) and patch dipole (C m), and whether a glow is lit.
    pub charge: f64,
    pub dipole: f64,
    pub glow: bool,
}

/// Every force (N) and torque (N m) on the ball in one step, and what set them.
#[derive(Clone, Copy, Debug, Default, PartialEq)]
struct Report {
    /// The field at the ball (V/m), and the voltage across the gas (V).
    pub field: [f64; 3],
    pub across: f64,
    pub glow: bool,
    /// The glow's current and the coil's current (A), and the magnetic field at the ball (T).
    pub glow_current: f64,
    pub current: f64,
    pub b: [f64; 3],
    pub weight: [f64; 3],
    pub buoyancy: [f64; 3],
    pub polarisation: [f64; 3],
    pub coulomb: [f64; 3],
    pub dipole_force: [f64; 3],
    pub image: [f64; 3],
    pub ion_drag: [f64; 3],
    pub gas_drag: [f64; 3],
    pub lorentz: [f64; 3],
    pub magnetisation: [f64; 3],
    pub field_gravity: [f64; 3],
    pub thread_torque: f64,
    pub damping_torque: f64,
    pub alignment: f64,
    pub dipole_torque: f64,
    pub vortex_torque: f64,
}

fn norm2(v: [f64; 3]) -> f64 {
    v[0] * v[0] + v[1] * v[1] + v[2] * v[2]
}

fn cross(a: [f64; 3], b: [f64; 3]) -> [f64; 3] {
    [a[1] * b[2] - a[2] * b[1], a[2] * b[0] - a[0] * b[2], a[0] * b[1] - a[1] * b[0]]
}

fn main() {
    let rig = Rig::new(Settings::default());
    let e = rig.unit_field([0.0; 3]).map(|x| x * 50e3);
    println!("US 12,417,870 at 50 kV: field at the ball ({:.3e}, {:.3e}, {:.3e}) V/m, tilted {:.1}° by the leads", e[0], e[1], e[2], e[0].hypot(e[1]).atan2(-e[2]).to_degrees());
    println!("  supply side {:.2} pF, {:.2} mJ stored", rig.capacitance() * 1e12, rig.energy() * 2.5e9 * 1e3);
    let tip = rig.tip_field();
    println!("  the bare copper tips reach {:.1} MV/m: corona in air above {:.1} kV (Peek)", tip * 50e3 / 1e6, peek_onset(WIRE, 1.0) / tip / 1e3);
    let b = rig.settings.earth_field;
    println!("  Earth's field at Lujiang {LUJIANG:?}: ({:.1}, {:.1}, {:.1}) µT", b[0] * 1e6, b[1] * 1e6, b[2] * 1e6);

    // Switch on for 4 s, then off for 4 s.
    let mut s = State { on: true, ..rig.at_rest() };
    println!("switched on, then off at 4 s ({} Pa, ball {:.1} nC, patches {:.1} nC across its radius, 1% out of round):", rig.settings.pressure, s.charge * 1e9, s.dipole / rig.settings.ball.radius * 1e9);
    for step in 1..=8000 {
        if step == 4001 {
            s.on = false;
        }
        let r = rig.step(&mut s);
        if step % 1000 == 0 {
            println!("  {:.0} s: supply {:5.1} kV, turned {:+6.1}°, swung {:5.1} µm, glow {}", s.t, s.supply / 1e3, s.twist.to_degrees(), s.position[0].hypot(s.position[1]) * 1e6, s.glow);
        }
        if step == 3000 {
            print_forces(&r, rig.settings.ball.mass());
        }
    }

    // The tests that would settle it, and the gas.
    let neutral = Ball { charge: 0.0, dipole: 0.0, aspect: 1.0, ..Ball::polyethylene() };
    for (name, settings) in [
        ("a neutral, round ball", Settings { ball: neutral, ..Settings::default() }),
        ("the PTFE ball", Settings { ball: Ball::ptfe(), ..Settings::default() }),
        ("the symmetric drive", Settings { drive: Drive::Symmetric, ..Settings::default() }),
        ("at 10 Pa", Settings { pressure: 10.0, ..Settings::default() }),
    ] {
        let rig = Rig::new(settings);
        let mut s = State { on: true, ..rig.at_rest() };
        let (mut turn, mut swing) = (0.0f64, 0.0f64);
        for _ in 0..4000 {
            rig.step(&mut s);
            turn = turn.max(s.twist.abs());
            swing = swing.max(s.position[0].hypot(s.position[1]));
        }
        let gas = rig.breakdown().map_or("no breakdown".to_string(), |v| format!("breaks down past {v:.0} V"));
        println!("{name:<22} turns up to {:5.1}°, swings up to {:6.1} µm ({gas})", turn.to_degrees(), swing * 1e6);
    }

    // Zhang's vortex field with his three constants: in Earth's field alone, supply off.
    for &(name, k) in &ZHANG_KMQ {
        let rig = Rig::new(Settings { vortex: Some(vortex_constant(k)), ..Settings::default() });
        let r = rig.run(&mut rig.at_rest(), DT);
        println!("Zhang's vortex, k_mq {name}: {:.1e} N m, {:.1e} times what the thread resists per radian", r.vortex_torque, r.vortex_torque / rig.settings.torsion);
    }

    // The weighing rig (Zhang & Xu): a nylon ring (radii 20 and 50 mm, 50 mm tall) around a wire
    // with a 50 mm gap at 55 kV, reported 0.83–1.20 g lighter with the top live and 0.76–1.08 g
    // heavier the other way round. The wires are modelled as stacks of rings of their own radius.
    println!("the weighing rig, the polarisation pull on the nylon ring (reported: −1.0 g, +0.9 g):");
    for (name, top, bottom) in [("top live", 55e3, 0.0), ("bottom live", 0.0, 55e3), ("±27.5 kV", 27.5e3, -27.5e3)] {
        let lift = ring_pull(top, bottom) / G0 * 1e3; // grams lifted per unit of ε_r − 1
        let grams = |x: f64| if x.abs() < 5e-3 { 0.0 } else { -x }; // no "−0.00"
        println!("  {name:<12} weighs {:+.2} to {:+.2} g", grams(lift * 2.6), grams(lift * 0.722));
        if top == -bottom {
            assert!(lift.abs() < 1e-6); // the live-end test: symmetric drive, no pull
        } else {
            assert!(lift.abs() > 0.03 && lift.abs() < 0.2 && lift.signum() == if top > 0.0 { 1.0 } else { -1.0 });
        }
    }
}

fn print_forces(r: &Report, mass: f64) {
    let w = mass * G0;
    let size = |v: [f64; 3]| (v[0] * v[0] + v[1] * v[1] + v[2] * v[2]).sqrt() / w;
    println!("  every force at 3 s, as a fraction of the weight (field {:.3e} V/m, {:.0} V across the gas, coil current {:.1} µA, B {:.1} µT):", size(r.field) * w, r.across, r.current * 1e6, r.b[2] * 1e6);
    for (name, v) in [
        ("buoyancy", r.buoyancy),
        ("polarisation", r.polarisation),
        ("Coulomb", r.coulomb),
        ("charge dipole", r.dipole_force),
        ("image in the wall", r.image),
        ("ion drag", r.ion_drag),
        ("gas drag", r.gas_drag),
        ("Lorentz", r.lorentz),
        ("diamagnetism", r.magnetisation),
        ("field energy (GR)", r.field_gravity),
    ] {
        println!("    {name:<18} {:.1e}", size(v));
    }
    println!("  torques (N m): thread {:.1e}, damping {:.1e}, shape {:.1e}, charge patches {:.1e}, vortex {:.1e}; glow {} ({:.1e} A)", r.thread_torque, r.damping_torque, r.alignment, r.dipole_torque, r.vortex_torque, r.glow, r.glow_current);
    assert!(r.weight[2] < 0.0);
}

/// The upward polarisation force (N) on the nylon ring, `½ε₀(ε_r − 1)∫(E²_top − E²_bottom)dA`, per
/// unit of `ε_r − 1`: multiply by 2.6 for the dilute limit and by 0.72 (`1 − 1/ε_r`) for a field
/// normal to the ring's faces, the two ends of the range for nylon (ε_r = 3.6).
fn ring_pull(top: f64, bottom: f64) -> f64 {
    let (a, n) = (0.75e-3, 200);
    let mut rings: Vec<Ring> = stack(a, 0.025 + a, 2.0 * a, n);
    rings.extend(stack(a, -0.025 - a, -2.0 * a, n));
    let volts: Vec<f64> = (0..2 * n).map(|i| if i < n { top } else { bottom }).collect();
    let q = charges(&rings, a, &volts);
    let e2 = |rho: f64, z: f64| {
        let [x, y] = field(&rings, &q, rho, z);
        x * x + y * y
    };
    let steps = 60;
    let dr = 0.03 / steps as f64;
    let integral: f64 = (0..=steps)
        .map(|i| {
            let rho = 0.02 + dr * i as f64;
            let w = if i == 0 || i == steps { 0.5 } else { 1.0 };
            w * (e2(rho, 0.025) - e2(rho, -0.025)) * 2.0 * PI * rho * dr
        })
        .sum();
    0.5 * metrology::EPS_0 * integral
}

#[cfg(test)]
mod tests {
    use super::*;

    fn close(a: f64, b: f64, rel: f64) -> bool {
        (a / b - 1.0).abs() < rel
    }

    fn horizontal(p: [f64; 3]) -> f64 {
        p[0].hypot(p[1])
    }

    // Run with the supply on for `seconds`, returning the largest twist and swing seen.
    fn switch_on(settings: Settings, seconds: f64) -> (Rig, State, f64, f64) {
        let rig = Rig::new(settings);
        let mut s = State { on: true, ..rig.at_rest() };
        let (mut twist, mut swing) = (0.0f64, 0.0f64);
        for _ in 0..(seconds / DT) as usize {
            rig.step(&mut s);
            twist = twist.max(s.twist.abs());
            swing = swing.max(horizontal(s.position));
        }
        (rig, s, twist, swing)
    }

    #[test]
    fn the_field_matches_an_independent_solution() {
        // NumPy/SciPy, written separately (research/external/…-checks.py): 0.996 21, 0.525 03,
        // −5.229 82 V/m per volt at the ball. The leads tilt the field 12° off the vertical.
        let rig = Rig::new(Settings::default());
        let [x, y, z] = rig.unit_field([0.0; 3]);
        assert!(close(x, 0.996_212, 1e-5) && close(y, 0.525_028, 1e-5) && close(z, -5.229_824, 1e-5), "{x} {y} {z}");
        assert!(close(rig.capacitance(), 9.4622e-12, 1e-4) && close(rig.energy(), 4.7311e-12, 1e-4));
        assert!(rig.breakdown().is_none()); // 1 Pa: no Townsend breakdown fits in 10 cm
    }

    #[test]
    fn earth_field_at_zhangs_lab() {
        // WMM2025 at Lujiang: 50.0 µT, dipping 48.1°, declination 5.7° west.
        let [x, y, z] = earth_field();
        let h = x.hypot(y);
        assert!(close((h * h + z * z).sqrt(), 50.0e-6, 2e-3));
        assert!(((-z).atan2(h).to_degrees() - 48.1).abs() < 0.1);
        assert!(((-y).atan2(x).to_degrees() + 5.69).abs() < 0.05, "{}", (-y).atan2(x).to_degrees());
    }

    #[test]
    fn switched_off_nothing_moves() {
        let rig = Rig::new(Settings::default());
        let mut s = rig.at_rest();
        let r = rig.run(&mut s, 2.0);
        assert!(s.position == [0.0; 3] && s.twist == 0.0 && s.supply == 0.0);
        assert!(r.coulomb == [0.0; 3] && r.alignment == 0.0 && r.dipole_torque == 0.0);
    }

    #[test]
    fn a_charged_ball_turns_and_swings() {
        // 0.2 nC with a 0.1 nC patch: the patch's dipole lines up with the leads' sideways field,
        // rocking the ball through tens of degrees every 3 s; the charge swings it 0.2–0.3 mm.
        let (_, s, twist, swing) = switch_on(Settings::default(), 5.0);
        assert!(twist > 0.5 && twist < 1.5, "{twist}");
        assert!(swing > 1e-4 && swing < 5e-4, "{swing}");
        assert!(!s.glow);
        // The forces behind it: the Coulomb force is a few per cent of the weight, the
        // polarisation pull 20 times less; the patch's torque is 1 000 times the shape's.
        let rig = Rig::new(Settings::default());
        let mut s = State { on: true, ..rig.at_rest() };
        let r = rig.run(&mut s, 0.5);
        let w = -r.weight[2];
        assert!(close(r.coulomb[2] / w, -0.0152, 0.02), "{}", r.coulomb[2] / w);
        assert!(horizontal(r.polarisation) < horizontal(r.coulomb) / 10.0);
        assert!(r.dipole_torque.abs() > 100.0 * r.alignment.abs());
        // Gravity of the field energy: 10⁻³¹ N, 10⁻²⁸ of the weight; magnetism 10⁻¹⁸ N.
        assert!(r.field_gravity[2] / w < 1e-27 && r.lorentz.iter().all(|f| f.abs() < 1e-16));
    }

    #[test]
    fn a_neutral_round_ball_barely_moves() {
        let ball = Ball { charge: 0.0, dipole: 0.0, aspect: 1.0, ..Ball::polyethylene() };
        let (_, _, twist, swing) = switch_on(Settings { ball, ..Settings::default() }, 5.0);
        // Only the polarisation pull acts: a swing of about 10 µm and no turn at all.
        assert!(twist < 1e-12, "{twist}");
        assert!(swing > 3e-6 && swing < 3e-5, "{swing}");
    }

    #[test]
    fn a_glow_clamps_the_field_and_drains_the_ball() {
        // At 10 Pa the gas breaks down while the supply rises, the walls take up all but 574 V,
        // and the field at the ball falls from 2.7×10⁵ to 3×10³ V/m.
        let settings = Settings { pressure: 10.0, ..Settings::default() };
        let rig = Rig::new(settings.clone());
        assert!(close(rig.breakdown().unwrap(), 574.3, 1e-3), "{:?}", rig.breakdown());
        let (rig, s, _, _) = switch_on(settings, 2.0);
        let r = rig.run(&mut s.clone(), DT);
        assert!(close(r.across, 574.3, 1e-3) && !r.glow, "{}", r.across);
        assert!(close(norm2(r.field).sqrt(), 574.3 * 5.3497, 1e-3));
        // The glow drained the 0.2 nC to the floating charge of a 3 eV glow, −6.7 pC.
        let floating = 4.0 * PI * EPS_0 * Ball::polyethylene().radius * floating_potential(3.0, N2_ION_MASS);
        assert!(close(s.charge, floating, 1e-6) && close(floating, -6.68e-12, 1e-2), "{floating:e}");
    }

    #[test]
    fn the_symmetric_drive() {
        // ±25 kV: the same field at the ball within 1%, but the live side's pull drops and turns
        // over: the leads keep it from vanishing.
        let single = Rig::new(Settings::default());
        let sym = Rig::new(Settings { drive: Drive::Symmetric, ..Settings::default() });
        let (a, b) = (single.unit_field([0.0; 3]), sym.unit_field([0.0; 3]));
        assert!(close(a[2], b[2], 0.01));
        let pull = |rig: &Rig| {
            let mut s = State { on: true, ..State::default() };
            rig.run(&mut s, 0.5).polarisation[2]
        };
        let (ps, pm) = (pull(&single), pull(&sym));
        assert!(ps > 0.0 && pm < 0.0 && pm.abs() < ps, "{ps:e} {pm:e}");
    }

    #[test]
    fn zhangs_vortex_field_with_his_constants() {
        let f: Vec<f64> = ZHANG_KMQ.iter().map(|&(_, k)| vortex_constant(k)).collect();
        // With the "calibrated" constant the vortex twists the thread by 10⁻¹³ rad.
        let rig = Rig::new(Settings { vortex: Some(f[0]), ball: Ball { charge: 0.0, dipole: 0.0, aspect: 1.0, ..Ball::polyethylene() }, ..Settings::default() });
        let mut s = State { on: true, ..State::default() };
        let r = rig.run(&mut s, 0.1);
        assert!((r.vortex_torque / rig.settings.torsion).abs() < 1e-12);
        // With Paper E's, Earth's own field drives it with the supply off: 10²⁰ N m on a ball
        // whose thread resists 10⁻¹⁰ N m per radian.
        let rig = Rig::new(Settings { vortex: Some(f[1]), ..Settings::default() });
        let r = rig.step(&mut rig.at_rest());
        assert!(r.vortex_torque.abs() > 1e20, "{:e}", r.vortex_torque);
    }

    // web/gravity-patent.html runs the same model in the browser and checks itself against these
    // on load; change both together.
    #[test]
    fn reference_values_for_the_web_model() {
        let rig = Rig::new(Settings::default());
        let u = rig.unit_field([0.0; 3]);
        for (a, b) in u.iter().zip([0.996_211_955_909_649_1, 0.525_028_041_381_924, -5.229_823_828_310_781]) {
            assert!(close(*a, b, 1e-12));
        }
        assert!(close(rig.capacitance(), 9.462_167_488_509_768e-12, 1e-12) && close(rig.energy(), 4.731_083_744_254_884e-12, 1e-12));
        for (a, b) in rig.settings.earth_field.iter().zip([3.324_258_878_254_768e-5, 3.313_419_922_402_510_4e-6, -3.721_905_735_339_772_4e-5]) {
            assert!(close(*a, b, 1e-12));
        }
        assert!(close(rig.settings.torsion, 1.214_551_168_410_677_5e-10, 1e-12));
        assert!(close(rig.tip_field(), 363.016_540_111_795_3, 1e-9), "{}", rig.tip_field());
        let mut s = State { on: true, ..rig.at_rest() };
        rig.run(&mut s, 1.0);
        assert!(close(s.twist, 0.928_723_630_577_296_3, 1e-9) && close(s.spin, 0.429_074_213_809_598_7, 1e-9));
        assert!(close(s.position[0], 7.407_758_897_706_475e-5, 1e-9) && close(s.position[1], 4.120_759_664_600_089e-5, 1e-9));
        let rig = Rig::new(Settings { pressure: 10.0, ..Settings::default() });
        let mut s = State { on: true, ..rig.at_rest() };
        rig.run(&mut s, 1.0);
        assert!(close(s.wall, 49_425.712_880_056_68, 1e-12) && close(s.charge, -6.681_877_402_210_595e-12, 1e-9));
        assert!(close(s.twist, 2.804_066_201_676_511e-5, 1e-6));
    }

    #[test]
    fn the_ptfe_ball_is_the_mass_test() {
        // The same charge on the same size of ball: the denser ball is pushed the same, so it
        // accelerates 950/2200 as much. A gravitational field would accelerate both alike.
        let pe = Rig::new(Settings::default());
        let ptfe = Rig::new(Settings { ball: Ball::ptfe(), ..Settings::default() });
        let push = |rig: &Rig| {
            let mut s = State { on: true, supply: 50e3, charge: 0.2e-9, ..State::default() };
            rig.step(&mut s).coulomb[0] / rig.settings.ball.mass()
        };
        assert!(close(push(&ptfe) / push(&pe), 950.0 / 2200.0, 1e-9));
    }

    #[test]
    fn the_patent_rig_runs() {
        main();
    }

    #[test]
    fn in_two_seconds_the_ball_turns_and_swings() {
        let rig = Rig::new(Settings::default());
        let mut s = State { on: true, ..rig.at_rest() };
        let (mut turn, mut swing) = (0.0f64, 0.0f64);
        for _ in 0..2000 {
            let r = rig.step(&mut s);
            assert!(r.field[2] <= 0.0); // the field at the ball points down, toward the grounded coil
            turn = turn.max(s.twist.abs());
            swing = swing.max(s.position[0].hypot(s.position[1]));
        }
        assert!(turn.to_degrees() > 20.0 && swing < 1e-3);
    }
}