math_fun/

lib.rs

//! Special functions for science and engineering problems.
//!
//! Provides several mathematical functions that often appear in many different
//! disciplines of science and engineering.
//!
//! The goal of this package is to provide simple-to-use, pure-rust implementations
//! without many dependencies.
//! Rather than trying to exhaustively provide as many functions as possible or
//! to cover all possible argument types and ranges, implementing widely used functions
//! in an efficient way is the first priority.
//!
//! Three types of input arguments are possible for indicating the position at which
//! the function is evaluated.
//! 1. Point evaluation: a single position
//! e.g., `sph_bessel_kind1_ordern_arg_real(order: usize, x: f64)`,
//! 2. Range evaluation: Three floating point arguments, `start`, `end`, and `step`
//! to indicate the range (inclusive of `start` and exclusive of `end`).
//! `step` can be negative, in which case `end` should be less than `start`.
//! e.g., `sph_bessel_kind1_ordern_arg_real_ranged(order: usize, start: f64, end: f64, step: f64)`
//! 3. Iterator evaluation: a collection of positions that implements IntoIterator
//! e.g., `sph_bessel_kind1_ordern_arg_real_iterable(order: usize, x_list: impl IntoIterator<Item =
//! f64>`

#![warn(missing_docs)]

use std::num::FpCategory::*;

/// Accept the three arguments indicating a range, and return a vector of that range.
/// It currently panics if wrong arguments are entered.
pub fn range_to_vec(start: f64, end: f64, step: f64) -> Vec<f64> {
    let diff: f64 = end - start;

    let sds_to_vec = |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).collect::<Vec<f64>>();

    match (start.classify(), diff.classify(), step.classify()) {
        (Normal | Zero | Subnormal, Normal | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !end.is_sign_negative() &&
                diff.is_sign_negative() == step.is_sign_negative() {
                sds_to_vec(start, diff, step)
            } else {
                panic!("Both 'start' and 'end' should be non-negative, and the sign of 'step' must match the sign of 'end - start'.");
            }
        },
        (_, _, Zero) => panic!("'step' cannot be zero."),
        (_, Zero, _) => panic!("'start' and 'end' should be different."),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite."),
        _ => panic!("Improper arguments.")
    }
}

/// Spherical Bessel function of the first kind, order = 0
/// j_n(x) = \sqrt{\pi/{2x}}J_{n+0.5}(x)
pub fn sph_bessel_kind1_order0_arg_real(x: f64) -> f64 {
    match x.classify() {
        Normal | Subnormal => x.sin()/x,
        Zero => 1.0,
        Infinite => 0.0,
        Nan => f64::NAN,
    }

}

/// Spherical Bessel function of the first kind, order = 0, a range input
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind1_order0_arg_real_ranged(start: f64, end: f64, step: f64) -> Vec<f64> {
    let diff: f64 = end - start;
    let default_expr = |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize)
        .map(|n| start + (n as f64)*step).map(|x| x.sin()/x).collect::<Vec<f64>>();

    match (start.classify(), diff.classify(), step.classify()) {
        (Normal | Subnormal, Normal | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !end.is_sign_negative() &&
                diff.is_sign_negative() == step.is_sign_negative() {
                default_expr(start, diff, step)
            } else {
                panic!("Both 'start' and 'end' should be non-negative, and the sign of 'step' must match the sign of 'end - start'.");
            }
        },
        (Zero, Normal | Subnormal, Normal | Subnormal) => {
            let mut result = Vec::with_capacity(((diff/step).ceil() as usize) + 1);

            if end.is_sign_positive() && step.is_sign_positive() {
                result.push(1.0);
                result.extend_from_slice(&default_expr(start + step, diff, step));
                result
            } else {
                panic!("Both 'end' and 'step' should be positive if 'start' is zero.");
            }
        },
        (_, _, Zero) => panic!("'step' cannot be zero."),
        (_, Zero, _) => panic!("'start' and 'end' should be different."),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite"),
        _ => panic!("Improper arguments.")
    }

}

