integrate 0.3.1

// extern crate test;
use std::f64::consts::PI;

use num::{ToPrimitive, Unsigned};

/// Returns the $k$-th positive zero of the Bessel function $J_0(x)$.
///
/// For $k \le 20$ the zero is read from a precomputed table of high-precision
/// values.  For $k > 20$ the zero is approximated via McMahon's asymptotic
/// expansion (see [`mcmahon_expansion`]).
///
/// This function is used internally to seed the node computation for
/// Gauss-Legendre quadrature rules of large order.
pub fn bessel_j0_zeros<U: Unsigned + ToPrimitive>(k: U) -> f64 {
    const J_Z: [f64; 20] = [
        2.404_825_557_695_773,
        5.520_078_110_286_311,
        8.653_727_912_911_013,
        1.179_153_443_901_428_1E1,
        1.493_091_770_848_778_7E1,
        1.807_106_396_791_092_4E1,
        2.121_163_662_987_926E1,
        2.435_247_153_074_930_2E1,
        2.749_347_913_204_025_3E1,
        3.063_460_646_843_197_6E1,
        3.377_582_021_357_357E1,
        3.691_709_835_366_404_5E1,
        4.005_842_576_462_824E1,
        4.319_979_171_317_673E1,
        4.634_118_837_166_181_5E1,
        4.948_260_989_739_781_5E1,
        5.262_405_184_111_5E1,
        5.576_551_075_501_998E1,
        5.890_698_392_608_094E1,
        6.204_846_919_022_716_6E1,
    ];

    let r: f64;
    let r2: f64;
    let mut z: f64;

    if J_Z.len() < k.to_usize().unwrap() {
        z = PI * (k.to_f64().unwrap() - 0.25E+00);
        r = 1.0E+00 / z;
        r2 = r * r;
        z = mcmahon_expansion(z, r, r2);
    } else {
        z = J_Z[k.to_usize().unwrap() - 1];
    }

    z
}

/// Applies McMahon's asymptotic expansion to refine the $k$-th zero of $J_0(x)$.
///
/// The starting approximation is $z = \pi(k - \tfrac{1}{4})$; the expansion
/// adds correction terms as a power series in $r = 1/z$ and $r^2$, improving
/// accuracy for large $k$ (i.e. large zeros).
///
/// Intended for internal use within the Gauss-Legendre node computation.
fn mcmahon_expansion(z: f64, r: f64, r2: f64) -> f64 {
    z + r
        * (0.125E+00
            + r2 * (-8.072_916_666_666_667E-2
                + r2 * (2.460_286_458_333_333_4E-1
                    + r2 * (-1.824_438_767_206_101
                        + r2 * (2.533_641_479_734_390_6E1
                            + r2 * (-5.676_444_121_351_834E2
                                + r2 * (1.869_047_652_823_206_6E4
                                    + r2 * (-8.493_535_802_991_488E5
                                        + r2 * 5.092_254_624_022_268E7))))))))
}

/// Applies Hankel's asymptotic expansion to approximate $J_1(x_k)^2$, where
/// $x_k$ is the $k$-th zero of $J_0$.
///
/// The argument `x` is $1/(k - 1/4)$ (a small parameter for large $k$) and
/// `x2 = x * x`.  The result is the squared value of $J_1$ evaluated at that
/// zero, used as a weight denominator in Gauss-Legendre quadrature.
///
/// Intended for internal use within the Gauss-Legendre node computation.
fn hankel_expansion(x: f64, x2: f64) -> f64 {
    x * (2.026_423_672_846_755_5E-1
        + x2 * x2
            * (-3.033_804_297_112_902_7E-4
                + x2 * (1.989_243_642_459_693E-4
                    + x2 * (-2.289_699_027_721_116_6E-4
                        + x2 * (4.337_107_191_307_463E-4
                            + x2 * (-1.236_323_497_271_754E-3
                                + x2 * (4.961_014_232_688_831_4E-3
                                    + x2 * (-2.668_373_937_023_237_7E-2
                                        + x2 * 1.853_953_982_063_456_2E-1))))))))
}

