rexafs 0.2.5

//! Modified Bessel function of the first kind, order zero, for FFT windows.
//!
//! [`bessel_i0`] uses the Cephes polynomial coefficients and interval mapping
//! documented in the [SciPy Cephes source](https://github.com/scipy/scipy/blob/v1.15.3/scipy/special/xsf/cephes/i0.h).
//! That source attributes the approximation to Stephen L. Moshier, Cephes Math
//! Library release 2.8. Coefficients are retained here without refitting.
//! The function and its argument are dimensionless. It supplies the shape of
//! the Kaiser–Bessel window; it does not set the Fourier amplitude normalization.

/// Chebyshev coefficients for `exp(-x) * I0(x)` on `0 <= x <= 8`.
const A: &[f64] = &[
    -4.415_341_646_479_339_5E-18,
    3.330_794_518_822_238_4E-17,
    -2.431_279_846_547_955E-16,
    1.715_391_285_555_133E-15,
    -1.168_533_287_799_345_1E-14,
    7.676_185_498_604_936E-14,
    -4.856_446_783_111_929E-13,
    2.955_052_663_129_64E-12,
    -1.726_826_291_441_556E-11,
    9.675_809_035_373_237E-11,
    -5.189_795_601_635_263E-10,
    2.659_823_724_682_386_6E-9,
    -1.300_025_009_986_248E-8,
    6.046_995_022_541_919E-8,
    -2.670_793_853_940_612E-7,
    1.117_387_539_120_103_7E-6,
    -4.416_738_358_458_750_5E-6,
    1.644_844_807_072_889_6E-5,
    -5.754_195_010_082_104E-5,
    1.885_028_850_958_416_5E-4,
    -5.763_755_745_385_824E-4,
    1.639_475_616_941_335_7E-3,
    -4.324_309_995_050_576E-3,
    1.054_646_039_459_499_8E-2,
    -2.373_741_480_589_947E-2,
    4.930_528_423_967_071E-2,
    -9.490_109_704_804_764E-2,
    1.716_209_015_222_087_7E-1,
    -3.046_826_723_431_984E-1,
    6.767_952_744_094_761E-1,
];

/// Chebyshev coefficients for `exp(-x) * sqrt(x) * I0(x)` on `x >= 8`,
/// evaluated after the reciprocal mapping `32 / x - 2`.
const B: &[f64] = &[
    -7.233_180_487_874_754E-18,
    -4.830_504_485_944_182E-18,
    4.465_621_420_296_76E-17,
    3.461_222_867_697_461E-17,
    -2.827_623_980_516_583_6E-16,
    -3.425_485_619_677_219E-16,
    1.772_560_133_056_526_3E-15,
    3.811_680_669_352_622_4E-15,
    -9.554_846_698_828_307E-15,
    -4.150_569_347_287_222E-14,
    1.540_086_217_521_41E-14,
    3.852_778_382_742_142_6E-13,
    7.180_124_451_383_666E-13,
    -1.794_178_531_506_806_2E-12,
    -1.321_581_184_044_771_3E-11,
    -3.149_916_527_963_241_6E-11,
    1.188_914_710_784_643_9E-11,
    4.940_602_388_224_97E-10,
    3.396_232_025_708_386_5E-9,
    2.266_668_990_498_178E-8,
    2.048_918_589_469_063_8E-7,
    2.891_370_520_834_756_7E-6,
    6.889_758_346_916_825E-5,
    3.369_116_478_255_694_3E-3,
    8.044_904_110_141_088E-1,
];

/// Evaluate a Chebyshev series using the Cephes recurrence without allocation.
///
/// If `array` contains `[a[n], ..., a[0]]`, returns
/// `a[0] / 2 + sum_{j=1..n} a[j] * T[j](x / 2)`, where `T[j]` is the
/// degree-j Chebyshev polynomial. The argument is dimensionless; the output
/// shares the coefficient units. The reversed coefficient order, halved
/// constant term and `x / 2` argument follow
/// [Cephes chbevl](https://netlib.org/cephes/doubldoc.html#chbevl).
///
/// The approximation tables in this module map their domains to `-2 <= x <= 2`.
/// This helper itself does not restrict `x` or check finiteness; large values can
/// overflow. An empty coefficient slice panics.
pub fn chvevl(x: f64, array: &[f64]) -> f64 {
    let mut b0: f64 = array[0];
    let mut b1: f64 = 0.0;
    let mut b2: f64 = 0.0;

    for p0 in &array[1..] {
        b2 = b1;
        b1 = b0;
        b0 = x * b1 - b2 + *p0;
    }

    0.5 * (b0 - b2)
}

/// Approximate the dimensionless modified Bessel function `I0(x)` for real `x`.
///
/// The defining series is `sum_{j=0..infinity} (x*x/4)^j / (j!)^2`, with `j`
/// a nonnegative integer and `I0(0) = 1`; see
/// [DLMF equation 10.25.2](https://dlmf.nist.gov/10.25.E2) at order zero.
/// Evaluation uses the even symmetry `I0(-x) = I0(x)` and separate Chebyshev
/// approximations below and above `abs(x) = 8`, followed by exponential scaling.
/// Increasing a positive argument increases the Kaiser–Bessel shape contribution.
///
/// This is the unscaled function. Large finite arguments can overflow the
/// intermediate exponential even when the final mathematical value is finite;
/// NaN and infinite inputs produce NaN. No validation or error return is provided.
pub fn bessel_i0(mut x: f64) -> f64 {
    let y: f64;

    if x < 0.0 {
        x = -x;
    }

    if x <= 8.0 {
        y = (&x / 2.0) - 2.0;
        return x.exp() * chvevl(y, A);
    }

    x.exp() * chvevl(32.0 / x - 2.0, B) / x.sqrt()
}

#[cfg(test)]
mod tests {
    use itertools::Itertools;
    use std::vec;

    use super::*;
    use crate::xafs::tests::TEST_TOL;
    use approx::assert_abs_diff_eq;

    #[test]
    fn test_bessel_i0() {
        let x: Vec<f64> = (0..=10).map(|x| x as f64).collect();
        let y = x.iter().map(|x| bessel_i0(*x)).collect_vec();

        // Expected values from scipy.special.i0
        let expected = [
            1.0,
            1.2660658777520082,
            2.279585302336067,
            4.880792585865024,
            11.30192195213633,
            27.239871823604442,
            67.23440697647797,
            168.59390851028968,
            427.56411572180474,
            1093.5883545113745,
            2815.716628466254,
        ];

        y.iter()
            .zip(expected.iter())
            .for_each(|(y, expected)| assert_abs_diff_eq!(y, expected, epsilon = TEST_TOL));
    }
}