/// Spherical Bessel function of the first kind, order = 0, an iterable input
pub fn sph_bessel_kind1_order0_arg_real_iterable(x_list: impl IntoIterator<Item = f64>) -> Vec<f64> {
    x_list.into_iter().map(|x: f64| match x.classify() {
        Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { x.sin()/x },
        Zero => 1.0,
        Infinite => 0.0,
        _ => f64::NAN
    }).collect::<Vec<f64>>()
}

/// Spherical Bessel function of the first kind, order = n
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind1_ordern_arg_real(order: usize, x: f64) -> f64 {
    let default_expr = match order {
        0   => |x: f64| x.sin()/x,
        1   => |x: f64| x.sin()/x.powi(2) - x.cos()/x,
        2   => |x: f64| x.sin()*(3.0/x.powi(3)-1.0/x) - x.cos()*(3.0/x.powi(2)),
        _   => panic!("Only orders from 0 to 2 are currently implemented.")
    };

    match x.classify() {
        Normal | Subnormal => default_expr(x),
        Zero => match order {
            0 => 1.0,
            _ => 0.0
        },
        Infinite => 0.0,
        Nan => f64::NAN,
    }
}

/// Spherical Bessel function of the first kind, order = n, a range input
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind1_ordern_arg_real_ranged(order: usize, start: f64, end: f64, step: f64) -> Vec<f64> {
    // To Do: Handle the case of 'list too long' (number of points > usize.MAX)
    let diff: f64 = end - start;

    let default_expr = match order {
        0   =>  {
            |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).map(|x| x.sin()/x).collect::<Vec<f64>>()
        },
        1   =>  {
            |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).map(|x| x.sin()/x.powi(2) - x.cos()/x).collect::<Vec<f64>>()
        },
        2   =>  {
            |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).map(|x| x.sin()*(3.0/x.powi(3)-1.0/x) - x.cos()*(3.0/x.powi(2))).collect::<Vec<f64>>()
        },
        _   =>  {
            panic!("Only orders from 0 to 2 are currently implemented.");
        }
    };

    match (start.classify(), diff.classify(), step.classify()) {
        (Normal | Subnormal, Normal | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !end.is_sign_negative() && 
                diff.is_sign_negative() == step.is_sign_negative() {
                default_expr(start, diff, step)
            } else {
                panic!("Both 'start' and 'end' should be non-negative, and the sign of 'step' must match the sign of 'end - start'.");
            }
        },
        (Zero, Normal | Subnormal, Normal | Subnormal) => {
            let mut result = Vec::with_capacity(((diff/step).ceil() as usize) + 1);
                
            if end.is_sign_positive() && step.is_sign_positive() { 
                result.push(match order {
                        0       => 1.0,
                        1..=2   => 0.0,
                        _       => panic!("Only orders from 0 to 2 are currently implemented.")
                });
                result.extend_from_slice(&default_expr(start + step, diff, step));
                result
            } else {
                panic!("Both 'end' and 'step' should be positive if 'start' is zero.");
            }
        },
        (_, _, Zero) => panic!("'step' cannot be zero."),
        (_, Zero, _) => panic!("'start' and 'end' should be different."),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite."),
        _ => panic!("Improper arguments.")
    }
}

/// Spherical Bessel function of the first kind, order = n, an iterable input
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind1_ordern_arg_real_iterable(order: usize, x_list: impl IntoIterator<Item = f64>) -> Vec<f64> {
    match order {
        0   =>  {
            x_list.into_iter().map(|x| match x.classify() {
                Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { x.sin()/x },
                Zero => 1.0,
                Infinite => 0.0,
                _ => f64::NAN
            }).collect::<Vec<f64>>()
        },
        1   =>  {
            x_list.into_iter().map(|x| match x.classify() {
                Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { x.sin()/x.powi(2) - x.cos()/x },
                Zero | Infinite => 0.0,
                _ => f64::NAN
            }).collect::<Vec<f64>>()
        },
        2   =>  {
            x_list.into_iter().map(|x| match x.classify() {
                Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { x.sin()*(3.0/x.powi(3)-1.0/x) - x.cos()*(3.0/x.powi(2)) },
                Zero | Infinite => 0.0,
                _ => f64::NAN
            }).collect::<Vec<f64>>()
        },
        _   =>  {
            panic!("Only orders from 0 to 2 are currently implemented.");
        }
    }
}