/// Apply $J_1$ Bessel function to the $k$ zeros of the $J_0(x)$ Bessel function.       
/// Then computes their squared value.
///
/// # Notes
///
/// Note that the first 21 first values are tabulated.  After that, they are computed using Hankel expansion.
pub fn bessel_j1_squared<U: Unsigned + ToPrimitive>(k: U) -> f64 {
    const J_1: &[f64; 21] = &[
        2.695_141_239_419_169E-1,
        1.157_801_385_822_036_9E-1,
        7.368_635_113_640_822E-2,
        5.403_757_319_811_628E-2,
        4.266_142_901_724_309E-2,
        3.524_210_349_099_61E-2,
        3.002_107_010_305_467_3E-2,
        2.614_739_149_530_809E-2,
        2.315_912_182_469_139_3E-2,
        2.078_382_912_226_786E-2,
        1.885_045_066_931_767E-2,
        1.724_615_756_966_500_8E-2,
        1.589_351_810_592_36E-2,
        1.473_762_609_647_219E-2,
        1.373_846_514_538_711_7E-2,
        1.286_618_173_761_513_3E-2,
        1.209_805_154_862_679_7E-2,
        1.141_647_122_449_160_9E-2,
        1.080_759_279_118_020_4E-2,
        1.026_037_292_628_076_2E-2,
        9.765_897_139_791_051E-3,
    ];

    let x: f64;
    let x2: f64;
    let z: f64;

    if J_1.len() < k.to_usize().unwrap() {
        x = 1.0 / (k.to_f64().unwrap() - 0.25);
        x2 = x * x;

        z = hankel_expansion(x, x2);
    } else {
        z = J_1[k.to_usize().unwrap() - 1];
    }
    z
}

#[cfg(test)]
mod tests {

    use num_traits::Pow;

    use super::*;

    const EPSILON: f64 = 10e-5;

