spectroxide/
electron_temp.rs

1//! Electron temperature state (ρ_e = T_e/T_z).
2//!
3//! This module provides the [`ElectronTemperature`] struct used by the solver
4//! to hold the current T_e/T_z ratio. The production T_e update is performed
5//! in `solver::compute_hubble_coefficients` using a perturbative form of the
6//! Compton equilibrium:
7//!
8//!   Δρ_eq = ΔI₄/(4 G₃) − ΔG₃/G₃ × (I₄/(4G₃))
9//!
10//! computed from Δn only. The full form ρ_eq = I₄/(4G₃) (below) has
11//! ~0.1% numerical error from near-cancellation that swamps the O(10⁻⁵)
12//! physical signal — do not use it in the solver. It is retained here only
13//! as a verification tool for off-path consistency checks and tests.
14//!
15//! References:
16//! - Chluba & Sunyaev (2012), MNRAS 419, 1294 [Eq. 15-18]
17
18use crate::spectrum::compton_equilibrium_ratio;
19
20/// State of the electron temperature solver.
21#[derive(Debug, Clone)]
22pub struct ElectronTemperature {
23    /// Current T_e/T_z ratio
24    pub rho_e: f64,
25}
26
27impl Default for ElectronTemperature {
28    fn default() -> Self {
29        ElectronTemperature { rho_e: 1.0 }
30    }
31}
32
33impl ElectronTemperature {
34    /// θ_e from a precomputed θ_z value (cosmology-aware).
35    ///
36    /// Pass `cosmo.theta_z(z)` so a non-default T_CMB is honoured.
37    #[inline]
38    pub fn theta_e_with(&self, theta_z_val: f64) -> f64 {
39        self.rho_e * theta_z_val
40    }
41
42    /// Set ρ_e from the full Compton-equilibrium form I₄/(4G₃).
43    ///
44    /// **Not used by the production solver.** The full form has ~0.1%
45    /// numerical error from near-cancellation of the two integrals, which
46    /// swamps the O(10⁻⁵) physical distortion signal. Retained only as a
47    /// reference for off-path tests and verification; the solver uses the
48    /// perturbative update in `solver::compute_hubble_coefficients` instead.
49    pub fn update_equilibrium(&mut self, x_grid: &[f64], n_full: &[f64]) {
50        self.rho_e = compton_equilibrium_ratio(x_grid, n_full);
51    }
52}
53
54#[cfg(test)]
55mod tests {
56    use super::*;
57    use crate::grid::FrequencyGrid;
58    use crate::spectrum::planck;
59
60    #[test]
61    fn test_equilibrium_for_planck() {
62        let grid = FrequencyGrid::log_uniform(1e-4, 50.0, 5000);
63        let n_pl: Vec<f64> = grid.x.iter().map(|&x| planck(x)).collect();
64
65        let mut te = ElectronTemperature::default();
66        te.update_equilibrium(&grid.x, &n_pl);
67
68        assert!(
69            (te.rho_e - 1.0).abs() < 1e-3,
70            "ρ_e = {}, expected 1.0 for Planck",
71            te.rho_e
72        );
73    }
74
75    // test_theta_e_with_scaling removed: theta_e_with(θ_z) is defined as
76    // rho_e * θ_z, so asserting (1.05 * θ_z).abs() < 1e-30 was tautological.
77
78    /// Verify ρ_e for a Bose-Einstein distribution with known μ.
79    ///
80    /// For n_BE(x, μ) = 1/(e^{x+μ}-1), the Compton equilibrium ratio is:
81    ///   ρ_e = I₄/(4G₃) where I₄ = ∫x⁴ n(1+n)dx, G₃ = ∫x³ n dx.
82    /// For μ > 0: spectrum is harder than Planck → ρ_e > 1.
83    /// For μ < 0: spectrum is softer → ρ_e < 1.
84    #[test]
85    fn test_equilibrium_for_bose_einstein() {
86        let grid = FrequencyGrid::log_uniform(1e-4, 50.0, 10000);
87
88        // Only test positive μ: n_BE(x, μ) = 1/(e^{x+μ}-1) is well-defined for μ > 0.
89        // Negative μ causes a pole at x = |μ| which makes the integral diverge.
90        for &mu in &[1e-4, 1e-3, 5e-3] {
91            let n_be: Vec<f64> = grid
92                .x
93                .iter()
94                .map(|&x| 1.0 / ((x + mu).exp() - 1.0))
95                .collect();
96
97            let mut te = ElectronTemperature::default();
98            te.update_equilibrium(&grid.x, &n_be);
99
100            // Independent numerical integration (proper trapezoidal for n(1+n))
101            let mut g3 = 0.0;
102            let mut i4 = 0.0;
103            for i in 1..grid.n {
104                let dx = grid.x[i] - grid.x[i - 1];
105                let x_mid = 0.5 * (grid.x[i] + grid.x[i - 1]);
106                let n_mid = 0.5 * (n_be[i] + n_be[i - 1]);
107                let n_l = n_be[i - 1];
108                let n_r = n_be[i];
109                let nn1_mid = 0.5 * (n_l * (1.0 + n_l) + n_r * (1.0 + n_r));
110                g3 += x_mid.powi(3) * n_mid * dx;
111                i4 += x_mid.powi(4) * nn1_mid * dx;
112            }
113            let rho_expected = i4 / (4.0 * g3);
114
115            let rel = (te.rho_e - rho_expected).abs() / rho_expected;
116            assert!(
117                rel < 1e-6,
118                "μ={mu}: ρ_e={:.10} vs expected {rho_expected:.10}, err={rel:.2e}",
119                te.rho_e
120            );
121
122            // Physical direction check: μ > 0 → fewer low-x photons → harder spectrum → ρ_e > 1
123            assert!(
124                te.rho_e > 1.0,
125                "μ={mu}>0 should give ρ_e>1, got {:.10}",
126                te.rho_e
127            );
128
129            // Larger μ → more deviation from 1 (use tolerance; at small μ the signal is ~μ²)
130            if mu > 1e-4 {
131                let rho_small = {
132                    let n_small: Vec<f64> = grid
133                        .x
134                        .iter()
135                        .map(|&x| 1.0 / ((x + 1e-4).exp() - 1.0))
136                        .collect();
137                    crate::spectrum::compton_equilibrium_ratio(&grid.x, &n_small)
138                };
139                assert!(
140                    te.rho_e >= rho_small - 1e-9,
141                    "Larger μ should give larger ρ_e: μ={mu} → {:.10} < μ=1e-4 → {rho_small:.10}",
142                    te.rho_e
143                );
144            }
145        }
146    }
147}
spectroxide/electron_temp.rs

spectroxide/
electron_temp.rs