/// Spherical Bessel function of the second kind, order = 0
/// j_n(x) = \sqrt{\pi/{2x}}Y_{n+0.5}(x)
pub fn sph_bessel_kind2_order0_arg_real(x: f64) -> f64 {
    match x.classify() {
        Normal | Subnormal => -x.cos()/x,
        Zero => f64::INFINITY,
        Infinite => 0.0,
        Nan => f64::NAN,
    }

}

/// Spherical Bessel function of the second kind, order = 0, a range input
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind2_order0_arg_real_ranged(start: f64, end: f64, step: f64) -> Vec<f64> {
    let diff: f64 = end - start;
    let default_expr = |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize)
        .map(|n| start + (n as f64)*step).map(|x| -x.cos()/x).collect::<Vec<f64>>();

    match (start.classify(), diff.classify(), step.classify()) {
        (Normal | Subnormal, Normal | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !end.is_sign_negative() &&
                diff.is_sign_negative() == step.is_sign_negative() {
                default_expr(start, diff, step)
            } else {
                panic!("Both 'start' and 'end' should be non-negative, and the sign of 'step' must match the sign of 'end - start'.");
            }
        },
        (Zero, Normal | Subnormal, Normal | Subnormal) => {
            let mut result = Vec::with_capacity(((diff/step).ceil() as usize) + 1);

            if end.is_sign_positive() && step.is_sign_positive() {
                result.push(f64::NEG_INFINITY);
                result.extend_from_slice(&default_expr(start + step, diff, step));
                result
            } else {
                panic!("Both 'end' and 'step' should be positive if 'start' is zero.");
            }
        },
        (_, _, Zero) => panic!("'step' cannot be zero."),
        (_, Zero, _) => panic!("'start' and 'end' should be different."),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite"),
        _ => panic!("Improper arguments.")
    }

}

/// Spherical Bessel function of the second kind, order = 0, an iterable input
pub fn sph_bessel_kind2_order0_arg_real_iterable(x_list: impl IntoIterator<Item = f64>) -> Vec<f64> {
    x_list.into_iter().map(|x: f64| match x.classify() {
        Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { -x.cos()/x },
        Zero => f64::NEG_INFINITY,
        Infinite => 0.0,
        _ => f64::NAN
    }).collect::<Vec<f64>>()
}

/// Spherical Bessel function of the second kind, order = n
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind2_ordern_arg_real(order: usize, x: f64) -> f64 {
    let default_expr = match order {
        0   => |x: f64| -x.cos()/x,
        1   => |x: f64| -x.cos()/x.powi(2) - x.sin()/x,
        2   => |x: f64| x.cos()*(-3.0/x.powi(3)+1.0/x) - x.sin()*(3.0/x.powi(2)),
        _   => panic!("Only orders from 0 to 2 are currently implemented")
    };

    match x.classify() {
        Normal | Subnormal => default_expr(x),
        Zero => f64::NEG_INFINITY,
        Infinite => 0.0,
        Nan => f64::NAN,
    }

}

