gravimetry 0.5.0

Satellite gravimetry: mission and sensor noise (electrostatic and cold-atom accelerometers, laser ranging, gradiometers), closed-loop simulation, GRACE processing (technical notes, smoothing, stripe removal, water height), earthquake gravity, a Kalman smoother and atmospheric waves.
Documentation
//! Cold-atom (light-pulse) interferometers as space accelerometers and gradiometers, at the
//! quantum-projection-noise limit.
//!
//! In microgravity the Mach–Zehnder π/2–π–π/2 sequence uses double diffraction, so the phase is
//! `Δφ = 2·k_eff·a·T²` (Lévèque et al. 2021, arXiv:2011.03382; Trimeche et al. 2019,
//! arXiv:1903.09828). Per-shot phase noise is `1/(C·√N)`, and the white noise density is the
//! per-shot noise times `√T_c`.

/// Two-photon Raman wave number `k_eff = 4π/λ` for ⁸⁷Rb on the D2 line (780.241 nm), m⁻¹.
pub const K_EFF_RB87: f64 = 4.0 * core::f64::consts::PI / 780.241e-9;

/// One atom interferometer (or each of a gradiometer's two).
#[derive(Clone, Copy, Debug, PartialEq)]
pub struct AtomInterferometer {
    /// Time between pulses `T`, s (the interferometer lasts `2T`).
    pub t: f64,
    /// Fringe contrast `C`, 0–1.
    pub contrast: f64,
    /// Atoms detected per shot `N`.
    pub atoms: f64,
    /// Cycle time `T_c`, s: `2T + T_prep` when shots run one after another, `T_prep` when
    /// interleaved.
    pub cycle: f64,
}

impl AtomInterferometer {
    /// Interferometer phase for acceleration `a` (m/s²), radians: `2·k_eff·a·T²`.
    pub fn phase(&self, a: f64) -> f64 {
        2.0 * K_EFF_RB87 * a * self.t * self.t
    }

    /// Quantum-projection phase noise per shot, `1/(C·√N)`, radians.
    pub fn phase_noise(&self) -> f64 {
        1.0 / (self.contrast * self.atoms.sqrt())
    }

    /// Acceleration noise per shot, m/s².
    pub fn acceleration_noise(&self) -> f64 {
        self.phase_noise() / (2.0 * K_EFF_RB87 * self.t * self.t)
    }

    /// White acceleration noise density, m s⁻² Hz⁻¹ᐟ²; use it as
    /// [`Accelerometer::AtomInterferometer`](crate::missions::Accelerometer::AtomInterferometer).
    pub fn acceleration_asd(&self) -> f64 {
        self.acceleration_noise() * self.cycle.sqrt()
    }

    /// Gravity-gradient noise per shot of a gradiometer made of two such interferometers a
    /// `baseline` (m) apart sharing the same lasers, s⁻²: `√2·σ_φ/(2·k_eff·D·T²)`.
    pub fn gradient_noise(&self, baseline: f64) -> f64 {
        std::f64::consts::SQRT_2 * self.phase_noise() / (2.0 * K_EFF_RB87 * baseline * self.t * self.t)
    }

    /// White gravity-gradient noise density of that gradiometer, s⁻² Hz⁻¹ᐟ².
    pub fn gradient_asd(&self, baseline: f64) -> f64 {
        self.gradient_noise(baseline) * self.cycle.sqrt()
    }
}

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

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

    /// Published instrument parameter sets, recomputed from first principles.
    #[test]
    fn published_instruments() {
        // CARIOQA (Lévèque et al. 2022): T = 1 s, C = 0.8, N = 1e5, T_c = 3.5 s. The paper quotes
        // 3e-10 (conservative) and 8e-11 with N = 1e6; QPN alone gives 2.3e-10 and 7.3e-11.
        let carioqa = AtomInterferometer { t: 1.0, contrast: 0.8, atoms: 1e5, cycle: 3.5 };
        assert!(close(carioqa.phase_noise(), 3.95e-3, 0.01));
        assert!(close(carioqa.acceleration_asd(), 2.3e-10, 0.01));
        assert!(close(AtomInterferometer { atoms: 1e6, ..carioqa }.acceleration_asd(), 7.3e-11, 0.01));
        // GRICE (Lévèque et al. 2021): T = 0.5 s, C = 0.43, N = 5e5, T_c = 2 s → 6e-10.
        let grice = AtomInterferometer { t: 0.5, contrast: 0.43, atoms: 5e5, cycle: 2.0 };
        assert!(close(grice.acceleration_asd(), 5.8e-10, 0.01));
        // Trimeche et al. 2019 gradiometer: T = 5 s, N = 1e6, ideal contrast, D = 0.5 m → 3.5 mE
        // per shot, and a 1087 rad differential phase for T_zz = 2.7e-6 s⁻².
        let tri = AtomInterferometer { t: 5.0, contrast: 1.0, atoms: 1e6, cycle: 1.0 };
        assert!(close(tri.gradient_noise(0.5), 3.5e-12, 0.01));
        assert!(close(tri.phase(2.7e-6 * 0.5), 1087.0, 0.001));
    }
}