    // Computed using scipy special function jn_zeros
    // https://docs.scipy.org/doc/scipy/reference/generated/scipy.special.jn_zeros.html
    const BESSEL_J0_FIRST_1000_ZEROS: [f64; 1_000] = [
        2.40482556e+00,
        5.52007811e+00,
        8.65372791e+00,
        1.17915344e+01,
        1.49309177e+01,
        1.80710640e+01,
        2.12116366e+01,
        2.43524715e+01,
        2.74934791e+01,
        3.06346065e+01,
        3.37758202e+01,
        3.69170984e+01,
        4.00584258e+01,
        4.31997917e+01,
        4.63411884e+01,
        4.94826099e+01,
        5.26240518e+01,
        5.57655108e+01,
        5.89069839e+01,
        6.20484692e+01,
        6.51899648e+01,
        6.83314693e+01,
        7.14729816e+01,
        7.46145006e+01,
        7.77560256e+01,
        8.08975559e+01,
        8.40390908e+01,
        8.71806298e+01,
        9.03221726e+01,
        9.34637188e+01,
        9.66052680e+01,
        9.97468199e+01,
        1.02888374e+02,
        1.06029931e+02,
        1.09171490e+02,
        1.12313050e+02,
        1.15454613e+02,
        1.18596177e+02,
        1.21737742e+02,
        1.24879309e+02,
        1.28020877e+02,
        1.31162446e+02,
        1.34304017e+02,
        1.37445588e+02,
        1.40587160e+02,
        1.43728734e+02,
        1.46870308e+02,
        1.50011882e+02,
        1.53153458e+02,
        1.56295034e+02,
        1.59436611e+02,
        1.62578189e+02,
        1.65719767e+02,
        1.68861345e+02,
        1.72002925e+02,
        1.75144504e+02,
        1.78286084e+02,
        1.81427665e+02,
        1.84569246e+02,
        1.87710827e+02,
        1.90852409e+02,
        1.93993991e+02,
        1.97135573e+02,
        2.00277156e+02,
        2.03418739e+02,
        2.06560322e+02,
        2.09701906e+02,
        2.12843490e+02,
        2.15985074e+02,
        2.19126658e+02,
        2.22268243e+02,
        2.25409827e+02,
        2.28551412e+02,
        2.31692998e+02,
        2.34834583e+02,
        2.37976169e+02,
        2.41117755e+02,
        2.44259341e+02,
        2.47400927e+02,
        2.50542513e+02,
        2.53684100e+02,
        2.56825686e+02,
        2.59967273e+02,
        2.63108860e+02,
        2.66250447e+02,
        2.69392034e+02,
        2.72533621e+02,
        2.75675209e+02,
        2.78816796e+02,
        2.81958384e+02,
        2.85099972e+02,
        2.88241560e+02,
        2.91383148e+02,
        2.94524736e+02,
        2.97666324e+02,
        3.00807912e+02,
        3.03949500e+02,
        3.07091089e+02,
        3.10232677e+02,
        3.13374266e+02,
        3.16515855e+02,
        3.19657444e+02,
        3.22799032e+02,
        3.25940621e+02,
        3.29082210e+02,
        3.32223799e+02,
        3.35365388e+02,
        3.38506978e+02,
        3.41648567e+02,
        3.44790156e+02,
        3.47931746e+02,
        3.51073335e+02,
        3.54214925e+02,
        3.57356514e+02,
        3.60498104e+02,
        3.63639693e+02,
        3.66781283e+02,
        3.69922873e+02,
        3.73064463e+02,
        3.76206053e+02,
        3.79347642e+02,
        3.82489232e+02,
        3.85630822e+02,
        3.88772412e+02,
        3.91914002e+02,
        3.95055593e+02,
        3.98197183e+02,
        4.01338773e+02,
        4.04480363e+02,
        4.07621953e+02,
        4.10763544e+02,
        4.13905134e+02,
        4.17046724e+02,
        4.20188315e+02,
        4.23329905e+02,
        4.26471496e+02,
        4.29613086e+02,
        4.32754677e+02,
        4.35896267e+02,
        4.39037858e+02,
        4.42179449e+02,
        4.45321039e+02,
        4.48462630e+02,
        4.51604221e+02,
        4.54745811e+02,
        4.57887402e+02,
        4.61028993e+02,
        4.64170584e+02,
        4.67312175e+02,
        4.70453766e+02,
        4.73595356e+02,
        4.76736947e+02,
        4.79878538e+02,
        4.83020129e+02,
        4.86161720e+02,
        4.89303311e+02,
        4.92444902e+02,
        4.95586493e+02,
        4.98728084e+02,
        5.01869675e+02,
        5.05011267e+02,
        5.08152858e+02,
        5.11294449e+02,
        5.14436040e+02,
        5.17577631e+02,
        5.20719222e+02,
        5.23860814e+02,
        5.27002405e+02,
        5.30143996e+02,
        5.33285587e+02,
        5.36427179e+02,
        5.39568770e+02,
        5.42710361e+02,
        5.45851953e+02,
        5.48993544e+02,
        5.52135135e+02,
        5.55276727e+02,
        5.58418318e+02,
        5.61559909e+02,
        5.64701501e+02,
        5.67843092e+02,
        5.70984684e+02,
        5.74126275e+02,
        5.77267867e+02,
        5.80409458e+02,
        5.83551050e+02,
        5.86692641e+02,
        5.89834233e+02,
        5.92975824e+02,
        5.96117416e+02,
        5.99259007e+02,
        6.02400599e+02,
        6.05542190e+02,
        6.08683782e+02,
        6.11825374e+02,
        6.14966965e+02,
        6.18108557e+02,
        6.21250148e+02,
        6.24391740e+02,
        6.27533332e+02,
        6.30674923e+02,
        6.33816515e+02,
        6.36958107e+02,
        6.40099698e+02,
        6.43241290e+02,
        6.46382882e+02,
        6.49524474e+02,
        6.52666065e+02,
        6.55807657e+02,
        6.58949249e+02,
        6.62090841e+02,
        6.65232432e+02,
        6.68374024e+02,
        6.71515616e+02,
        6.74657208e+02,
        6.77798799e+02,
        6.80940391e+02,
        6.84081983e+02,
        6.87223575e+02,
        6.90365167e+02,
        6.93506759e+02,
        6.96648350e+02,
        6.99789942e+02,
        7.02931534e+02,
        7.06073126e+02,
        7.09214718e+02,
        7.12356310e+02,
        7.15497902e+02,
        7.18639493e+02,
        7.21781085e+02,
        7.24922677e+02,
        7.28064269e+02,
        7.31205861e+02,
        7.34347453e+02,
        7.37489045e+02,
        7.40630637e+02,
        7.43772229e+02,
        7.46913821e+02,
        7.50055413e+02,
        7.53197005e+02,
        7.56338597e+02,
        7.59480189e+02,
        7.62621781e+02,
        7.65763373e+02,
        7.68904965e+02,
        7.72046557e+02,
        7.75188149e+02,
        7.78329741e+02,
        7.81471333e+02,
        7.84612925e+02,
        7.87754517e+02,
        7.90896109e+02,
        7.94037701e+02,
        7.97179293e+02,
        8.00320885e+02,
        8.03462477e+02,
        8.06604069e+02,
        8.09745661e+02,
        8.12887253e+02,
        8.16028845e+02,
        8.19170437e+02,
        8.22312029e+02,
        8.25453621e+02,
        8.28595213e+02,
        8.31736805e+02,
        8.34878397e+02,
        8.38019990e+02,
        8.41161582e+02,
        8.44303174e+02,
        8.47444766e+02,
        8.50586358e+02,
        8.53727950e+02,
        8.56869542e+02,
        8.60011134e+02,
        8.63152726e+02,
        8.66294319e+02,
        8.69435911e+02,
        8.72577503e+02,
        8.75719095e+02,
        8.78860687e+02,
        8.82002279e+02,
        8.85143871e+02,
        8.88285464e+02,
        8.91427056e+02,
        8.94568648e+02,
        8.97710240e+02,
        9.00851832e+02,
        9.03993424e+02,
        9.07135017e+02,
        9.10276609e+02,
        9.13418201e+02,
        9.16559793e+02,
        9.19701385e+02,
        9.22842977e+02,
        9.25984570e+02,
        9.29126162e+02,
        9.32267754e+02,
        9.35409346e+02,
        9.38550938e+02,
        9.41692531e+02,
        9.44834123e+02,
        9.47975715e+02,
        9.51117307e+02,
        9.54258900e+02,
        9.57400492e+02,
        9.60542084e+02,
        9.63683676e+02,
        9.66825268e+02,
        9.69966861e+02,
        9.73108453e+02,
        9.76250045e+02,
        9.79391637e+02,
        9.82533230e+02,
        9.85674822e+02,
        9.88816414e+02,
        9.91958006e+02,
        9.95099599e+02,
        9.98241191e+02,
        1.00138278e+03,
        1.00452438e+03,
        1.00766597e+03,
        1.01080756e+03,
        1.01394915e+03,
        1.01709074e+03,
        1.02023234e+03,
        1.02337393e+03,
        1.02651552e+03,
        1.02965711e+03,
        1.03279871e+03,
        1.03594030e+03,
        1.03908189e+03,
        1.04222348e+03,
        1.04536508e+03,
        1.04850667e+03,
        1.05164826e+03,
        1.05478985e+03,
        1.05793144e+03,
        1.06107304e+03,
        1.06421463e+03,
        1.06735622e+03,
        1.07049781e+03,
        1.07363941e+03,
        1.07678100e+03,
        1.07992259e+03,
        1.08306418e+03,
        1.08620578e+03,
        1.08934737e+03,
        1.09248896e+03,
        1.09563055e+03,
        1.09877214e+03,
        1.10191374e+03,
        1.10505533e+03,
        1.10819692e+03,
        1.11133851e+03,
        1.11448011e+03,
        1.11762170e+03,
        1.12076329e+03,
        1.12390488e+03,
        1.12704648e+03,
        1.13018807e+03,
        1.13332966e+03,
        1.13647125e+03,
        1.13961284e+03,
        1.14275444e+03,
        1.14589603e+03,
        1.14903762e+03,
        1.15217921e+03,
        1.15532081e+03,
        1.15846240e+03,
        1.16160399e+03,
        1.16474558e+03,
        1.16788718e+03,
        1.17102877e+03,
        1.17417036e+03,
        1.17731195e+03,
        1.18045355e+03,
        1.18359514e+03,
        1.18673673e+03,
        1.18987832e+03,
        1.19301991e+03,
        1.19616151e+03,
        1.19930310e+03,
        1.20244469e+03,
        1.20558628e+03,
        1.20872788e+03,
        1.21186947e+03,
        1.21501106e+03,
        1.21815265e+03,
        1.22129425e+03,
        1.22443584e+03,
        1.22757743e+03,
        1.23071902e+03,
        1.23386062e+03,
        1.23700221e+03,
        1.24014380e+03,
        1.24328539e+03,
        1.24642699e+03,
        1.24956858e+03,
        1.25271017e+03,
        1.25585176e+03,
        1.25899336e+03,
        1.26213495e+03,
        1.26527654e+03,
        1.26841813e+03,
        1.27155972e+03,
        1.27470132e+03,
        1.27784291e+03,
        1.28098450e+03,
        1.28412609e+03,
        1.28726769e+03,
        1.29040928e+03,
        1.29355087e+03,
        1.29669246e+03,
        1.29983406e+03,
        1.30297565e+03,
        1.30611724e+03,
        1.30925883e+03,
        1.31240043e+03,
        1.31554202e+03,
        1.31868361e+03,
        1.32182520e+03,
        1.32496680e+03,
        1.32810839e+03,
        1.33124998e+03,
        1.33439157e+03,
        1.33753317e+03,
        1.34067476e+03,
        1.34381635e+03,
        1.34695794e+03,
        1.35009954e+03,
        1.35324113e+03,
        1.35638272e+03,
        1.35952431e+03,
        1.36266591e+03,
        1.36580750e+03,
        1.36894909e+03,
        1.37209068e+03,
        1.37523228e+03,
        1.37837387e+03,
        1.38151546e+03,
        1.38465705e+03,
        1.38779864e+03,
        1.39094024e+03,
        1.39408183e+03,
        1.39722342e+03,
        1.40036501e+03,
        1.40350661e+03,
        1.40664820e+03,
        1.40978979e+03,
        1.41293138e+03,
        1.41607298e+03,
        1.41921457e+03,
        1.42235616e+03,
        1.42549775e+03,
        1.42863935e+03,
        1.43178094e+03,
        1.43492253e+03,
        1.43806412e+03,
        1.44120572e+03,
        1.44434731e+03,
        1.44748890e+03,
        1.45063049e+03,
        1.45377209e+03,
        1.45691368e+03,
        1.46005527e+03,
        1.46319686e+03,
        1.46633846e+03,
        1.46948005e+03,
        1.47262164e+03,
        1.47576323e+03,
        1.47890483e+03,
        1.48204642e+03,
        1.48518801e+03,
        1.48832960e+03,
        1.49147120e+03,
        1.49461279e+03,
        1.49775438e+03,
        1.50089597e+03,
        1.50403757e+03,
        1.50717916e+03,
        1.51032075e+03,
        1.51346234e+03,
        1.51660394e+03,
        1.51974553e+03,
        1.52288712e+03,
        1.52602871e+03,
        1.52917031e+03,
        1.53231190e+03,
        1.53545349e+03,
        1.53859508e+03,
        1.54173668e+03,
        1.54487827e+03,
        1.54801986e+03,
        1.55116145e+03,
        1.55430305e+03,
        1.55744464e+03,
        1.56058623e+03,
        1.56372782e+03,
        1.56686942e+03,
        1.57001101e+03,
        1.57315260e+03,
        1.57629419e+03,
        1.57943579e+03,
        1.58257738e+03,
        1.58571897e+03,
        1.58886056e+03,
        1.59200216e+03,
        1.59514375e+03,
        1.59828534e+03,
        1.60142693e+03,
        1.60456853e+03,
        1.60771012e+03,
        1.61085171e+03,
        1.61399330e+03,
        1.61713490e+03,
        1.62027649e+03,
        1.62341808e+03,
        1.62655967e+03,
        1.62970127e+03,
        1.63284286e+03,
        1.63598445e+03,
        1.63912604e+03,
        1.64226764e+03,
        1.64540923e+03,
        1.64855082e+03,
        1.65169241e+03,
        1.65483401e+03,
        1.65797560e+03,
        1.66111719e+03,
        1.66425878e+03,
        1.66740038e+03,
        1.67054197e+03,
        1.67368356e+03,
        1.67682515e+03,
        1.67996675e+03,
        1.68310834e+03,
        1.68624993e+03,
        1.68939152e+03,
        1.69253312e+03,
        1.69567471e+03,
        1.69881630e+03,
        1.70195789e+03,
        1.70509949e+03,
        1.70824108e+03,
        1.71138267e+03,
        1.71452426e+03,
        1.71766586e+03,
        1.72080745e+03,
        1.72394904e+03,
        1.72709063e+03,
        1.73023223e+03,
        1.73337382e+03,
        1.73651541e+03,
        1.73965700e+03,
        1.74279860e+03,
        1.74594019e+03,
        1.74908178e+03,
        1.75222337e+03,
        1.75536497e+03,
        1.75850656e+03,
        1.76164815e+03,
        1.76478974e+03,
        1.76793134e+03,
        1.77107293e+03,
        1.77421452e+03,
        1.77735611e+03,
        1.78049771e+03,
        1.78363930e+03,
        1.78678089e+03,
        1.78992248e+03,
        1.79306408e+03,
        1.79620567e+03,
        1.79934726e+03,
        1.80248885e+03,
        1.80563045e+03,
        1.80877204e+03,
        1.81191363e+03,
        1.81505522e+03,
        1.81819682e+03,
        1.82133841e+03,
        1.82448000e+03,
        1.82762159e+03,
        1.83076319e+03,
        1.83390478e+03,
        1.83704637e+03,
        1.84018796e+03,
        1.84332956e+03,
        1.84647115e+03,
        1.84961274e+03,
        1.85275433e+03,
        1.85589593e+03,
        1.85903752e+03,
        1.86217911e+03,
        1.86532071e+03,
        1.86846230e+03,
        1.87160389e+03,
        1.87474548e+03,
        1.87788708e+03,
        1.88102867e+03,
        1.88417026e+03,
        1.88731185e+03,
        1.89045345e+03,
        1.89359504e+03,
        1.89673663e+03,
        1.89987822e+03,
        1.90301982e+03,
        1.90616141e+03,
        1.90930300e+03,
        1.91244459e+03,
        1.91558619e+03,
        1.91872778e+03,
        1.92186937e+03,
        1.92501096e+03,
        1.92815256e+03,
        1.93129415e+03,
        1.93443574e+03,
        1.93757733e+03,
        1.94071893e+03,
        1.94386052e+03,
        1.94700211e+03,
        1.95014370e+03,
        1.95328530e+03,
        1.95642689e+03,
        1.95956848e+03,
        1.96271007e+03,
        1.96585167e+03,
        1.96899326e+03,
        1.97213485e+03,
        1.97527644e+03,
        1.97841804e+03,
        1.98155963e+03,
        1.98470122e+03,
        1.98784281e+03,
        1.99098441e+03,
        1.99412600e+03,
        1.99726759e+03,
        2.00040918e+03,
        2.00355078e+03,
        2.00669237e+03,
        2.00983396e+03,
        2.01297555e+03,
        2.01611715e+03,
        2.01925874e+03,
        2.02240033e+03,
        2.02554193e+03,
        2.02868352e+03,
        2.03182511e+03,
        2.03496670e+03,
        2.03810830e+03,
        2.04124989e+03,
        2.04439148e+03,
        2.04753307e+03,
        2.05067467e+03,
        2.05381626e+03,
        2.05695785e+03,
        2.06009944e+03,
        2.06324104e+03,
        2.06638263e+03,
        2.06952422e+03,
        2.07266581e+03,
        2.07580741e+03,
        2.07894900e+03,
        2.08209059e+03,
        2.08523218e+03,
        2.08837378e+03,
        2.09151537e+03,
        2.09465696e+03,
        2.09779855e+03,
        2.10094015e+03,
        2.10408174e+03,
        2.10722333e+03,
        2.11036492e+03,
        2.11350652e+03,
        2.11664811e+03,
        2.11978970e+03,
        2.12293129e+03,
        2.12607289e+03,
        2.12921448e+03,
        2.13235607e+03,
        2.13549766e+03,
        2.13863926e+03,
        2.14178085e+03,
        2.14492244e+03,
        2.14806404e+03,
        2.15120563e+03,
        2.15434722e+03,
        2.15748881e+03,
        2.16063041e+03,
        2.16377200e+03,
        2.16691359e+03,
        2.17005518e+03,
        2.17319678e+03,
        2.17633837e+03,
        2.17947996e+03,
        2.18262155e+03,
        2.18576315e+03,
        2.18890474e+03,
        2.19204633e+03,
        2.19518792e+03,
        2.19832952e+03,
        2.20147111e+03,
        2.20461270e+03,
        2.20775429e+03,
        2.21089589e+03,
        2.21403748e+03,
        2.21717907e+03,
        2.22032066e+03,
        2.22346226e+03,
        2.22660385e+03,
        2.22974544e+03,
        2.23288703e+03,
        2.23602863e+03,
        2.23917022e+03,
        2.24231181e+03,
        2.24545340e+03,
        2.24859500e+03,
        2.25173659e+03,
        2.25487818e+03,
        2.25801978e+03,
        2.26116137e+03,
        2.26430296e+03,
        2.26744455e+03,
        2.27058615e+03,
        2.27372774e+03,
        2.27686933e+03,
        2.28001092e+03,
        2.28315252e+03,
        2.28629411e+03,
        2.28943570e+03,
        2.29257729e+03,
        2.29571889e+03,
        2.29886048e+03,
        2.30200207e+03,
        2.30514366e+03,
        2.30828526e+03,
        2.31142685e+03,
        2.31456844e+03,
        2.31771003e+03,
        2.32085163e+03,
        2.32399322e+03,
        2.32713481e+03,
        2.33027640e+03,
        2.33341800e+03,
        2.33655959e+03,
        2.33970118e+03,
        2.34284277e+03,
        2.34598437e+03,
        2.34912596e+03,
        2.35226755e+03,
        2.35540915e+03,
        2.35855074e+03,
        2.36169233e+03,
        2.36483392e+03,
        2.36797552e+03,
        2.37111711e+03,
        2.37425870e+03,
        2.37740029e+03,
        2.38054189e+03,
        2.38368348e+03,
        2.38682507e+03,
        2.38996666e+03,
        2.39310826e+03,
        2.39624985e+03,
        2.39939144e+03,
        2.40253303e+03,
        2.40567463e+03,
        2.40881622e+03,
        2.41195781e+03,
        2.41509940e+03,
        2.41824100e+03,
        2.42138259e+03,
        2.42452418e+03,
        2.42766577e+03,
        2.43080737e+03,
        2.43394896e+03,
        2.43709055e+03,
        2.44023214e+03,
        2.44337374e+03,
        2.44651533e+03,
        2.44965692e+03,
        2.45279852e+03,
        2.45594011e+03,
        2.45908170e+03,
        2.46222329e+03,
        2.46536489e+03,
        2.46850648e+03,
        2.47164807e+03,
        2.47478966e+03,
        2.47793126e+03,
        2.48107285e+03,
        2.48421444e+03,
        2.48735603e+03,
        2.49049763e+03,
        2.49363922e+03,
        2.49678081e+03,
        2.49992240e+03,
        2.50306400e+03,
        2.50620559e+03,
        2.50934718e+03,
        2.51248877e+03,
        2.51563037e+03,
        2.51877196e+03,
        2.52191355e+03,
        2.52505514e+03,
        2.52819674e+03,
        2.53133833e+03,
        2.53447992e+03,
        2.53762152e+03,
        2.54076311e+03,
        2.54390470e+03,
        2.54704629e+03,
        2.55018789e+03,
        2.55332948e+03,
        2.55647107e+03,
        2.55961266e+03,
        2.56275426e+03,
        2.56589585e+03,
        2.56903744e+03,
        2.57217903e+03,
        2.57532063e+03,
        2.57846222e+03,
        2.58160381e+03,
        2.58474540e+03,
        2.58788700e+03,
        2.59102859e+03,
        2.59417018e+03,
        2.59731177e+03,
        2.60045337e+03,
        2.60359496e+03,
        2.60673655e+03,
        2.60987814e+03,
        2.61301974e+03,
        2.61616133e+03,
        2.61930292e+03,
        2.62244452e+03,
        2.62558611e+03,
        2.62872770e+03,
        2.63186929e+03,
        2.63501089e+03,
        2.63815248e+03,
        2.64129407e+03,
        2.64443566e+03,
        2.64757726e+03,
        2.65071885e+03,
        2.65386044e+03,
        2.65700203e+03,
        2.66014363e+03,
        2.66328522e+03,
        2.66642681e+03,
        2.66956840e+03,
        2.67271000e+03,
        2.67585159e+03,
        2.67899318e+03,
        2.68213477e+03,
        2.68527637e+03,
        2.68841796e+03,
        2.69155955e+03,
        2.69470115e+03,
        2.69784274e+03,
        2.70098433e+03,
        2.70412592e+03,
        2.70726752e+03,
        2.71040911e+03,
        2.71355070e+03,
        2.71669229e+03,
        2.71983389e+03,
        2.72297548e+03,
        2.72611707e+03,
        2.72925866e+03,
        2.73240026e+03,
        2.73554185e+03,
        2.73868344e+03,
        2.74182503e+03,
        2.74496663e+03,
        2.74810822e+03,
        2.75124981e+03,
        2.75439140e+03,
        2.75753300e+03,
        2.76067459e+03,
        2.76381618e+03,
        2.76695777e+03,
        2.77009937e+03,
        2.77324096e+03,
        2.77638255e+03,
        2.77952415e+03,
        2.78266574e+03,
        2.78580733e+03,
        2.78894892e+03,
        2.79209052e+03,
        2.79523211e+03,
        2.79837370e+03,
        2.80151529e+03,
        2.80465689e+03,
        2.80779848e+03,
        2.81094007e+03,
        2.81408166e+03,
        2.81722326e+03,
        2.82036485e+03,
        2.82350644e+03,
        2.82664803e+03,
        2.82978963e+03,
        2.83293122e+03,
        2.83607281e+03,
        2.83921440e+03,
        2.84235600e+03,
        2.84549759e+03,
        2.84863918e+03,
        2.85178078e+03,
        2.85492237e+03,
        2.85806396e+03,
        2.86120555e+03,
        2.86434715e+03,
        2.86748874e+03,
        2.87063033e+03,
        2.87377192e+03,
        2.87691352e+03,
        2.88005511e+03,
        2.88319670e+03,
        2.88633829e+03,
        2.88947989e+03,
        2.89262148e+03,
        2.89576307e+03,
        2.89890466e+03,
        2.90204626e+03,
        2.90518785e+03,
        2.90832944e+03,
        2.91147103e+03,
        2.91461263e+03,
        2.91775422e+03,
        2.92089581e+03,
        2.92403741e+03,
        2.92717900e+03,
        2.93032059e+03,
        2.93346218e+03,
        2.93660378e+03,
        2.93974537e+03,
        2.94288696e+03,
        2.94602855e+03,
        2.94917015e+03,
        2.95231174e+03,
        2.95545333e+03,
        2.95859492e+03,
        2.96173652e+03,
        2.96487811e+03,
        2.96801970e+03,
        2.97116129e+03,
        2.97430289e+03,
        2.97744448e+03,
        2.98058607e+03,
        2.98372766e+03,
        2.98686926e+03,
        2.99001085e+03,
        2.99315244e+03,
        2.99629404e+03,
        2.99943563e+03,
        3.00257722e+03,
        3.00571881e+03,
        3.00886041e+03,
        3.01200200e+03,
        3.01514359e+03,
        3.01828518e+03,
        3.02142678e+03,
        3.02456837e+03,
        3.02770996e+03,
        3.03085155e+03,
        3.03399315e+03,
        3.03713474e+03,
        3.04027633e+03,
        3.04341792e+03,
        3.04655952e+03,
        3.04970111e+03,
        3.05284270e+03,
        3.05598429e+03,
        3.05912589e+03,
        3.06226748e+03,
        3.06540907e+03,
        3.06855067e+03,
        3.07169226e+03,
        3.07483385e+03,
        3.07797544e+03,
        3.08111704e+03,
        3.08425863e+03,
        3.08740022e+03,
        3.09054181e+03,
        3.09368341e+03,
        3.09682500e+03,
        3.09996659e+03,
        3.10310818e+03,
        3.10624978e+03,
        3.10939137e+03,
        3.11253296e+03,
        3.11567455e+03,
        3.11881615e+03,
        3.12195774e+03,
        3.12509933e+03,
        3.12824092e+03,
        3.13138252e+03,
        3.13452411e+03,
        3.13766570e+03,
        3.14080730e+03,
    ];
    // Compute Bessel J1 of first 100 roots of Bessel J0
    // Computed using scipy special function j1
    const BESSEL_J1_FIRST_100: [f64; 100] = [
        0.5191474972894669,
        -0.34026480655836816,
        0.27145229992838193,
        -0.23245983136472478,
        0.20654643307799597,
        -0.18772880304043946,
        0.17326589422922986,
        -0.16170155068925002,
        0.15218121377059457,
        -0.14416597768637324,
        0.13729694340850299,
        -0.13132462666866793,
        0.12606949712727336,
        -0.12139862477175016,
        0.11721119889066538,
        -0.1134291926164298,
        0.10999114304627802,
        -0.10684788825471286,
        0.10395957286936207,
        -0.10129349893394328,
        0.09882255380119992,
        -0.09652404046467991,
        0.09437879398467641,
        -0.09237050482355331,
        0.09048519416295768,
        -0.0887108024409698,
        0.0870368633240976,
        -0.08545424291091484,
        0.08395492928345757,
        -0.0825318613083098,
        0.08117878831953207,
        -0.07989015430874276,
        0.07866100171930493,
        -0.07748689103965987,
        0.07636383321829138,
        -0.07528823255205501,
        0.0742568381822715,
        -0.07326670270620797,
        0.07231514670236977,
        -0.071399728196232,
        0.0705182162733357,
        -0.06966856819003345,
        0.06884890944684943,
        -0.068057516381685,
        0.06729280091473991,
        -0.06655329713771117,
        0.06583764948942708,
        -0.065144602300793,
        0.06447299052550669,
        -0.06382173150081485,
        0.063189817605714,
        -0.06257630970331138,
        0.06198033127024816,
        -0.06140106312970013,
        0.06083773871596136,
        -0.060289639808346895,
        0.059756092680416394,
        -0.059236464617564454,
        0.0587301607620443,
        -0.058236621249651344,
        0.057755318606728306,
        -0.05728575537997613,
        0.05682746197485656,
        -0.05637999468123287,
        0.05594293386737833,
        -0.05551588232564262,
        0.05509846375495184,
        -0.05469032136696416,
        0.0542911166041467,
        -0.05390052795930568,
        0.05351824988721487,
        -0.0531439917999698,
        0.05277747713856061,
        -0.05241844251392233,
        0.052066636911401065,
        -0.05172182095317582,
        0.051383766213711116,
        -0.05105225458379304,
        0.05072707767912493,
        -0.050408036289839954,
        0.05009493986762599,
        -0.0497876060474634,
        0.04948586020124814,
        -0.04918953502081814,
        0.048898470128121006,
        -0.0486125117104598,
        0.04833151217893186,
        -0.04805532984833819,
        0.04778382863698613,
        -0.047516877784940494,
        0.04725435158939828,
        -0.046996129155970165,
        0.04674209416475137,
        -0.04649213465015328,
        0.046246142793549716,
        -0.04600401472786514,
        0.04576565035330103,
        -0.04553095316345705,
        0.04529983008116179,
        -0.04507219130337834,
    ];

    #[test]
    fn test_bessel_j0() {
        for k in 1..=1_000_usize {
            // computing root using method cited in Ignace Bogaert paper
            let x: f64 = bessel_j0_zeros::<usize>(k);

            // reference value
            let j0x: f64 = BESSEL_J0_FIRST_1000_ZEROS[k - 1];

            println!(
                "bessel_j0_zero: {} \t scipy: {} \t difference: {}",
                x,
                j0x,
                (x - j0x).abs()
            );

            assert!(x - j0x < EPSILON);
        }
    }

    #[test]
    fn test_bessel_j1_squared() {
        for k in 1..100_usize {
            // computing bessel J0 root using method cited in Ignace Bogaert paper
            let x: f64 = bessel_j1_squared(k);

            // reference value
            let j1 = BESSEL_J1_FIRST_100[k - 1];
            let j1_squared: f64 = j1.pow(2);

            if k < 10 {
                println!(
                    "bessel_j1_squared: {} \t scipy: {} \t difference: {}",
                    x,
                    j1_squared,
                    (x - j1_squared).abs()
                );
            }

            assert!(x - j1_squared < EPSILON);
        }
    }
}