/// Spherical Bessel function of the second kind, order = n, a range input
/// It currently panics if wrong arguments are entered.
/// In the future the return type may be changed to Results to deal with errors graciously.
pub fn sph_bessel_kind2_ordern_arg_real_ranged(order: usize, start: f64, end: f64, step: f64) -> Vec<f64> {
    // To Do: Handle the case of 'list too long' (number of points > usize.MAX)
    let diff: f64 = end - start;

    let default_expr = match order {
        0   =>  {
            |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).map(|x| -x.cos()/x).collect::<Vec<f64>>()
        },
        1   =>  {
            |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).map(|x| -x.cos()/x.powi(2) - x.sin()/x).collect::<Vec<f64>>()
        },
        2   =>  {
            |start: f64, diff: f64, step: f64| (0 .. (diff/step).ceil() as usize).map(|n| start + (n as f64)*step).map(|x| x.cos()*(-3.0/x.powi(3)+1.0/x) - x.sin()*(3.0/x.powi(2))).collect::<Vec<f64>>()
        },
        _   =>  {
            panic!("Only orders from 0 to 2 are currently implemented.");
        },
    };
 
    match (start.classify(), diff.classify(), step.classify()) {
        (Normal | Subnormal, Normal | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !end.is_sign_negative() && 
                diff.is_sign_negative() == step.is_sign_negative() {
                default_expr(start, diff, step)
            } else {
                panic!("Both 'start' and 'end' should be non-negative, and the sign of 'step' must match the sign of 'end - start'.");
            }
        },
        (Zero, Normal | Subnormal, Normal | Subnormal) => {
            let mut result = Vec::with_capacity(((diff/step).ceil() as usize) + 1);
                
            if end.is_sign_positive() && step.is_sign_positive() { 
                result.push(match order {
                        0..=2   => f64::NEG_INFINITY,
                        _       => panic!("Only orders from 0 to 2 are currently implemented.")
                });
                result.extend_from_slice(&default_expr(start+step, diff, step));
                result
            } else {
                panic!("Both 'end' and 'step' should be positive if 'start' is zero.");
            }
        },
        (_, _, Zero) => panic!("'step' cannot be zero."),
        (_, Zero, _) => panic!("'start' and 'end' should be different."),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite."),
        _ => panic!("Improper arguments.")
    }
}

/// Spherical Bessel function of the second kind, order = n, an iterable input
/// It currently panics if wrong arguments are entered.
pub fn sph_bessel_kind2_ordern_arg_real_iterable(order: usize, x_list: impl IntoIterator<Item = f64>) -> Vec<f64> {
    match order {
        0   =>  {
            x_list.into_iter().map(|x| match x.classify() {
                Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { -x.cos()/x },
                Zero => f64::NEG_INFINITY,
                Infinite => 0.0,
                _ => f64::NAN
            }).collect::<Vec<f64>>()
        },
        1   =>  {
            x_list.into_iter().map(|x| match x.classify() {
                Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { -x.cos()/x.powi(2) - x.sin()/x },
                Zero => f64::NEG_INFINITY,
                Infinite => 0.0,
                _ => f64::NAN
            }).collect::<Vec<f64>>()
        },
        2   =>  {
            x_list.into_iter().map(|x| match x.classify() {
                Normal | Subnormal => if x.is_sign_negative() { f64::NAN } else { x.cos()*(-3.0/x.powi(3)+1.0/x) - x.sin()*(3.0/x.powi(2)) },
                Zero => f64::NEG_INFINITY,
                Infinite => 0.0,
                _ => f64::NAN
            }).collect::<Vec<f64>>()
        },
        _   =>  {
            panic!("Only orders from 0 to 2 are currently implemented.");
        }
    }
}

#[cfg(test)]
mod tests {
    use super::*;
    use std::time::Instant;
    use plotly::{
        common::Mode,
        layout::{Axis, Layout},
        Plot, Scatter, ImageFormat
    };

    #[test]
    fn benchmark_sph_bessel_kind1_order0_1million() {
        let start_time = Instant::now();

        for n in 0..1000000 {
            // sph_bessel_1st_0_unchecked(0.5);
            // sph_bessel_1st_0_unchecked((n as f64)/1.0e5 + 0.1);
            sph_bessel_kind1_order0_arg_real((n as f64)/1.0e5 + 0.1);
        }

        let elapsed_time = start_time.elapsed();
        println!("Elapsed time: {:?}", elapsed_time);
    }

    #[test]
    fn benchmark_sph_bessel_kind1_order0_range_10million() {
        let start_time = Instant::now();

        for n in 0..1000 {
            sph_bessel_kind1_ordern_arg_real_ranged(2, (n as f64)/2.0e5 + 0.1, 10000.0, 1.0);
        }

        let elapsed_time = start_time.elapsed();
        println!("Elapsed time: {:?}", elapsed_time);
    }

