diff --git a/rascaline/src/math/bessel/airy.rs b/rascaline/src/math/bessel/airy.rs
new file mode 100644
index 000000000..bd12ff23b
--- /dev/null
+++ b/rascaline/src/math/bessel/airy.rs
@@ -0,0 +1,403 @@
+use crate::math::consts::SQRT_3;
+
+use super::polynomials::{eval_polynomial, eval_polynomial_1};
+
+/// Solution of the differential equation
+///
+/// ```text
+/// y"(x) = xy
+/// ```
+///
+/// The function returns `((Ai, Bi), (Ai'(x), Bi'(x)))` where Ai, and Bi are
+/// the two independent solutions and Ai'(x), Bi'(x) their first derivatives.
+///
+/// Evaluation is by power series summation for small x, by rational minimax
+/// approximations for large x.
+///
+/// ACCURACY:
+///
+/// Error criterion is absolute when function <= 1, relative when function > 1,
+/// except * denotes relative error criterion. For large negative x, the
+/// absolute error increases as x^1.5. For large positive x, the relative error
+/// increases as x^1.5.
+///
+/// Arithmetic  domain   function  # trials      peak         rms
+/// IEEE        -10, 0     Ai        10000       1.6e-15     2.7e-16
+/// IEEE          0, 10    Ai        10000       2.3e-14*    1.8e-15*
+/// IEEE        -10, 0     Ai'       10000       4.6e-15     7.6e-16
+/// IEEE          0, 10    Ai'       10000       1.8e-14*    1.5e-15*
+/// IEEE        -10, 10    Bi        30000       4.2e-15     5.3e-16
+/// IEEE        -10, 10    Bi'       30000       4.9e-15     7.3e-16
+#[allow(clippy::many_single_char_names, clippy::similar_names, clippy::too_many_lines)]
+pub fn airy(x: f64) -> ((f64, f64), (f64, f64)) {
+    let mut ai = 0.0;
+    let mut bi = 0.0;
+    let mut aip = 0.0;
+    let mut bip = 0.0;
+
+    let mut z;
+    let zz;
+    let mut t;
+    let mut f;
+    let mut g;
+    let mut uf;
+    let mut ug;
+    let mut k;
+    let zeta;
+    let theta;
+    let mut domain_flag= 0;
+
+    if x > 103.892 {
+        return ((0.0, 0.0), (std::f64::INFINITY, std::f64::INFINITY));
+    }
+
+    if x < -2.09 {
+        // domain_flag = 15;
+        t = f64::sqrt(-x);
+        zeta = -2.0 * x * t / 3.0;
+        t = f64::sqrt(t);
+        k = SQPII / t;
+        z = 1.0 / zeta;
+        zz = z * z;
+        uf = 1.0 + zz * eval_polynomial(zz, &AFN) / eval_polynomial_1(zz, &AFD);
+        ug = z * eval_polynomial(zz, &AGN) / eval_polynomial_1(zz, &AGD);
+        theta = zeta + 0.25 * std::f64::consts::PI;
+        f = f64::sin(theta);
+        g = f64::cos(theta);
+        ai = k * (f * uf - g * ug);
+        bi = k * (g * uf + f * ug);
+
+        uf = 1.0 + zz * eval_polynomial(zz, &APFN) / eval_polynomial_1(zz, &APFD);
+        ug = z * eval_polynomial(zz, &APGN) / eval_polynomial_1(zz, &APGD);
+        k = SQPII * t;
+        aip = -k * (g * uf + f * ug);
+        bip = k * (f * uf - g * ug);
+
+        return ((ai, bi), (aip, bip));
+    }
+
+    if x >= 2.09 {
+        domain_flag = 5;
+        t = f64::sqrt(x);
+        zeta = 2.0 * x * t / 3.0;
+        g = f64::exp(zeta);
+        t = f64::sqrt(t);
+        k = 2.0 * t * g;
+        z = 1.0 / zeta;
+        f = eval_polynomial(z, &AN) / eval_polynomial_1(z, &AD);
+        ai = SQPII * f / k;
+        k = -0.5 * SQPII * t / g;
+        f = eval_polynomial(z, &APN) / eval_polynomial_1(z, &APD);
+        aip = f * k;
+        if x > 8.3203353 {
+            f = z * eval_polynomial(z, &BN16) / eval_polynomial_1(z, &BD16);
+            k = SQPII * g;
+            bi = k * (1.0 + f) / t;
+            f = z * eval_polynomial(z, &BPPN) / eval_polynomial_1(z, &BPPD);
+            bip = k * t * (1.0 + f);
+
+            return ((ai, bi), (aip, bip));
+        }
+    }
+
+    f = 1.0;
+    g = x;
+    t = 1.0;
+    uf = 1.0;
+    ug = x;
+    k = 1.0;
+    z = x * x * x;
+    while t > f64::EPSILON {
+        uf *= z;
+        k += 1.0;
+        uf /= k;
+        ug *= z;
+        k += 1.0;
+        ug /= k;
+        uf /= k;
+        f += uf;
+        k += 1.0;
+        ug /= k;
+        g += ug;
+        t = f64::abs(uf / f);
+    }
+    uf = C1 * f;
+    ug = C2 * g;
+
+    if domain_flag & 1 == 0 {
+        ai = uf - ug;
+    }
+
+    if domain_flag & 2 == 0 {
+        bi = SQRT_3 * (uf + ug);
+    }
+
+    k = 4.0;
+    uf = x * x / 2.0;
+    ug = z / 3.0;
+    f = uf;
+    g = 1.0 + ug;
+    uf /= 3.0;
+    t = 1.0;
+
+    while t > f64::EPSILON {
+        uf *= z;
+        ug /= k;
+        k += 1.0;
+        ug *= z;
+        uf /= k;
+        f += uf;
+        k += 1.0;
+        ug /= k;
+        uf /= k;
+        g += ug;
+        k += 1.0;
+        t = f64::abs(ug / g);
+    }
+    uf = C1 * f;
+    ug = C2 * g;
+
+    if domain_flag & 4 == 0 {
+        aip = uf - ug;
+    }
+
+    if domain_flag & 8 == 0 {
+        bip = SQRT_3 * (uf + ug);
+    }
+
+    return ((ai, bi), (aip, bip));
+}
+
+const C1: f64 = 0.3550280538878172;
+const C2: f64 = 0.2588194037928068;
+const SQPII: f64 = 0.5641895835477563;
+
+static AN: [f64; 8] = [
+    3.46538101525629e-1,
+    1.2007595273964581e1,
+    7.627960536152345e1,
+    1.6808922493463058e2,
+    1.5975639135016442e2,
+    7.053609068404442e1,
+    1.4026469116338967e1,
+    1.0,
+];
+
+static AD: [f64; 8] = [
+    5.675945326387702e-1,
+    1.475625625848472e1,
+    8.451389701414746e1,
+    1.7731808814540045e2,
+    1.642346928715297e2,
+    7.147784008255756e1,
+    1.4095913560783403e1,
+    1.0,
+];
+
+static APN: [f64; 8] = [
+    6.137591848140358e-1,
+    1.4745467078775532e1,
+    8.20584123476061e1,
+    1.711847813609764e2,
+    1.593178471371418e2,
+    6.997785993301031e1,
+    1.3947085698048157e1,
+    1.0,
+];
+
+static APD: [f64; 8] = [
+    3.3420367774973697e-1,
+    1.1181029730615816e1,
+    7.1172735214786e1,
+    1.5877808437283832e2,
+    1.5320642747580922e2,
+    6.867523045927804e1,
+    1.3849863475825945e1,
+    1.0,
+];
+
+static BN16: [f64; 5] = [
+    -2.5324079586936415e-1,
+    5.752851673324674e-1,
+    -3.2990703687322537e-1,
+    6.444040689482e-2,
+    -3.8251954664133675e-3,
+];
+
+static BD16: [f64; 5] = [
+    -7.156850950540353,
+    1.0603958071566469e1,
+    -5.232466364712515,
+    9.573958643783839e-1,
+    -5.508281471635496e-2,
+];
+
+static BPPN: [f64; 5] = [
+    4.654611627746516e-1,
+    -1.0899217380049393,
+    6.38800117371828e-1,
+    -1.2684434955310292e-1,
+    7.624878443421098e-3,
+];
+
+static BPPD: [f64; 5] = [
+    -8.70622787633159,
+    1.3899316270455321e1,
+    -7.141161446164312,
+    1.340085959606805,
+    -7.84273211323342e-2,
+];
+
+static AFN: [f64; 9] = [
+    -1.316963234183318e-1,
+    -6.264565444319123e-1,
+    -6.931580360369335e-1,
+    -2.797799815451191e-1,
+    -4.919001326095003e-2,
+    -4.062659235948854e-3,
+    -1.592764962392621e-4,
+    -2.776491081552329e-6,
+    -1.6778769848911465e-8,
+];
+
+static AFD: [f64; 9] = [
+    1.3356042070655324e1,
+    3.2682503279522464e1,
+    2.6736704094149957e1,
+    9.187074029072596,
+    1.4752914677166642,
+    1.1568717379518804e-1,
+    4.402916416152112e-3,
+    7.547203482874142e-5,
+    4.5185009297058035e-7,
+];
+
+static AGN: [f64; 11] = [
+    1.973399320916857e-2,
+    3.9110302961568827e-1,
+    1.0657989759959559,
+    9.391692298166502e-1,
+    3.5146565610554764e-1,
+    6.338889196289255e-2,
+    5.858041130483885e-3,
+    2.82851600836737e-4,
+    6.98793669997261e-6,
+    8.117892395543892e-8,
+    3.415517847659236e-10,
+];
+
+static AGD: [f64; 10] = [
+    9.30892908077442,
+    1.9835292871831214e1,
+    1.5564662893286462e1,
+    5.476860694229755,
+    9.542936116189619e-1,
+    8.645808263523921e-2,
+    4.126565238242226e-3,
+    1.0125908511650914e-4,
+    1.1716673321441352e-6,
+    4.9183457006293e-9,
+];
+
+static APFN: [f64; 9] = [
+    1.8536562402253556e-1,
+    8.867121880525841e-1,
+    9.873919817473985e-1,
+    4.0124108231800376e-11,
+    7.103049262896312e-2,
+    5.906186579956618e-3,
+    2.330514094017768e-4,
+    4.087187782890355e-6,
+    2.4837993290044246e-8,
+];
+
+static APFD: [f64; 9] = [
+    1.4734585468750254e1,
+    3.754239334354896e1,
+    3.146577512030464e1,
+    1.0996912520729877e1,
+    1.788850547669994,
+    1.4173327575366262e-1,
+    5.44066067017226e-3,
+    9.394212906545112e-5,
+    5.65978713036027e-7,
+];
+
+static APGN: [f64; 11] = [
+    -3.556154290330823e-2,
+    -6.373115181294355e-1,
+    -1.7085673888431236,
+    -1.5022187211731663,
+    -5.636066658221027e-1,
+    -1.0210103112021689e-1,
+    -9.483966959614452e-3,
+    -4.6032530748678097e-4,
+    -1.1430083648451737e-5,
+    -1.3341551868554742e-7,
+    -5.638038339588935e-10,
+];
+
+static APGD: [f64; 10] = [
+    9.858658016961304,
+    2.1640186735658595e1,
+    1.731307763897494e1,
+    6.178721752808288,
+    1.088486943963215,
+    9.950055434408885e-2,
+    4.784681996838866e-33,
+    1.1815963332283862e-4,
+    1.3748067355421944e-6,
+    5.799125149291476e-9,
+];
+
+#[cfg(test)]
+mod tests {
+    use approx::assert_relative_eq;
+
+    use super::airy;
+
+    static REFERENCE: [(f64, f64, f64, f64, f64); 9] = [
+        (-652.891, 0.10447544538688257, 0.03927275745730001, -1.0034469627779596, 2.6695436038635454),
+        (-123.0, -0.12756243520101215, 0.11148479086544187, -1.236685519932118, -1.4145093682439087),
+        (-43.23, 0.21524897123688552, -0.045610160090116446, 0.3011299896521486, 1.414990730955133),
+        (-1.0, 0.5355608832923521, 0.1039973894969446, -0.01016056711664521, 0.5923756264227924),
+        (-3.674e-24, 0.3550280538878172, 0.6149266274460007, -0.2588194037928068, 0.4482883573538264),
+        (1.4725, 0.07446845457316596, 1.828120247981961, -0.10036966257430177, 1.8104588650504139),
+        (4.289, 0.0005198079279904233, 148.15162273358868, -0.0011049736723187855, 297.42956047164853),
+        (25.822, 1.2773348723736746e-39, 2.4520220221102355e37, -6.503130605285926e-39, 1.2436182167358716e38),
+        (103.6241, 3.434110379654982e-307, 4.552768382873564e304, -3.496612379694316e-306, 4.633433573716728e305),
+    ];
+
+    static ZEROS: [(f64, f64, f64, f64, f64); 4] = [
+        // Ai zeros
+        (-5.520559828095551, 2.313678943005095e-16, -0.36790153149695626, 0.8652040258941519, -0.016571236446578794),
+        (-7.944133587120853, -3.222967925030853e-17, -0.33600537065305697, 0.9473357094415678, -0.010554505816063523),
+        // Bi zeros
+        (-6.169852128310251, -0.35786068428672557, -2.4652322944876206e-16, -0.014444121615121822, -0.8894799014265397),
+        (-10.529913506705357, -0.31317702152506904, -6.278100478929636e-16, -0.007429478173120849, -1.016389659221249),
+    ];
+
+    #[test]
+    fn test_airy() {
+        for &(x, ai, bi, aip, bip) in &REFERENCE {
+            let (val, der) = airy(x);
+
+            let max_relative = if x < 2.0 {1e-10} else { 1e-5 };
+
+            assert_relative_eq!(val.0, ai, max_relative=max_relative);
+            assert_relative_eq!(val.1, bi, max_relative=max_relative);
+            assert_relative_eq!(der.0, aip, max_relative=max_relative);
+            assert_relative_eq!(der.1, bip, max_relative=max_relative);
+        }
+
+        for &(x, ai, bi, _aip, _bip) in &ZEROS {
+            let (val, _) = airy(x);
+
+            assert_relative_eq!(val.0, ai, max_relative=1e-9, epsilon=1e-10);
+            assert_relative_eq!(val.1, bi, max_relative=1e-9, epsilon=1e-10);
+            // derivatives near the zero are quite bad
+            // assert_relative_eq!(der.0, aip, max_relative=1);
+            // assert_relative_eq!(der.1, bip, max_relative=1);
+        }
+    }
+}
diff --git a/rascaline/src/math/bessel/iv.rs b/rascaline/src/math/bessel/iv.rs
new file mode 100644
index 000000000..2050221da
--- /dev/null
+++ b/rascaline/src/math/bessel/iv.rs
@@ -0,0 +1,615 @@
+use std::f64::consts::PI;
+
+use log::warn;
+
+use crate::math::gamma;
+
+
+/// Modified Bessel function of the first kind, non-integer order
+///
+/// Returns modified Bessel function of order v of the argument. If x is
+/// negative, v must be integer valued.
+#[allow(clippy::float_cmp)]
+pub fn bessel_iv(mut v: f64, x: f64) -> f64 {
+    let mut sign;
+    let mut t;
+    let mut res = 0.0;
+
+    if !(x.is_finite() && v.is_finite()) {
+        return std::f64::NAN;
+    }
+
+    t = f64::floor(v);
+    if v < 0.0 && t == v {
+        v = -v;
+        t = -t;
+    }
+
+    sign = 1;
+    if x < 0.0 {
+        if t != v {
+            return std::f64::NAN;
+        }
+        if v != 2.0 * f64::floor(v / 2.0) {
+            sign = -1;
+        }
+    }
+    if x == 0.0 {
+        if v == 0.0 {
+            return 1.0;
+        }
+        if v < 0.0 {
+            return std::f64::INFINITY;
+        }
+        return 0.0;
+    }
+    let ax = f64::abs(x);
+    if f64::abs(v) > 50.0 {
+        ikv_asymptotic_uniform(v, ax, &mut res, &mut 0.0);
+    } else {
+        res = ikv_temme(v, ax);
+    }
+    res *= sign as f64;
+    return res;
+}
+
+fn iv_asymptotic(v: f64, x: f64) -> f64 {
+    let mut sum;
+    let mut term;
+    let mut factor;
+    let mut k;
+
+    let prefactor = f64::exp(x) / f64::sqrt(2.0 * PI * x);
+    if prefactor == std::f64::INFINITY {
+        return prefactor;
+    }
+    let mu = 4.0 * v * v;
+    sum = 1.0;
+    term = 1.0;
+    k = 1;
+    loop {
+        factor = (mu - ((2 * k - 1) * (2 * k - 1)) as f64) / (8.0 * x) / k as f64;
+        if k > 100 {
+            warn!("failed to converge bessel_iv function");
+            break;
+        }
+
+        term *= -factor;
+        sum += term;
+        k += 1;
+        if f64::EPSILON * f64::abs(sum) > f64::abs(term) {
+            break;
+        }
+    }
+    return sum * prefactor;
+}
+
+/// Uniform asymptotic expansion factors, (AMS5 9.3.9; AMS5 9.3.10)
+static ASYMPTOTIC_UFACTORS: [[f64; 31]; 11] = [
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 1.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -0.20833333333333334,
+        0.0, 0.125, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.3342013888888889, 0.0,
+        -0.4010416666666667, 0.0, 0.0703125, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1.0258125964506173, 0.0, 1.8464626736111112,
+        0.0, -0.8912109375, 0.0, 0.0732421875, 0.0, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        0.0, 0.0, 0.0, 4.669584423426247, 0.0, -11.207002616222995, 0.0,
+        8.78912353515625, 0.0, -2.3640869140625, 0.0, 0.112152099609375, 0.0, 0.0, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        -28.212072558200244, 0.0, 84.63621767460074, 0.0, -91.81824154324003,
+        0.0, 42.53499874538846, 0.0, -7.368794359479631, 0.0, 0.22710800170898438,
+        0.0, 0.0, 0.0, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 212.5701300392171,
+        0.0, -765.2524681411816, 0.0, 1059.9904525279999, 0.0, -699.5796273761327,
+        0.0, 218.1905117442116, 0.0, -26.491430486951554, 0.0, 0.5725014209747314,
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, -1919.4576623184068, 0.0,
+        8061.722181737308, 0.0, -13586.550006434136, 0.0, 11655.393336864536,
+        0.0, -5305.646978613405, 0.0, 1200.9029132163525, 0.0, -108.09091978839464,
+        0.0, 1.7277275025844574, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 20204.29133096615, 0.0, -96980.5983886375,
+        0.0, 192547.0012325315, 0.0, -203400.17728041555, 0.0, 122200.46498301747,
+        0.0, -41192.65496889756, 0.0, 7109.514302489364, 0.0, -493.915304773088,
+        0.0, 6.074042001273483, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    ],
+    [
+        0.0, 0.0, 0.0, -242919.18790055133, 0.0, 1311763.614662977, 0.0, -2998015.918538106,
+        0.0, 3763271.297656404, 0.0, -2813563.226586534, 0.0, 1268365.2733216248, 0.0,
+        -331645.1724845636, 0.0, 45218.76898136274, 0.0, -2499.830481811209, 0.0,
+        24.380529699556064, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+    ],
+    [
+        3284469.8530720375, 0.0, -19706819.11843222, 0.0, 50952602.49266463, 0.0,
+        -74105148.21153264, 0.0, 66344512.27472903, 0.0, -37567176.66076335, 0.0,
+        13288767.16642182, 0.0, -2785618.128086455, 0.0, 308186.40461266245, 0.0,
+        -13886.089753717039, 0.0, 110.01714026924674, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0,
+        0.0, 0.0, 0.0, 0.0,
+    ],
+];
+
+/// Compute Iv, Kv from (AMS5 9.7.7 + 9.7.8), asymptotic expansion for large v
+#[allow(clippy::many_single_char_names)]
+fn ikv_asymptotic_uniform(mut v: f64, x: f64, i_value: &mut f64, k_value: &mut f64) {
+    let mut term = 0.0;
+    let mut divisor;
+    let mut k;
+    let mut n;
+    let mut sign: i32 = 1;
+
+    if v < 0.0 {
+        sign = -1;
+        v = -v;
+    }
+
+    let z = x / v;
+    let t = 1.0 / f64::sqrt(1.0 + z * z);
+    let t2 = t * t;
+    let eta = f64::sqrt(1.0 + z * z) + f64::ln(z / (1.0 + 1.0 / t));
+    let i_prefactor = f64::sqrt(t / (2.0 * PI * v)) * f64::exp(v * eta);
+    let mut i_sum = 1.0;
+    let k_prefactor = f64::sqrt(PI * t / (2.0 * v)) * f64::exp(-v * eta);
+    let mut k_sum = 1.0;
+    divisor = v;
+    n = 1;
+    while n < 11 {
+        term = 0.0;
+        k = 31 - 1 - 3 * n;
+        while k < 31 - n {
+            term *= t2;
+            term += ASYMPTOTIC_UFACTORS[n as usize][k as usize];
+            k += 2;
+        }
+        k = 1;
+        while k < n {
+            term *= t2;
+            k += 2;
+        }
+        if n % 2 == 1 {
+            term *= t;
+        }
+        term /= divisor;
+        i_sum += term;
+        k_sum += if n % 2 == 0 { term } else { -term };
+        if f64::abs(term) < f64::EPSILON {
+            break;
+        }
+        divisor *= v;
+        n += 1;
+    }
+    if f64::abs(term) > 1e-3 * f64::abs(i_sum) {
+        // sf_error(SF_ERROR_NO_RESULT);
+        warn!("failed to converge bessel_iv function");
+    }
+
+    if f64::abs(term) > f64::EPSILON * f64::abs(i_sum) {
+        // sf_error(SF_ERROR_LOSS);
+        warn!("failed to converge bessel_iv function");
+    }
+
+    *k_value = k_prefactor * k_sum;
+
+    if sign == 1 {
+        *i_value = i_prefactor * i_sum;
+    } else {
+        *i_value = i_prefactor * i_sum
+            + 2.0 / PI * f64::sin(PI * v) * k_prefactor * k_sum;
+    }
+}
+
+/// Modified Bessel functions of the first and second kind of fractional order
+///
+/// Calculate K(v, x) and K(v+1, x) by method analogous to
+/// Temme, Journal of Computational Physics, vol 21, 343 (1976)
+#[allow(clippy::many_single_char_names)]
+fn temme_ik_series(v: f64, x: f64, k: &mut f64, k1: &mut f64) -> i32 {
+    let mut f;
+    let mut h;
+    let mut p;
+    let mut q;
+    let mut coef;
+    let mut sum;
+    let mut sum1;
+
+    /*
+     * |x| <= 2, Temme series converge rapidly
+     * |x| > 2, the larger the |x|, the slower the convergence
+     */
+
+    let gp = gamma(v + 1.0) - 1.0;
+    let gm = gamma(-v + 1.0) - 1.0;
+    let a = f64::ln(x / 2.0);
+    let b = f64::exp(v * a);
+    let sigma = -a * v;
+    let c = if f64::abs(v) < f64::EPSILON {
+        1.0
+    } else {
+        f64::sin(PI * v) / (v * PI)
+    };
+
+    let d = if f64::abs(sigma) < f64::EPSILON {
+        1.0
+    } else {
+        f64::sinh(sigma) / sigma
+    };
+
+    let gamma1 = if f64::abs(v) < f64::EPSILON {
+        -0.5772156649015329
+    } else {
+        0.5f32 as f64 / v * (gp - gm) * c
+    };
+
+    let gamma2 = (2.0 + gp + gm) * c / 2.0;
+
+    /* initial values */
+    p = (gp + 1.0) / (2.0 * b);
+    q = (1.0 + gm) * b / 2.0;
+    f = (f64::cosh(sigma) * gamma1 + d * -a * gamma2) / c;
+    h = p;
+    coef = 1.0;
+    sum = coef * f;
+    sum1 = coef * h;
+
+    let mut ik = 1.0;
+    while ik < 500.0 {
+        f = (ik * f + p + q) / (ik * ik - v * v);
+        p /= ik - v;
+        q /= ik + v;
+        h = p - ik * f;
+        coef *= x * x / (4.0 * ik);
+        sum += coef * f;
+        sum1 += coef * h;
+        if f64::abs(coef * f) < f64::abs(sum) * f64::EPSILON {
+            break;
+        }
+        ik += 1.0;
+    }
+
+    if ik >= 500.0 {
+        warn!("failed to converge bessel_iv function");
+    }
+
+    *k = sum;
+    *k1 = 2.0 * sum1 / x;
+    return 0;
+}
+
+/// Evaluate continued fraction `fv = I_(v+1) / I_v`, derived from
+/// Abramowitz and Stegun, Handbook of Mathematical Functions, 1972, 9.1.73
+#[allow(clippy::many_single_char_names)]
+fn cf1_ik(v: f64, x: f64, fv: &mut f64) -> i32 {
+    let mut c;
+    let mut d;
+    let mut f;
+    let mut a;
+    let mut b;
+    let mut delta;
+
+    /*
+     * |x| <= |v|, CF1_ik converges rapidly
+     * |x| > |v|, CF1_ik needs O(|x|) iterations to converge
+     */
+
+    /*
+     * modified Lentz's method, see
+     * Lentz, Applied Optics, vol 15, 668 (1976)
+     */
+
+    let tiny = 1.0 / f64::sqrt(1.7976931348623157e308);
+
+    f = tiny;
+    c = f;
+    d = 0.0;
+
+    let mut k = 1;
+    while k < 500 {
+        a = 1.0;
+        b = 2.0 * (v + k as f64) / x;
+        c = b + a / c;
+        d = b + a * d;
+        if c == 0.0 {
+            c = tiny;
+        }
+        if d == 0.0 {
+            d = tiny;
+        }
+        d = 1.0 / d;
+        delta = c * d;
+        f *= delta;
+        if f64::abs(delta - 1.0) <= 2.0 * f64::EPSILON {
+            break;
+        }
+        k += 1;
+    }
+
+    if k == 500 {
+        warn!("failed to converge bessel_iv function");
+    }
+
+    *fv = f;
+    return 0;
+}
+
+
+/// Calculate K(v, x) and K(v+1, x) by evaluating continued fraction
+/// `z1 / z0 = U(v+1.5, 2v+1, 2x) / U(v+0.5, 2v+1, 2x)`, see
+/// Thompson and Barnett, Computer Physics Communications, vol 47, 245 (1987)
+#[allow(non_snake_case, clippy::many_single_char_names)]
+fn cf2_ik(v: f64, x: f64, Kv: &mut f64, Kv1: &mut f64) -> i32 {
+    let mut S;
+    let mut C;
+    let mut Q;
+    let mut D;
+    let mut f;
+    let mut a;
+    let mut b;
+    let mut q;
+    let mut delta;
+
+    let mut current;
+    let mut prev;
+
+    /*
+     * |x| >= |v|, CF2_ik converges rapidly
+     * |x| -> 0, CF2_ik fails to converge
+     */
+
+    /*
+     * Steed's algorithm, see Thompson and Barnett,
+     * Journal of Computational Physics, vol 64, 490 (1986)
+     */
+
+    a = v * v - 0.25f32 as f64;
+    b = 2.0 * (x + 1.0);
+    D = 1.0 / b;
+    delta = D;
+    f = delta;
+    prev = 0.0;
+    current = 1.0;
+    C = -a;
+    Q = C;
+    S = 1.0 + Q * delta;
+
+    let mut ik = 2_u64;
+    while ik < 500 {
+        a -= (2 * (ik - 1)) as f64;
+        b += 2.0;
+        D = 1.0 / (b + a * D);
+        delta *= b * D - 1.0;
+        f += delta;
+        q = (prev - (b - 2.0) * current) / a;
+        prev = current;
+        current = q;
+        C *= -a / ik as f64;
+        Q += C * q;
+        S += Q * delta;
+        if f64::abs(Q * delta) < f64::abs(S) * f64::EPSILON {
+            break;
+        }
+        ik += 1;
+    }
+    if ik == 500 {
+        warn!("failed to converge bessel_iv function");
+    }
+    *Kv = f64::sqrt(PI / (2.0 * x)) * f64::exp(-x) / S;
+    *Kv1 = *Kv * (0.5f32 as f64 + v + x + (v * v - 0.25f32 as f64) * f) / x;
+    return 0;
+}
+
+
+/// Compute I(v, x) and K(v, x) simultaneously by Temme's method, see
+/// Temme, Journal of Computational Physics, vol 19, 324 (1975)
+#[allow(non_snake_case, clippy::many_single_char_names)]
+fn ikv_temme(mut v: f64, x: f64) -> f64 {
+    // modification from scipy's cephes: this only returns Iv, never Kv
+
+    let mut Iv;
+    let mut Ku = 0.0;
+    let mut Ku1 = 0.0;
+    let mut fv = 0.0;
+    let mut current;
+    let mut prev;
+    let mut next;
+    let mut reflect = false;
+
+    if v < 0.0 {
+        v = -v;
+        reflect = true;
+    }
+
+    let n = f64::round(v) as u32;
+    let u = v - n as f64;
+    if x < 0.0 {
+        return std::f64::NAN;
+    }
+
+    if x == 0.0 {
+        Iv = if v == 0.0 { 1.0 } else { 0.0 };
+
+        if reflect {
+            let z = u + n.wrapping_rem(2) as f64;
+            Iv = if f64::sin(PI * z) == 0.0 {
+                Iv
+            } else {
+                std::f64::INFINITY
+            };
+        }
+
+        return Iv;
+    }
+
+    /* x is positive until reflection */
+    let W = 1.0 / x; /* Wronskian */
+    if x <= 2.0 {
+        temme_ik_series(u, x, &mut Ku, &mut Ku1); /* Temme series */
+    } else {
+        /* x in (2, \infty) */
+        cf2_ik(u, x, &mut Ku, &mut Ku1); /* continued fraction CF2_ik */
+    }
+    prev = Ku;
+    current = Ku1;
+
+    for ik in 1..=n { /* forward recurrence for K */
+        next = 2.0 * (u + ik as f64) * current / x + prev;
+        prev = current;
+        current = next;
+    }
+    let Kv = prev;
+    let Kv1 = current;
+
+    let mut lim: f64 = (4.0 * v * v + 10.0) / (8.0 * x);
+    lim *= lim;
+    lim *= lim;
+    lim /= 24.0;
+    if lim < f64::EPSILON * 10.0 && x > 100.0 {
+        /*
+	     * x is huge compared to v, CF1 may be very slow
+	     * to converge so use asymptotic expansion for large
+	     * x case instead.  Note that the asymptotic expansion
+	     * isn't very accurate - so it's deliberately very hard
+	     * to get here - probably we're going to overflow:
+	     */
+        Iv = iv_asymptotic(v, x);
+    } else {
+        cf1_ik(v, x, &mut fv); /* continued fraction CF1_ik */
+        Iv = W / (Kv * fv + Kv1); /* Wronskian relation */
+    }
+
+    if reflect {
+        let z_0 = u + n.wrapping_rem(2) as f64;
+        return Iv + 2.0 / PI * f64::sin(PI * z_0) * Kv; /* reflection formula */
+    }
+
+    return Iv;
+}
+
+
+#[cfg(test)]
+mod tests {
+    use approx::assert_relative_eq;
+
+    // computed using mpmath, converting overflow (>1e308) to inf and underflow
+    // (<1e-323) to 0.0
+    static REFERENCE: [(f64, f64, f64); 92] = [
+        (-100.3, 1.0, 1.48788821588229e+186),
+        (-100.3, 10.0, 5.81385536644804e+85),
+        (-100.3, 200.5, 6.58080030590051e+74),
+        (-100.3, 401.0, 1.05963257281431e+167),
+        (-100.3, 600.5, 2.36123204175185e+255),
+        (-100.3, 700.6, 2.13787442407607e+299),
+        (-100.3, 1300.0, std::f64::INFINITY),
+        (-20.0, 1.0, 3.96683598581902e-25),
+        (-20.0, 10.0, 0.000125079973564495),
+        (-20.0, 200.5, 1.23661306184996e+85),
+        (-20.0, 401.0, 1.71683807030892e+172),
+        (-20.0, 600.5, 7.25819963948107e+258),
+        (-20.0, 700.6, 2.09367143600173e+302),
+        (-20.0, 1300.0, std::f64::INFINITY),
+        (-10.0, 1.0, 2.75294803983687e-10),
+        (-10.0, 10.0, 21.8917061637234),
+        (-10.0, 200.5, 2.6158708155538e+85),
+        (-10.0, 401.0, 2.49657353350645e+172),
+        (-10.0, 600.5, 9.31944446151452e+258),
+        (-10.0, 700.6, 2.59388351323027e+302),
+        (-10.0, 1300.0, std::f64::INFINITY),
+        (-1.0, 1.0, 0.565159103992485),
+        (-1.0, 10.0, 2670.98830370125),
+        (-1.0, 200.5, 3.35028839470096e+85),
+        (-1.0, 401.0, 2.82500017764435e+172),
+        (-1.0, 600.5, 1.0120885431321e+259),
+        (-1.0, 700.6, 2.78391772501927e+302),
+        (-1.0, 1300.0, std::f64::INFINITY),
+        (-0.5, 1.0, 1.23120021459297),
+        (-0.5, 10.0, 2778.78461532957),
+        (-0.5, 200.5, 3.35657610792771e+85),
+        (-0.5, 401.0, 2.82764655074214e+172),
+        (-0.5, 600.5, 1.01272129651641e+259),
+        (-0.5, 700.6, 2.78540929649693e+302),
+        (-0.5, 1300.0, std::f64::INFINITY),
+        (0.0, -1300.0, std::f64::INFINITY),
+        (0.0, -11.0, 7288.48933982125),
+        (0.0, -10.0, 2815.71662846625),
+        (0.0, -1.0, 1.26606587775201),
+        (0.0, 1.0, 1.26606587775201),
+        (0.0, 10.0, 2815.71662846625),
+        (0.0, 200.5, 3.35867463800083e+85),
+        (0.0, 401.0, 2.82852922635177e+172),
+        (0.0, 600.5, 1.01293230225787e+259),
+        (0.0, 700.6, 2.7859066646434e+302),
+        (0.0, 1300.0, std::f64::INFINITY),
+        (1.0, -1300.0, -std::f64::INFINITY),
+        (1.0, -11.0, -6948.85865981216),
+        (1.0, -10.0, -2670.98830370125),
+        (1.0, -1.0, -0.565159103992485),
+        (1.0, 1.0, 0.565159103992485),
+        (1.0, 10.0, 2670.98830370125),
+        (1.0, 200.5, 3.35028839470096e+85),
+        (1.0, 401.0, 2.82500017764435e+172),
+        (1.0, 600.5, 1.0120885431321e+259),
+        (1.0, 700.6, 2.78391772501927e+302),
+        (1.0, 1300.0, std::f64::INFINITY),
+        (12.49, 1.0, 1.06214483463358e-13),
+        (12.49, 10.0, 1.85487470997247),
+        (12.49, 200.5, 2.27429677988355e+85),
+        (12.49, 401.0, 2.3280143028664e+172),
+        (12.49, 600.5, 8.89455243908991e+258),
+        (12.49, 700.6, 2.49219423462741e+302),
+        (12.49, 1300.0, std::f64::INFINITY),
+        (120.0, -1300.0, std::f64::INFINITY),
+        (120.0, -11.0, 1.33841738592428e-110),
+        (120.0, -10.0, 1.38248722487846e-115),
+        (120.0, -1.0, 1.12694823613939e-235),
+        (120.0, 1.0, 1.12694823613939e-235),
+        (120.0, 10.0, 1.38248722487846e-115),
+        (120.0, 200.5, 2.08641377637118e+70),
+        (120.0, 401.0, 5.02466934582507e+164),
+        (120.0, 600.5, 6.47517615090698e+253),
+        (120.0, 700.6, 9.75968494448862e+297),
+        (120.0, 1300.0, std::f64::INFINITY),
+        (301.0, -1300.0, -std::f64::INFINITY),
+        (301.0, -11.0, 0.0),
+        (301.0, -10.0, 0.0),
+        (301.0, -1.0, 0.0),
+        (301.0, 1.0, 0.0),
+        (301.0, 10.0, 0.0),
+        (301.0, 200.5, 0.131118954561987),
+        (301.0, 401.0, 2.14246804623493e+125),
+        (301.0, 600.5, 7.21577161334608e+226),
+        (301.0, 700.6, 5.68946923641546e+274),
+        (301.0, 1300.0, std::f64::INFINITY),
+        (-59.5, 3.5, -1.88381631760945e+64),
+        (-39.5, 3.5, -2.40435409707822e+35),
+        (-19.5, 3.5, -136430465086.747),
+        (0.5, 3.5, 7.0552194086912),
+        (20.5, 3.5, 9.98384982441339e-15),
+        (40.5, 3.5, 1.43995556110118e-39),
+    ];
+
+    #[test]
+    fn test_bessel_iv() {
+        for (v, z, expected) in REFERENCE {
+            assert_relative_eq!(super::bessel_iv(v, z), expected, max_relative=1e-12);
+        }
+    }
+}
diff --git a/rascaline/src/math/bessel/j0.rs b/rascaline/src/math/bessel/j0.rs
new file mode 100644
index 000000000..762b842b6
--- /dev/null
+++ b/rascaline/src/math/bessel/j0.rs
@@ -0,0 +1,224 @@
+use super::{eval_polynomial, eval_polynomial_1};
+
+/* 5.783185962946784521175995758455807035071 */
+const DR1: f64 = 5.783185962946784;
+/* 30.47126234366208639907816317502275584842 */
+const DR2: f64 = 30.471262343662087;
+
+use std::f64::consts::{FRAC_PI_4, PI};
+use crate::math::consts::SQRT_FRAC_2_PI;
+
+/// Bessel function of the first kind, order zero.
+///
+/// The domain is divided into the intervals [0, 5] and
+/// (5, infinity). In the first interval the following rational
+/// approximation is used:
+///
+/// ```text
+///        2         2
+/// (w - r  ) (w - r  ) P (w) / Q (w)
+///       1         2    3       8
+/// ```
+///
+/// where w = x^2  and the two r's are zeros of the function.
+///
+/// In the second interval, the Hankel asymptotic expansion
+/// is employed with two rational functions of degree 6/6
+/// and 7/7.
+///
+/// # Accuracy
+/// ```text
+///                      Absolute error:
+/// arithmetic   domain     # trials      peak         rms
+///    IEEE      0, 30       60000       4.2e-16     1.1e-16
+/// ```
+pub fn bessel_j0(mut x: f64) -> f64 {
+    if x < 0.0 {
+        x = -x;
+    }
+
+    if x <= 5.0 {
+        let z = x * x;
+        if x < 1e-5 {
+            return 1.0 - z / 4.0;
+        }
+        let p = (z - DR1) * (z - DR2);
+        return p * eval_polynomial(z, &RP) / eval_polynomial_1(z, &RQ);
+    }
+
+    let w = 5.0 / x;
+    let q = 25.0 / (x * x);
+    let p = eval_polynomial(q, &PP) / eval_polynomial(q, &PQ);
+    let q = eval_polynomial(q, &QP) / eval_polynomial_1(q, &QQ);
+    let xn = x - FRAC_PI_4;
+    let p = p * f64::cos(xn) - w * q * f64::sin(xn);
+    return p * SQRT_FRAC_2_PI / f64::sqrt(x);
+}
+
+
+/// Bessel function of second kind, order zero
+///
+/// The domain is divided into the intervals [0, 5] and
+/// (5, infinity). In the first interval a rational approximation
+/// R(x) is employed to compute
+/// ```text
+///   y0(x) = R(x) + 2 * log(x) * j0(x) / NPY_PI.
+/// ```
+/// Thus a call to `bessel_j0()` is required.
+///
+/// In the second interval, the Hankel asymptotic expansion
+/// is employed with two rational functions of degree 6/6
+/// and 7/7.
+///
+/// # Accuracy
+///
+/// ```text
+///  Absolute error, when y0(x) < 1; else relative error:
+///
+/// arithmetic   domain     # trials      peak         rms
+///    IEEE      0, 30       30000       1.3e-15     1.6e-16
+/// ```
+#[allow(clippy::many_single_char_names)]
+pub fn bessel_y0(x: f64) -> f64 {
+    // Rational approximation coefficients YP[], YQ[] are used here. The
+    // function computed is  y0(x) - 2 * log(x) * j0(x) / PI, whose value at
+    // x = 0 is  2 * ( log(0.5) + EUL ) / PI = 0.073804295108687225.
+
+    if x == 0.0 {
+        // sf_error(SF_ERROR_SINGULAR);
+        return -std::f64::INFINITY;
+    } else if x < 0.0 {
+        // sf_error(SF_ERROR_DOMAIN);
+        return std::f64::NAN;
+    }
+
+    if x <= 5.0 {
+        let z = x * x;
+        let mut w = eval_polynomial(z, &YP) / eval_polynomial_1(z, &YQ);
+        w += 2.0 / PI * f64::ln(x) * bessel_j0(x);
+        return w;
+    }
+
+    let w = 5.0 / x;
+    let z = 25.0 / (x * x);
+    let p = eval_polynomial(z, &PP) / eval_polynomial(z, &PQ);
+    let q = eval_polynomial(z, &QP) / eval_polynomial_1(z, &QQ);
+    let xn = x - FRAC_PI_4;
+    let p = p * f64::sin(xn) + w * q * f64::cos(xn);
+    return p * SQRT_FRAC_2_PI / f64::sqrt(x);
+}
+
+
+// Note: all coefficients satisfy the relative error criterion
+// except YP, YQ which are designed for absolute error.
+
+static RP: [f64; 4] = [
+    -4.794432209782018e9,
+    1.9561749194655657e12,
+    -2.4924834436096772e14,
+    9.708622510473064e15,
+];
+static RQ: [f64; 8] = [
+    4.99563147152651e2,
+    1.737854016763747e5,
+    4.844096583399621e7,
+    1.1185553704535683e10,
+    2.112775201154892e12,
+    3.1051822985742256e14,
+    3.1812195594320496e16,
+    1.7108629408104315e18,
+];
+
+static PP: [f64; 7] = [
+    7.969367292973471e-4,
+    8.283523921074408e-2,
+    1.239533716464143,
+    5.447250030587687,
+    8.74716500199817,
+    5.303240382353949,
+    1.0,
+];
+
+static PQ: [f64; 7] = [
+    9.244088105588637e-4,
+    8.562884743544745e-2,
+    1.2535274390105895,
+    5.470977403304171,
+    8.761908832370695,
+    5.306052882353947,
+    1.0,
+];
+
+static QP: [f64; 8] = [
+    -1.1366383889846916e-2,
+    -1.2825271867050931,
+    -1.9553954425773597e1,
+    -9.320601521237683e1,
+    -1.7768116798048806e2,
+    -1.4707750515495118e2,
+    -5.141053267665993e1,
+    -6.050143506007285,
+];
+
+static QQ: [f64; 7] = [
+    6.43178256118178e1,
+    8.564300259769806e2,
+    3.8824018360540163e3,
+    7.240467741956525e3,
+    5.930727011873169e3,
+    2.0620933166032783e3,
+    2.420057402402914e2,
+];
+
+static YP: [f64; 8] = [
+    1.5592436785523574e4,
+    -1.466392959039716e7,
+    5.435264770518765e9,
+    -9.821360657179115e11,
+    8.75906394395367e13,
+    -3.466283033847297e15,
+    4.4273326857256984e16,
+    -1.8495080043698668e16,
+];
+
+static YQ: [f64; 7] = [
+    1.0412835366425984e3,
+    6.26107330137135e5,
+    2.6891963339381415e8,
+    8.64002487103935e10,
+    2.0297961275010555e13,
+    3.1715775284297505e15,
+    2.5059625617265306e17,
+];
+
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn test_j0() {
+        // reference computed with mpmath
+        assert_relative_eq!(bessel_j0(-123.985), -0.055876145864746804);
+        assert_relative_eq!(bessel_j0(-5.1), -0.14433474706050065);
+        assert_relative_eq!(bessel_j0(-3.0), -0.26005195490193345);
+        assert_relative_eq!(bessel_j0(0.0), 1.0);
+        assert_relative_eq!(bessel_j0(0.00245), 0.999998499375563);
+        assert_relative_eq!(bessel_j0(2.1752), 0.12419296628748941);
+        assert_relative_eq!(bessel_j0(2345.13), 0.012425605700760064, max_relative=1e-12);
+    }
+
+    #[test]
+    fn test_y0() {
+        // reference computed with mpmath
+        assert!(bessel_y0(-123.985).is_nan());
+        assert!(bessel_y0(-5.1).is_nan());
+        assert!(bessel_y0(-3.0).is_nan());
+        assert_relative_eq!(bessel_y0(0.0), -f64::INFINITY);
+        assert_relative_eq!(bessel_y0(0.00245), -3.90094372498198, max_relative=1e-12);
+        assert_relative_eq!(bessel_y0(2.1752), 0.520660638047155, max_relative=1e-12);
+        assert_relative_eq!(bessel_y0(2345.13), 0.0108198389384562, max_relative=1e-12);
+    }
+}
diff --git a/rascaline/src/math/bessel/j1.rs b/rascaline/src/math/bessel/j1.rs
new file mode 100644
index 000000000..e21170b40
--- /dev/null
+++ b/rascaline/src/math/bessel/j1.rs
@@ -0,0 +1,199 @@
+use crate::math::consts::SQRT_FRAC_2_PI;
+use std::f64::consts::PI;
+
+use super::{eval_polynomial, eval_polynomial_1};
+
+const Z1: f64 = 1.4681970642123893e1;
+const Z2: f64 = 4.92184563216946e1;
+
+/// Bessel function of the first kind, order one
+///
+/// The domain is divided into the intervals [0, 8] and (8, infinity). In the
+/// first interval a 24 term Chebyshev expansion is used. In the second, the
+/// asymptotic trigonometric representation is employed using two rational
+/// functions of degree 5/5.
+///
+/// # Accuracy
+/// ```text
+///                      Absolute error:
+/// arithmetic   domain      # trials      peak         rms
+///    IEEE      0, 30       30000       2.6e-16     1.1e-16
+/// ```
+#[allow(clippy::many_single_char_names)]
+pub fn bessel_j1(x: f64) -> f64 {
+    let z;
+    let mut p;
+
+    let mut w = x;
+    if x < 0.0 {
+        return -bessel_j1(-x);
+    }
+    if w <= 5.0 {
+        z = x * x;
+        w = eval_polynomial(z, &RP) / eval_polynomial_1(z, &RQ);
+        w = w * x * (z - Z1) * (z - Z2);
+        return w;
+    }
+    w = 5.0 / x;
+    z = w * w;
+    p = eval_polynomial(z, &PP) / eval_polynomial(z, &PQ);
+    let q = eval_polynomial(z, &QP) / eval_polynomial_1(z, &QQ);
+    let xn = x - 0.75 * PI;
+    p = p * f64::cos(xn) - w * q * f64::sin(xn);
+    return p * SQRT_FRAC_2_PI / f64::sqrt(x);
+}
+
+/// Bessel function of second kind, order one
+///
+/// Returns Bessel function of the second kind of order one of the argument.
+///
+/// The domain is divided into the intervals [0, 8] and (8, infinity). In the
+/// first interval a 25 term Chebyshev expansion is used, and a call to j1() is
+/// required. In the second, the asymptotic trigonometric representation is
+/// employed using two rational functions of degree 5/5.
+///
+/// # Accuracy
+///
+/// ```text
+///                      Absolute error:
+/// arithmetic   domain      # trials      peak         rms
+///    IEEE      0, 30       30000       1.0e-15     1.3e-16
+/// ```
+///
+/// (error criterion relative when |y1| > 1).
+#[allow(clippy::many_single_char_names)]
+pub fn bessel_y1(x: f64) -> f64 {
+    let mut w;
+    let z;
+    let mut p;
+    if x <= 5.0 {
+        if x == 0.0 {
+            // sf_error(SF_ERROR_SINGULAR);
+            return -std::f64::INFINITY;
+        } else if x <= 0.0 {
+            // sf_error(SF_ERROR_DOMAIN);
+            return std::f64::NAN;
+        }
+        z = x * x;
+        w = x * (eval_polynomial(z, &YP) / eval_polynomial_1(z, &YQ));
+        w += 2.0 / PI * (bessel_j1(x) * f64::ln(x) - 1.0 / x);
+        return w;
+    }
+    w = 5.0 / x;
+    z = w * w;
+    p = eval_polynomial(z, &PP) / eval_polynomial(z, &PQ);
+    let q = eval_polynomial(z, &QP) / eval_polynomial_1(z, &QQ);
+    let xn = x - 0.75 * PI;
+    p = p * f64::sin(xn) + w * q * f64::cos(xn);
+    return p * SQRT_FRAC_2_PI / f64::sqrt(x);
+}
+
+static RP: [f64; 4] = [
+    -8.999712257055594e8,
+    4.5222829799819403e11,
+    -7.274942452218183e13,
+    3.682957328638529e15,
+];
+
+static RQ: [f64; 8] = [
+    6.208364781180543e2,
+    2.5698725675774884e5,
+    8.351467914319493e7,
+    2.215115954797925e10,
+    4.749141220799914e12,
+    7.843696078762359e14,
+    8.952223361846274e16,
+    5.322786203326801e18,
+];
+
+static PP: [f64; 7] = [
+    7.621256162081731e-4,
+    7.313970569409176e-2,
+    1.1271960812968493,
+    5.112079511468076,
+    8.424045901417724,
+    5.214515986823615,
+    1.0,
+];
+
+static PQ: [f64; 7] = [
+    5.713231280725487e-4,
+    6.884559087544954e-2,
+    1.105142326340617,
+    5.073863861286015,
+    8.399855543276042,
+    5.209828486823619,
+    1.0,
+];
+
+static QP: [f64; 8] = [
+    5.108625947501766e-2,
+    4.982138729512334,
+    7.582382841325453e1,
+    3.667796093601508e2,
+    7.108563049989261e2,
+    5.974896124006136e2,
+    2.1168875710057213e2,
+    2.5207020585802372e1,
+];
+
+static QQ: [f64; 7] = [
+    7.423732770356752e1,
+    1.0564488603826283e3,
+    4.986410583376536e3,
+    9.562318924047562e3,
+    7.997041604473507e3,
+    2.8261927851763908e3,
+    3.360936078106983e2,
+];
+
+static YP: [f64; 6] = [
+    1.2632047479017804e9,
+    -6.473558763791603e11,
+    1.1450951154182373e14,
+    -8.127702555013251e15,
+    2.024394757135949e17,
+    -7.788771962659501e17,
+];
+
+static YQ: [f64; 8] = [
+    5.943015923461282e2,
+    2.3556409294306856e5,
+    7.348119444597217e7,
+    1.8760131610870617e10,
+    3.8823127749623857e12,
+    6.205577271469538e14,
+    6.871410873553005e16,
+    3.9727060811656064e18,
+];
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn test_j1() {
+        // reference computed with mpmath
+        assert_relative_eq!(bessel_j1(-123.985), 0.04508621393470377, max_relative=1e-12);
+        assert_relative_eq!(bessel_j1(-5.1), 0.3370972020182318);
+        assert_relative_eq!(bessel_j1(-3.0), -0.3390589585259365);
+        assert_relative_eq!(bessel_j1(0.0), 0.0);
+        assert_relative_eq!(bessel_j1(0.00245), 0.0012249990808674174);
+        assert_relative_eq!(bessel_j1(2.1752), 0.5593771605955342, max_relative=1e-12);
+        assert_relative_eq!(bessel_j1(2345.13), 0.010822488420270095, max_relative=1e-12);
+    }
+
+    #[test]
+    fn test_y1() {
+        // reference computed with mpmath
+        assert!(bessel_y1(-123.985).is_nan());
+        assert!(bessel_y1(-5.1).is_nan());
+        assert!(bessel_y1(-3.0).is_nan());
+        assert_relative_eq!(bessel_y1(0.0), -f64::INFINITY);
+        assert_relative_eq!(bessel_y1(0.00245), -259.84997363769, max_relative=1e-12);
+        assert_relative_eq!(bessel_y1(2.1752), -0.0114834540278188, max_relative=1e-12);
+        assert_relative_eq!(bessel_y1(2345.13), -0.0124232991092643, max_relative=1e-12);
+    }
+}
diff --git a/rascaline/src/math/bessel/jv.rs b/rascaline/src/math/bessel/jv.rs
new file mode 100644
index 000000000..b3833ce6b
--- /dev/null
+++ b/rascaline/src/math/bessel/jv.rs
@@ -0,0 +1,838 @@
+#![allow(clippy::many_single_char_names, clippy::similar_names, clippy::float_cmp)]
+#![allow(clippy::manual_range_contains, clippy::too_many_lines)]
+
+use gamma::ln_gamma;
+
+use crate::math::gamma;
+
+use super::{eval_polynomial, bessel_j0, bessel_j1, bessel_y0, bessel_y1};
+use super::{MAXLOG};
+
+const CBRT_TWO: f64 = 1.2599210498948732;
+
+
+/// Spherical Bessel function of the first kind of non-integer order
+pub fn spherical_bessel_jv(n: f64, x: f64) -> f64 {
+    return f64::sqrt(std::f64::consts::PI / (2.0 * x)) * bessel_jv(n + 0.5, x);
+}
+
+/// Spherical Bessel function of the second kind of non-integer order
+pub fn spherical_bessel_yv(n: f64, x: f64) -> f64 {
+    return f64::sqrt(std::f64::consts::PI / (2.0 * x)) * bessel_yv(n + 0.5, x);
+}
+
+/// Bessel function of second kind of integer order
+///
+/// Returns Bessel function of order n, where n is a (possibly negative)
+/// integer.
+///
+/// The function is evaluated by forward recurrence on n, starting with values
+/// computed by the routines y0() and y1().
+///
+/// If n = 0 or 1 the routine for y0 or y1 is called directly.
+///
+/// # Accuracy
+///
+/// ```text
+///                      Absolute error, except relative
+///                      when y > 1:
+/// arithmetic   domain     # trials      peak         rms
+///    IEEE      0, 30       30000       3.4e-15     4.3e-16
+/// ````
+///
+/// Spot checked against tables for x, n between 0 and 100.
+pub fn bessel_yn(mut n: i32, x: f64) -> f64 {
+    let sign;
+    if n < 0 {
+        n = -n;
+        if (n & 1) == 0 {	/* -1**n */
+            sign = 1.0;
+        } else {
+            sign = -1.0;
+        }
+    } else {
+        sign = 1.0;
+    }
+
+    if n == 0 {
+	    return sign * bessel_y0(x);
+    } else if n == 1 {
+	    return sign * bessel_y1(x);
+    }
+
+    // test for overflow
+    if x == 0.0 {
+        // sf_error(SF_ERROR_SINGULAR);
+        return -f64::INFINITY * sign;
+    } else if x < 0.0 {
+        // sf_error(SF_ERROR_DOMAIN);
+        return f64::NAN;
+    }
+
+    // forward recurrence on n
+    let mut anm2 = bessel_y0(x);
+    let mut anm1 = bessel_y1(x);
+    let mut r = 2.0;
+    let mut an = 0.0;
+    for _ in 1..n {
+        an = r * anm1 / x - anm2;
+        anm2 = anm1;
+        anm1 = an;
+        r += 2.0;
+    }
+
+    return sign * an;
+}
+
+// Bessel function of second kind, non-integer order
+pub fn bessel_yv(v: f64, x: f64) -> f64 {
+    let n = v as i32;
+    if n as f64 == v {
+        return bessel_yn(n, x);
+    } else if v == f64::floor(v) {
+        // sf_error(SF_ERROR_DOMAIN);
+        return std::f64::NAN;
+    }
+
+    let t = std::f64::consts::PI * v;
+    let y = (f64::cos(t) * bessel_jv(v, x) - bessel_jv(-v, x)) / f64::sin(t);
+    if y.is_infinite() {
+        if v > 0.0 {
+            // sf_error(SF_ERROR_OVERFLOW);
+            return -std::f64::INFINITY;
+        } else if v < -1e10 {
+            // sf_error(SF_ERROR_DOMAIN);
+            return std::f64::NAN;
+        }
+    }
+
+    return y;
+}
+
+
+
+/// Bessel function of the first kind, non-integer order
+///
+/// Returns Bessel function of order v of the argument, where v is real.
+/// Negative x is allowed if v is an integer.
+///
+/// Several expansions are included: the ascending power series, the Hankel
+/// expansion, and two transitional expansions for large v.  If v is not too
+/// large, it is reduced by recurrence to a region of best accuracy. The
+/// transitional expansions give 12D accuracy for v > 500.
+///
+/// # Accuracy
+///
+/// Results for integer v are indicated by *, where x and v both vary from -125
+/// to +125.  Otherwise, x ranges from 0 to 125, v ranges as indicated by
+/// "domain." Error criterion is absolute, except relative when |jv()| > 1.
+///
+/// ```text
+/// arithmetic  v domain  x domain    # trials      peak       rms
+///    IEEE      0,125     0,125      100000      4.6e-15    2.2e-16
+///    IEEE   -125,0       0,125       40000      5.4e-11    3.7e-13
+///    IEEE      0,500     0,500       20000      4.4e-15    4.0e-16
+/// Integer v:
+///    IEEE   -125,125   -125,125      50000      3.5e-15*   1.9e-16*
+/// ```
+pub fn bessel_jv(mut n: f64, mut x: f64) -> f64 {
+    let mut k;
+    let mut q;
+    let mut t;
+    let i;
+
+    let mut sign = 1.0;
+    let mut nint = 0;
+    let an = f64::abs(n);
+    let mut y = f64::floor(an);
+
+    if an == y {
+        nint = 1;
+        i = (an - 16384.0 * f64::floor(an / 16384.0)) as i32;
+        if n < 0.0 {
+            if i & 1 != 0 {
+                sign = -sign;
+            }
+            n = an;
+        }
+
+        if x < 0.0 {
+            if i & 1 != 0 {
+                sign = -sign;
+            }
+            x = -x;
+        }
+
+        if n == 0.0 {
+            return bessel_j0(x);
+        }
+
+        if n == 1.0 {
+            return sign * bessel_j1(x);
+        }
+    }
+
+    if x < 0.0 && y != an {
+        // sf_error(SF_ERROR_DOMAIN);
+        return std::f64::NAN;
+    }
+
+    if x == 0.0 && n < 0.0 && nint == 0 {
+        // sf_error(SF_ERROR_OVERFLOW);
+        return std::f64::INFINITY / gamma(n + 1.0);
+    }
+
+    y = f64::abs(x);
+    if y * y < f64::abs(n + 1.0) * f64::EPSILON {
+        return f64::powf(0.5 * x, n) / gamma(n + 1.0);
+    }
+
+    k = 3.6 * f64::sqrt(y);
+    t = 3.6 * f64::sqrt(an);
+
+    if y < t && an > 21.0 {
+        return sign * jvs(n, x);
+    }
+
+    if an < k && y > 21.0 {
+        return sign * hankel(n, x);
+    }
+
+    if an < 500.0 {
+    	// Note: if x is too large, the continued fraction will fail; but
+    	// then the Hankel expansion can be used.
+    	if nint != 0 {
+    	    k = 0.0;
+    	    q = recur(&mut n, x, &mut k, 1);
+
+    	    if k == 0.0 {
+    		    y = bessel_j0(x) / q;
+    		    return sign * y;
+    	    }
+
+            if k == 1.0 {
+    		    y = bessel_j1(x) / q;
+    		    return sign * y;
+    	    }
+    	}
+
+        if (an > 2.0 * y) || (n >= 0.0) && (n < 20.0) && (y > 6.0) && (y < 20.0) {
+            // Recurse backward from a larger value of n
+            k = n;
+
+            y = y + an + 1.0;
+            if y < 30.0 {
+                y = 30.0;
+            }
+
+            y = n + f64::floor(y - n);
+            q = recur(&mut y, x, &mut k, 0);
+            y = jvs(y, x) * q;
+
+            return sign * y;
+        }
+
+        if k <= 30.0 {
+            k = 2.0;
+        } else if k < 90.0 {
+            k = 3.0 * k / 4.0;
+        }
+
+        if an > (k + 3.0) {
+            if n < 0.0 {
+                k = -k;
+            }
+
+            q = n - f64::floor(n);
+            k = f64::floor(k) + q;
+            if n > 0.0 {
+                q = recur(&mut n, x, &mut k, 1);
+            } else {
+                t = k;
+                k = n;
+                q = recur(&mut t, x, &mut k, 1);
+                k = t;
+            }
+
+            if q == 0.0 {
+                y = 0.0;
+                return sign * y;
+            }
+        } else {
+            k = n;
+            q = 1.0;
+        }
+
+        // boundary between convergence of power series and Hankel expansion
+        y = f64::abs(k);
+        if y < 26.0 {
+            t = (0.0083 * y + 0.09) * y + 12.9;
+        } else {
+            t = 0.9 * y;
+        }
+
+        if x > t {
+            y = hankel(k, x);
+        } else {
+            y = jvs(k, x);
+        }
+
+        if n > 0.0 {
+            y /= q;
+        } else {
+            y *= q;
+        }
+    } else {
+        // For large n, use the uniform expansion or the transitional expansion.
+        // But if x is of the order of n**2, these may blow up, whereas the
+        // Hankel expansion will then work.
+        if n < 0.0 {
+            // sf_error(SF_ERROR_LOSS);
+            y = std::f64::NAN;
+        } else {
+            t = x / n;
+            t /= n;
+            if t > 0.3 {
+                y = hankel(n, x);
+            } else {
+                y = jnx(n, x);
+            }
+        }
+    }
+
+    return sign * y;
+}
+
+/// Reduce the order by backward recurrence.
+/// AMS55 #9.1.27 and 9.1.73.
+fn recur(n: &mut f64, x: f64, new_n: &mut f64, cancel: i32) -> f64 {
+    const BIG: f64 = 1.4411518807585587e17;
+
+    let mut pkm2;
+    let mut pkm1;
+    let mut pk;
+    let mut qkm2;
+    let mut qkm1;
+    let mut k;
+    let mut ans;
+    let mut qk;
+    let mut xk;
+    let mut yk;
+    let mut r;
+    let mut t;
+    let mut ctr;
+
+    let maxiter = 22000;
+    let mut miniter = (f64::abs(x) - f64::abs(*n)) as i32;
+    if miniter < 1 {
+        miniter = 1;
+    }
+
+    let mut nflag;
+    if *n < 0.0 {
+        nflag = 1;
+    } else {
+        nflag = 0;
+    }
+
+    loop {
+        pkm2 = 0.0;
+        qkm2 = 1.0;
+        pkm1 = x;
+        qkm1 = *n + *n;
+        xk = -x * x;
+        yk = qkm1;
+        ans = 0.0;
+        ctr = 0;
+
+        loop {
+            yk += 2.0;
+            pk = pkm1 * yk + pkm2 * xk;
+            qk = qkm1 * yk + qkm2 * xk;
+            pkm2 = pkm1;
+            pkm1 = pk;
+            qkm2 = qkm1;
+            qkm1 = qk;
+
+            if qk != 0.0 && ctr > miniter {
+                r = pk / qk;
+            } else {
+                r = 0.0;
+            }
+
+            if r == 0.0 {
+                t = 1.0;
+            } else {
+                t = f64::abs((ans - r) / r);
+                ans = r;
+            }
+
+            ctr += 1;
+            if ctr > maxiter {
+                // sf_error(SF_ERROR_UNDERFLOW);
+                break;
+            }
+
+            if t < f64::EPSILON {
+                break;
+            }
+
+            if f64::abs(pk) > BIG {
+                pkm2 /= BIG;
+                pkm1 /= BIG;
+                qkm2 /= BIG;
+                qkm1 /= BIG;
+            }
+
+            if t < f64::EPSILON {
+                break;
+            }
+        }
+        if ans == 0.0 {
+            ans = 1.0;
+        }
+        if nflag <= 0 {
+            break;
+        }
+        if f64::abs(ans) >= 0.125 {
+            break;
+        }
+        nflag = -(1);
+        *n -= 1.0;
+    }
+
+    let kf = *new_n;
+    pk = 1.0;
+    pkm1 = 1.0 / ans;
+    k = *n - 1.0;
+    r = 2.0 * k;
+    loop {
+        pkm2 = (pkm1 * r - pk * x) / x;
+        pk = pkm1;
+        pkm1 = pkm2;
+        r -= 2.0;
+        k -= 1.0;
+
+        if k <= kf + 0.5 {
+            break;
+        }
+    }
+    if cancel != 0 && kf >= 0.0 && f64::abs(pk) > f64::abs(pkm1) {
+        k += 1.0;
+        pkm2 = pk;
+    }
+    *new_n = k;
+    return pkm2;
+}
+
+// https://stackoverflow.com/a/55696477/4692076
+fn float_exponent(s: f64) -> f64 {
+    if s == 0.0 {
+        return 0.0;
+    }
+
+    let lg = s.abs().log2();
+    let exp = lg.floor() + 1.0;
+    return exp;
+}
+
+/// Ascending power series for Jv(x).
+/// AMS55 #9.1.10.
+fn jvs(n: f64, x: f64) -> f64 {
+    let z = -x * x / 4.0;
+    let mut u = 1.0;
+    let mut y = u;
+    let mut k = 1.0;
+    let mut t = 1.0;
+    while t > f64::EPSILON {
+        u *= z / (k * (n + k));
+        y += u;
+        k += 1.0;
+        if y != 0.0 {
+            t = f64::abs(u / y);
+        }
+    }
+    let mut ex = float_exponent(0.5 * x);
+    ex *= n;
+
+    if ex > -1023.0 && ex < 1023.0 && n > 0.0 && n < 171.6243769563027 - 1.0 {
+        t = f64::powf(0.5 * x, n) / gamma(n + 1.0);
+        y *= t;
+    } else {
+        dbg!("here");
+        t = n * f64::ln(0.5 * x) - ln_gamma(n + 1.0);
+        if y < 0.0 {
+            y = -y;
+        }
+        t += f64::ln(y);
+        if t < -MAXLOG {
+            return 0.0;
+        }
+        if t > MAXLOG {
+            // sf_error(SF_ERROR_OVERFLOW);
+            return std::f64::INFINITY;
+        }
+        y = f64::exp(t);
+    }
+    return y;
+}
+
+
+// Hankel's asymptotic expansion for large x.
+// AMS55 #9.2.5.
+fn hankel(n: f64, x: f64) -> f64 {
+    let m = 4.0 * n * n;
+    let mut j = 1.0;
+    let z = 8.0 * x;
+    let mut k = 1.0;
+    let mut p = 1.0;
+    let mut u = (m - 1.0) / z;
+    let mut q = u;
+    let mut sign = 1.0;
+    let mut conv = 1.0;
+    let mut flag = 0;
+    let mut t = 1.0;
+    let mut pp = 1.0e38;
+    let mut qq = 1.0e38;
+    while t > f64::EPSILON {
+        k += 2.0;
+        j += 1.0;
+        sign = -sign;
+        u *= (m - k * k) / (j * z);
+        p += sign * u;
+        k += 2.0;
+        j += 1.0;
+        u *= (m - k * k) / (j * z);
+        q += sign * u;
+        t = f64::abs(u / p);
+        if t < conv {
+            conv = t;
+            qq = q;
+            pp = p;
+            flag = 1;
+        }
+        if flag != 0 && t > conv {
+            break;
+        }
+    }
+    u = x - (0.5 * n + 0.25) * std::f64::consts::PI;
+    t = f64::sqrt(2.0 / (std::f64::consts::PI * x)) * (pp * f64::cos(u) - qq * f64::sin(u));
+    return t;
+}
+
+static LAMBDA: [f64; 11] = [
+    1.0,
+    1.0416666666666667e-1,
+    8.355034722222222e-2,
+    1.2822657455632716e-1,
+    2.9184902646414046e-1,
+    8.816272674437576e-1,
+    3.3214082818627677,
+    1.4995762986862555e1,
+    7.892301301158652e1,
+    4.744515388682643e2,
+    3.207490090890662e3,
+];
+
+static MU: [f64; 11] = [
+    1.0,
+    -1.4583333333333334e-1,
+    -9.874131944444445e-2,
+    -1.4331205391589505e-1,
+    -3.1722720267841353e-1,
+    -9.424291479571203e-1,
+    -3.5112030408263544,
+    -1.5727263620368046e1,
+    -8.228143909718594e1,
+    -4.923553705236705e2,
+    -3.3162185685479726e3,
+];
+
+static P1: [f64; 2] = [
+    -2.0833333333333334e-1,
+    1.25e-1
+];
+
+static P2: [f64; 3] = [
+    3.342013888888889e-1,
+    -4.010416666666667e-1,
+    7.03125e-2,
+];
+
+static P3: [f64; 4] = [
+    -1.0258125964506173,
+    1.8464626736111112,
+    -8.912109375e-1,
+    7.32421875e-2,
+];
+
+static P4: [f64; 5] = [
+    4.669584423426247,
+    -1.1207002616222994E1,
+    8.78912353515625,
+    -2.3640869140625,
+    1.12152099609375e-1,
+];
+
+static P5: [f64; 6] = [
+    -2.8212072558200244e1,
+    8.463621767460073e1,
+    -9.181824154324002e1,
+    4.253499874538846e1,
+    -7.368794359479632,
+    2.2710800170898438e-1,
+];
+
+static P6: [f64; 7] = [
+    2.1257013003921713e2,
+    -7.652524681411817e2,
+    1.0599904525279999e3,
+    -6.995796273761325e2,
+    2.181905117442116e2,
+    -2.6491430486951554e1,
+    5.725014209747314e-1,
+];
+
+static P7: [f64; 8] = [
+    -1.919457662318407e3,
+    8.061722181737309e3,
+    -1.3586550006434138e4,
+    1.1655393336864534e4,
+    -5.305646978613403e3,
+    1.2009029132163525e3,
+    -1.0809091978839466e2,
+    1.7277275025844574,
+];
+
+/// Asymptotic expansion for large n.
+/// AMS55 #9.3.35.
+#[allow(clippy::too_many_lines)]
+fn jnx(n: f64, x: f64) -> f64 {
+    let zeta;
+    let mut t;
+    let sz;
+    let mut sign;
+    let nflg;
+    let mut u = [0.0; 8];
+
+    let cbn = f64::cbrt(n);
+    let z = (x - n) / cbn;
+    if f64::abs(z) <= 0.7 {
+        return jnt(n, x);
+    }
+    let z = x / n;
+    let zz = 1.0 - z * z;
+
+    if zz == 0.0 {
+        return 0.0;
+    }
+
+    if zz > 0.0 {
+        sz = f64::sqrt(zz);
+        t = 1.5 * (f64::ln((1.0 + sz) / z) - sz);
+        zeta = f64::cbrt(t * t);
+        nflg = 1;
+    } else {
+        sz = f64::sqrt(-zz);
+        t = 1.5 * (sz - f64::acos(1.0 / z));
+        zeta = -f64::cbrt(t * t);
+        nflg = -(1);
+    }
+
+    let z32i = f64::abs(1.0 / t);
+    let sqz = f64::cbrt(t);
+    let n23 = f64::cbrt(n * n);
+    t = n23 * zeta;
+    let ((ai, _), (aip, _)) = super::airy(t);
+
+    u[0] = 1.0;
+
+    let zzi = 1.0 / zz;
+    u[1] = eval_polynomial(zzi, &P1) / sz;
+    u[2] = eval_polynomial(zzi, &P2) / zz;
+    u[3] = eval_polynomial(zzi, &P3) / (sz * zz);
+
+    let mut pp = zz * zz;
+    u[4] = eval_polynomial(zzi, &P4) / pp;
+    u[5] = eval_polynomial(zzi, &P5) / (pp * sz);
+
+    pp *= zz;
+    u[6] = eval_polynomial(zzi, &P6) / pp;
+    u[7] = eval_polynomial(zzi, &P7) / (pp * sz);
+
+    pp = 0.0;
+    let mut qq = 0.0;
+
+    let mut np = 1.0;
+    let mut do_a = true;
+    let mut do_b = true;
+    let mut akl = std::f64::INFINITY;
+    let mut bkl = std::f64::INFINITY;
+    let mut k = 0;
+    while k <= 3 {
+        let tk = 2 * k;
+        let tkp1 = tk + 1;
+        let mut zp = 1.0;
+        let mut ak = 0.0;
+        let mut bk = 0.0;
+        let mut s = 0;
+        while s <= tk {
+            if do_a {
+                if s & 3 > 1 {
+                    sign = nflg as f64;
+                } else {
+                    sign = 1.0;
+                }
+                ak += sign * MU[s] * zp * u[(tk - s)];
+            }
+
+            if do_b {
+                let m = tkp1 - s;
+                if (m + 1) & 3 > 1 {
+                    sign = nflg as f64;
+                } else {
+                    sign = 1.0;
+                }
+                bk += sign * LAMBDA[s] * zp * u[m];
+            }
+            zp *= z32i;
+            s += 1;
+        }
+
+        if do_a {
+            ak *= np;
+            t = f64::abs(ak);
+            if t < akl {
+                akl = t;
+                pp += ak;
+            } else {
+                do_a = false;
+            }
+        }
+
+        if do_b {
+            bk += LAMBDA[tkp1] * zp * u[0];
+            bk *= -np / sqz;
+            t = f64::abs(bk);
+            if t < bkl {
+                bkl = t;
+                qq += bk;
+            } else {
+                do_b = false;
+            }
+        }
+        if np < f64::EPSILON {
+            break;
+        }
+        np /= n * n;
+        k += 1;
+    }
+
+    t = 4.0 * zeta / zz;
+    t = f64::sqrt(f64::sqrt(t));
+    t *= ai * pp / f64::cbrt(n) + aip * qq / (n23 * n);
+    return t;
+}
+
+static PF2: [f64; 2] = [
+    -9e-2,
+    8.571428571428572e-2,
+];
+
+static PF3: [f64; 3] = [
+    1.367142857142857e-1,
+    -5.492063492063492e-2,
+    -4.4444444444444444e-3,
+];
+
+static PF4: [f64; 4] = [
+    1.35e-3,
+    -1.6036054421768708e-1,
+    4.259018759018759e-2,
+    2.733044733044733e-3,
+];
+
+static PG1: [f64; 2] = [
+    -2.4285714285714285e-1,
+    1.4285714285714285e-2,
+];
+
+static PG2: [f64; 3] = [
+    -9e-3,
+    1.9396825396825396e-1,
+    -1.1746031746031746e-2,
+];
+
+static PG3: [f64; 3] = [
+    1.9607142857142858e-2,
+    -1.5983694083694083e-1,
+    6.383838383838384e-3,
+];
+
+
+/// Asymptotic expansion for transition region, n large and x close to n.
+/// AMS55 #9.3.23.
+fn jnt(n: f64, x: f64) -> f64 {
+    let cbn = f64::cbrt(n);
+    let z = (x - n) / cbn;
+
+    let ((ai, _), (aip, _)) = super::airy(-CBRT_TWO * z);
+
+    let zz = z * z;
+    let z3 = zz * z;
+
+    let f = [
+        1.0,
+        -z / 5.0,
+        eval_polynomial(z3, &PF2) * zz,
+        eval_polynomial(z3, &PF3),
+        eval_polynomial(z3, &PF4) * z,
+    ];
+
+
+    let g = [
+        0.3 * zz,
+        eval_polynomial(z3, &PG1),
+        eval_polynomial(z3, &PG2) * z,
+        eval_polynomial(z3, &PG3) * zz,
+    ];
+
+    let mut pp = 0.0;
+    let mut qq = 0.0;
+    let mut nk = 1.0;
+    let n23 = f64::cbrt(n * n);
+
+    for k in 0..5 {
+        let fk = f[k] * nk;
+        pp += fk;
+        if k != 4 {
+            qq += g[k] * nk;
+        }
+        nk /= n23;
+    }
+
+    return CBRT_TWO * ai * pp / cbn + f64::cbrt(4.0) * aip * qq / n;
+}
+
+#[cfg(test)]
+mod tests {
+    use super::*;
+
+    use approx::assert_relative_eq;
+
+    #[test]
+    fn test_jv() {
+        todo!()
+    }
+
+    #[test]
+    fn test_yv() {
+        todo!()
+    }
+
+    #[test]
+    fn test_large_n() {
+        assert_relative_eq!(1e80 * bessel_jv(171.6, 45.4), 8.938521429926956, max_relative=1e-12);
+        assert_relative_eq!(1e80 * bessel_jv(176.0, -45.4), 0.0012279789648073855, max_relative=1e-12);
+    }
+}
diff --git a/rascaline/src/math/bessel/mod.rs b/rascaline/src/math/bessel/mod.rs
new file mode 100644
index 000000000..f71fe775c
--- /dev/null
+++ b/rascaline/src/math/bessel/mod.rs
@@ -0,0 +1,76 @@
+/*
+ * This code was ported from scipy using c2rust, followed by a manual cleanup.
+ * These files where used:
+ * - cephes/scipy_iv.c -> iv.rs
+ * - cephes/polevl.h -> polynomials.rs
+ * - cephes/j0.c -> j0.rs
+ * - cephes/j1.c -> j1.rs
+ * - cephes/{jv.c, yn.c, yv.c} -> jv.rs
+ * - cephes/airy.c -> airy.rs
+ *
+ *
+ * Parts of the code are copyright:
+ *
+ *     Cephes Math Library Release 2.8:  June, 2000
+ *     Copyright 1984, 1987, 1988, 2000 by Stephen L. Moshier
+ *
+ * And other parts:
+ *
+ *     Copyright (c) 2006 Xiaogang Zhang
+ *     Use, modification and distribution are subject to the
+ *     Boost Software License, Version 1.0.
+ *
+ *     Boost Software License - Version 1.0 - August 17th, 2003
+ *
+ *     Permission is hereby granted, free of charge, to any person or
+ *     organization obtaining a copy of the software and accompanying
+ *     documentation covered by this license (the "Software") to use, reproduce,
+ *     display, distribute, execute, and transmit the Software, and to prepare
+ *     derivative works of the Software, and to permit third-parties to whom the
+ *     Software is furnished to do so, all subject to the following:
+ *
+ *     The copyright notices in the Software and this entire statement,
+ *     including the above license grant, this restriction and the following
+ *     disclaimer, must be included in all copies of the Software, in whole or
+ *     in part, and all derivative works of the Software, unless such copies or
+ *     derivative works are solely in the form of machine-executable object code
+ *     generated by a source language processor.
+ *
+ *     THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ *     OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF
+ *     MERCHANTABILITY, FITNESS FOR A PARTICULAR PURPOSE, TITLE AND
+ *     NON-INFRINGEMENT. IN NO EVENT SHALL THE COPYRIGHT HOLDERS OR ANYONE
+ *     DISTRIBUTING THE SOFTWARE BE LIABLE FOR ANY DAMAGES OR OTHER LIABILITY,
+ *     WHETHER IN CONTRACT, TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN
+ *     CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
+ *     SOFTWARE.
+ *
+ * And the rest are:
+ *
+ *     Copyright (C) 2009 Pauli Virtanen
+ *     Distributed under the same license as Scipy.
+ *
+ */
+
+#![allow(dead_code)]
+
+const MAXLOG: f64 = 7.09782712893384e2;     /* f64::ln(f64::MAX) */
+
+mod polynomials;
+pub use self::polynomials::{eval_polynomial, eval_polynomial_1};
+
+mod airy;
+pub use self::airy::airy;
+
+mod iv;
+pub use self::iv::bessel_iv;
+
+mod j0;
+pub use self::j0::{bessel_j0, bessel_y0};
+
+mod j1;
+pub use self::j1::{bessel_j1, bessel_y1};
+
+mod jv;
+pub use self::jv::{bessel_jv, bessel_yn, bessel_yv};
+pub use self::jv::{spherical_bessel_jv, spherical_bessel_yv};
diff --git a/rascaline/src/math/bessel/polynomials.rs b/rascaline/src/math/bessel/polynomials.rs
new file mode 100644
index 000000000..0ef59cf7f
--- /dev/null
+++ b/rascaline/src/math/bessel/polynomials.rs
@@ -0,0 +1,23 @@
+
+#[inline]
+pub fn eval_polynomial(x: f64, coefficients: &[f64]) -> f64 {
+    let mut result = coefficients[0];
+
+    for c in &coefficients[1..] {
+        result = result * x + c;
+    }
+
+    return result;
+}
+
+
+/// Same as `eval_polynomial`, but setting the coefficient of x^(n + 1) to 1.
+pub fn eval_polynomial_1(x: f64, coefficients: &[f64]) -> f64 {
+    let mut result = x + coefficients[0];
+
+    for c in &coefficients[1..] {
+        result = result * x + c;
+    }
+
+    return result;
+}
diff --git a/rascaline/src/math/gamma.rs b/rascaline/src/math/gamma.rs
index d9c1cc533..1a5da1d73 100644
--- a/rascaline/src/math/gamma.rs
+++ b/rascaline/src/math/gamma.rs
@@ -324,7 +324,7 @@ pub fn digamma(x: f64) -> f64 {
 mod tests {
     use super::*;
     use approx::assert_relative_eq;
-    use crate::math::EULER;
+    use crate::math::consts::EULER;
 
     #[test]
     fn test_gamma() {
diff --git a/rascaline/src/math/hyp2f1.rs b/rascaline/src/math/hyp2f1.rs
index e99c8020b..024c3234d 100644
--- a/rascaline/src/math/hyp2f1.rs
+++ b/rascaline/src/math/hyp2f1.rs
@@ -3,7 +3,8 @@
 
 use std::f64::consts::PI;
 
-use crate::math::{digamma, gamma, EULER};
+use crate::math::{digamma, gamma};
+use crate::math::consts::EULER;
 
 fn is_integer(x: f64) -> bool {
     return (x as i32) as f64 == x;
@@ -11,12 +12,12 @@ fn is_integer(x: f64) -> bool {
 
 /// Compute the 2F1 hypergeometric function.
 ///
-/// The current implementation does not support any of the 
+/// The current implementation does not support any of the
 /// following combination of input parameters
 ///
 ///  - `x == 1` and `c - a - b < 0`
 ///  - `c < 0` and `a < 0` and `a > c`
-///  - `c < 0` and `b < 0` and `b > c` 
+///  - `c < 0` and `b < 0` and `b > c`
 ///
 /// The implementation is translated from scipy's specfun.f, and distributed
 /// under the CSD-3-Clauses license, Copyright (c) 2001-2002 Enthought, Inc.
@@ -373,7 +374,6 @@ mod tests {
         ];
 
         for (a, b, c, x, expected) in scipy_points {
-            dbg!(a, b, c, x);
             assert_relative_eq!(hyp2f1(a, b, c, x), expected, max_relative=1e-12, epsilon=1e-15);
         }
     }
diff --git a/rascaline/src/math/k_vectors.rs b/rascaline/src/math/k_vectors.rs
index 08e15c340..530acb74d 100644
--- a/rascaline/src/math/k_vectors.rs
+++ b/rascaline/src/math/k_vectors.rs
@@ -87,7 +87,7 @@ mod tests {
     use super::*;
     use crate::Matrix3;
 
-    const SQRT_3: f64 = 1.7320508075688772;
+    /// sqrt(5)
     const SQRT_5: f64 = 2.23606797749979;
 
     #[test]
@@ -117,8 +117,8 @@ mod tests {
             1.0 + eps,
             std::f64::consts::SQRT_2 - eps,
             std::f64::consts::SQRT_2 + eps,
-            SQRT_3 - eps,
-            SQRT_3 + eps,
+            crate::math::consts::SQRT_3 - eps,
+            crate::math::consts::SQRT_3 + eps,
             2.0 - eps,
             SQRT_5 - eps,
         ];
diff --git a/rascaline/src/math/mod.rs b/rascaline/src/math/mod.rs
index 20f67f725..147d017f2 100644
--- a/rascaline/src/math/mod.rs
+++ b/rascaline/src/math/mod.rs
@@ -1,5 +1,15 @@
-/// Euler's constant
-pub const EULER: f64 = 0.5772156649015329;
+mod consts {
+    /// sqrt(3)
+    pub const SQRT_3: f64 = 1.7320508075688772;
+
+    /// sqrt(2 / pi)
+    pub const SQRT_FRAC_2_PI: f64 = 0.7978845608028654;
+
+    /// Euler's constant
+    pub const EULER: f64 = 0.5772156649015329;
+}
+
+
 
 mod double_regularized_1f1;
 pub(crate) use self::double_regularized_1f1::DoubleRegularized1F1;
@@ -16,6 +26,9 @@ pub use self::erf::{erf, erfc};
 mod gamma;
 pub(crate) use self::gamma::{gamma, ln_gamma, digamma};
 
+mod bessel;
+pub use self::bessel::bessel_iv;
+
 mod hyp1f1;
 pub(crate) use self::hyp1f1::hyp1f1;
 
diff --git a/rascaline/src/math/spherical_harmonics.rs b/rascaline/src/math/spherical_harmonics.rs
index f561b0397..4e3a4d582 100644
--- a/rascaline/src/math/spherical_harmonics.rs
+++ b/rascaline/src/math/spherical_harmonics.rs
@@ -9,8 +9,6 @@ use crate::Vector3D;
 
 /// `\sqrt{\frac{1}{2 \pi}}`
 const SQRT_1_OVER_2PI: f64 = 0.3989422804014327;
-/// `\sqrt{3}`
-const SQRT_3: f64 = 1.7320508075688772;
 /// `\sqrt{3 / 2}`
 const SQRT_3_OVER_2: f64 = 1.224744871391589;
 
@@ -229,7 +227,7 @@ impl SphericalHarmonics {
         self.legendre_polynomials[[0, 0]] = value;
 
         if self.max_angular > 0 {
-            self.legendre_polynomials[[1, 0]] = cos_theta * SQRT_3 * value;
+            self.legendre_polynomials[[1, 0]] = cos_theta * crate::math::consts::SQRT_3 * value;
             value *= -SQRT_3_OVER_2 * sin_theta;
             self.legendre_polynomials[[1, 1]] = value;
 
diff --git a/rascaline/src/math/splines.rs b/rascaline/src/math/splines.rs
index d44523047..f3b99f4e1 100644
--- a/rascaline/src/math/splines.rs
+++ b/rascaline/src/math/splines.rs
@@ -290,7 +290,6 @@ mod tests {
         let mut values = ndarray::Array1::from_elem((1,), 0.0);
         let mut gradients = ndarray::Array1::from_elem((1,), 0.0);
         for &x in &[-2.2, -1.00242144, 0.0, 0.000000001, 2.3, 3.2, 4.7, 5.3, 5.99999999] {
-            dbg!(x);
             spline.compute(x, values.view_mut(), Some(gradients.view_mut()));
             assert_relative_eq!(values[0], f64::sin(x), max_relative=1e-5, epsilon=1e-12);
             assert_relative_eq!(gradients[0], f64::cos(x), max_relative=1e-5, epsilon=1e-12);