Function decompose_nonlinear_be

pub fn decompose_nonlinear_be(
    x_grid: &[f64],
    delta_n: &[f64],
    x_min: f64,
    x_max: f64,
) -> DistortionParams

Bianchini & Fabbian (2022) nonlinear fit: μ inside the BE exponential.

Reference: Bianchini & Fabbian (2022), arXiv:2206.02762, Eqs. (1)–(4).

Model: Δn_model(x; μ, δ, y) = [n_pl(x/(1+δ)) − n_pl(x)] + [n_BE(x+μ) − n_pl(x)] + y · Y_SZ(x) with δ ≡ ΔT/T₀. Fits (μ, δ, y) by Levenberg-Marquardt on the band [x_min, x_max] with a trapezoidal inner product matching decompose_gram_schmidt.

Initial guess: bootstrap from decompose_gram_schmidt (converted via δ_BF = δ_GS + μ/β_μ). This gives the linearised optimum for free; the LM iterations only refine the O(μ²) nonlinear correction.

In the small-(μ, δ, y) limit the model reduces to a linear fit of Δn ≈ δ·G(x) + μ·(−G(x)/x) + y·Y_SZ(x), which spans the same 3-D subspace as the CJ2014 basis (Y_SZ, M, G) since M(x) = G(x)/β_μ − G(x)/x. The two methods therefore give the SAME μ and y to O(μ²), but a DIFFERENT ΔT/T: a pure B&F BE distortion with chemical potential μ_BF has ΔT/T = 0 in the B&F parameterisation and ΔT/T = −μ_BF/β_μ in CJ2014. Concretely: δ_BF = δ_CJ + μ/β_μ.