    #[test]
    fn benchmark_sph_bessel_kind1_order0_iterable_10million() {
        let start_time = Instant::now();

        for n in 0..1000 {
            sph_bessel_kind1_order0_arg_real_iterable((0..10000).map(|m| ((n + m) as f64)/3.0e3)
                .collect::<Vec<f64>>());
        }

        let elapsed_time = start_time.elapsed();
        println!("Elapsed time: {:?}", elapsed_time);
    }


    #[test]
    fn plot_result() {
        let x_start: f64 = 0.0;
        let x_end: f64 = 10.0;
        let x_step: f64 = 0.05;

        let x_vec: Vec<f64> = range_to_vec(x_start, x_end, x_step);

        let sph_bessel2_0: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(0, x_start, x_end, x_step);
        let sph_bessel2_1: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(1, x_start, x_end, x_step);
        let sph_bessel2_2: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(2, x_start, x_end, x_step);

        let trace_sb2_0 = Scatter::new(x_vec.clone(), sph_bessel2_0)
            .mode(Mode::Lines)
            .name("Order 0");
        let trace_sb2_1 = Scatter::new(x_vec.clone(), sph_bessel2_1)
            .mode(Mode::Lines)
            .name("Order 1");
        let trace_sb2_2 = Scatter::new(x_vec.clone(), sph_bessel2_2)
            .mode(Mode::Lines)
            .name("Order 2");

        let mut plot_sb = Plot::new();
        plot_sb.add_trace(trace_sb2_0);
        plot_sb.add_trace(trace_sb2_1);
        plot_sb.add_trace(trace_sb2_2);

        let layout_sb = Layout::new()
            .title("Spherical Bessel functions of the second kind".into())
            .x_axis(Axis::new().range(vec![x_start, x_end]))
            .y_axis(Axis::new().range(vec![-2.5, 1.0]));
        plot_sb.set_layout(layout_sb);

        plot_sb.write_image("./sph_bessel_kind2.png", ImageFormat::PNG, 600, 400, 1.0);
    }

    #[test]
    fn test_sph_bessel_kind2() {
        println!("{:?}", sph_bessel_kind2_order0_arg_real(0.1*f64::MIN_POSITIVE));
        println!("{:?}", sph_bessel_kind2_order0_arg_real_ranged(0.0, 2.5*f64::MIN_POSITIVE, 0.5*f64::MIN_POSITIVE));
        println!("{:?}", sph_bessel_kind2_ordern_arg_real_ranged(2, 0.0, 2.5*f64::MIN_POSITIVE, 0.5*f64::MIN_POSITIVE));
    }

    #[test]
    fn test_sph_bessel_kind1_ord0_arg_iterable() {
        let x_start: f64 = 0.0;
        let x_end: f64 = 10.0;
        let x_step: f64 = 0.1;

        let x_vec: Vec<f64> = range_to_vec(x_start, x_end, x_step);

        let sb10_r: Vec<f64> = sph_bessel_kind1_order0_arg_real_ranged(x_start, x_end, x_step);
        let sb10_i: Vec<f64> = sph_bessel_kind1_order0_arg_real_iterable(x_vec.clone());

        let trace_sb10_r = Scatter::new(x_vec.clone(), sb10_r)
            .mode(Mode::Lines)
            .name("Ranged");
        let trace_sb10_i = Scatter::new(x_vec.clone(), sb10_i)
            .mode(Mode::Markers)
            .name("Iterated");

        let mut plot_sb = Plot::new();
        plot_sb.add_trace(trace_sb10_r);
        plot_sb.add_trace(trace_sb10_i);

        let layout_sb = Layout::new()
            .title("Spherical Bessel functions of the first kind".into())
            .x_axis(Axis::new().range(vec![x_start, x_end]))
            .y_axis(Axis::new().range(vec![-0.5, 1.5]));
        plot_sb.set_layout(layout_sb);

        plot_sb.write_image("./sph_bessel_kind10.png", ImageFormat::PNG, 600, 400, 1.0);
    }

    #[test]
    fn test_sph_bessel_kind1_ordn_arg_iterable() {
        let x_start: f64 = 0.0;
        let x_end: f64 = 10.0;
        let x_step: f64 = 0.1;

        let x_vec: Vec<f64> = range_to_vec(x_start, x_end, x_step);

        let sb11_r: Vec<f64> = sph_bessel_kind1_ordern_arg_real_ranged(1, x_start, x_end, x_step);
        let sb11_i: Vec<f64> = sph_bessel_kind1_ordern_arg_real_iterable(1, x_vec.clone());

        let trace_sb11_r = Scatter::new(x_vec.clone(), sb11_r)
            .mode(Mode::Lines)
            .name("Ranged");
        let trace_sb11_i = Scatter::new(x_vec.clone(), sb11_i)
            .mode(Mode::Markers)
            .name("Iterated");

        let mut plot_sb = Plot::new();
        plot_sb.add_trace(trace_sb11_r);
        plot_sb.add_trace(trace_sb11_i);

        let layout_sb = Layout::new()
            .title("Spherical Bessel functions of the first kind".into())
            .x_axis(Axis::new().range(vec![x_start, x_end]))
            .y_axis(Axis::new().range(vec![-0.5, 1.5]));
        plot_sb.set_layout(layout_sb);

        plot_sb.write_image("./sph_bessel_kind11.png", ImageFormat::PNG, 600, 400, 1.0);
    }

    #[test]
    fn test_sph_bessel_kind2_ord0_arg_iterable() {
        let x_start: f64 = 0.0;
        let x_end: f64 = 10.0;
        let x_step: f64 = 0.1;

        let x_vec: Vec<f64> = range_to_vec(x_start, x_end, x_step);

        let sb20_r: Vec<f64> = sph_bessel_kind2_order0_arg_real_ranged(x_start, x_end, x_step);
        let sb20_i: Vec<f64> = sph_bessel_kind2_order0_arg_real_iterable(x_vec.clone());

        let trace_sb20_r = Scatter::new(x_vec.clone(), sb20_r)
            .mode(Mode::Lines)
            .name("Ranged");
        let trace_sb20_i = Scatter::new(x_vec.clone(), sb20_i)
            .mode(Mode::Markers)
            .name("Iterated");

        let mut plot_sb = Plot::new();
        plot_sb.add_trace(trace_sb20_r);
        plot_sb.add_trace(trace_sb20_i);

        let layout_sb = Layout::new()
            .title("Spherical Bessel functions of the second kind".into())
            .x_axis(Axis::new().range(vec![x_start, x_end]))
            .y_axis(Axis::new().range(vec![-2.5, 1.0]));
        plot_sb.set_layout(layout_sb);

        plot_sb.write_image("./sph_bessel_kind20.png", ImageFormat::PNG, 600, 400, 1.0);
    }

    #[test]
    fn test_sph_bessel_kind2_ordn_arg_iterable() {
        let x_start: f64 = 0.0;
        let x_end: f64 = 10.0;
        let x_step: f64 = 0.1;

        let x_vec: Vec<f64> = range_to_vec(x_start, x_end, x_step);

        let sb21_r: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(1, x_start, x_end, x_step);
        let sb21_i: Vec<f64> = sph_bessel_kind2_ordern_arg_real_iterable(1, x_vec.clone());

        let trace_sb21_r = Scatter::new(x_vec.clone(), sb21_r)
            .mode(Mode::Lines)
            .name("Ranged");
        let trace_sb21_i = Scatter::new(x_vec.clone(), sb21_i)
            .mode(Mode::Markers)
            .name("Iterated");

        let mut plot_sb = Plot::new();
        plot_sb.add_trace(trace_sb21_r);
        plot_sb.add_trace(trace_sb21_i);

        let layout_sb = Layout::new()
            .title("Spherical Bessel functions of the second kind".into())
            .x_axis(Axis::new().range(vec![x_start, x_end]))
            .y_axis(Axis::new().range(vec![-2.5, 1.0]));
        plot_sb.set_layout(layout_sb);

        plot_sb.write_image("./sph_bessel_kind21.png", ImageFormat::PNG, 600, 400, 1.0);
    }

}