//! 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.
//!
//! In addition to functions that take in a single argument to evaluate the function at,
//! `sph_bessel_kind1_ordern_arg_real(order: usize, x: f64)`,
//! the same functions that take in three arguments, `start`, `stop`, and `step`,
//! to return a vector that contains the evaluation of the function over a range,
//! `sph_bessel_kind1_ordern_arg_real_ranged(order: usize, start: f64, stop: f64, step: 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.
/// In the future the return type may be changed to Results to deal with errors graciously.
pub fn range_to_vec(start: f64, stop: f64, step: f64) -> Vec<f64> {
        let default_expr = |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() as usize).map(|n| start + (n as f64)*step).collect::<Vec<f64>>();

        match (start.classify(), stop.classify(), step.classify()) {
        (Normal | Zero | Subnormal, Normal | Zero | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !stop.is_sign_negative() &&
                (stop - start).is_sign_negative() == step.is_sign_negative() {
                default_expr(start, stop, step)
            } else {
                panic!("Both 'start' and 'stop' should be non-negative, and the sign of 'step' must match the sign of 'stop - start'.");
            }
        },
        (_, _, Zero) => panic!("'step' cannot be zero"),
        (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 => x.sin()/x,
        Zero | Subnormal => 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.
/// In the future the return type may be changed to Results to deal with errors graciously.
pub fn sph_bessel_kind1_order0_arg_real_ranged(start: f64, stop: f64, step: f64) -> Vec<f64> {
    // let default_expr = |items| items.into_iter().map(|x| x.sin()/x).collect::<Vec<f64>>();
    let default_expr = |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() as usize).map(|n| start + (n as f64)*step).map(|x| x.sin()/x).collect::<Vec<f64>>();

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

            result.push(0.0);
            if stop.is_sign_positive() && step.is_sign_positive() { 
                result.extend_from_slice(&default_expr(start, stop, step));
                result
            } else {
                panic!("Both 'stop' and 'step' should be non-negative");
            }
        },
        (Zero | Subnormal, Zero | Subnormal, _) => vec![0.0],
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite"),
        _ => panic!("Improper arguments")
    }
}

/// Spherical Bessel function of the first kind, order = n
/// 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_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 => default_expr(x),
        Zero | Subnormal => 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.
/// In the future the return type may be changed to Results to deal with errors graciously.
pub fn sph_bessel_kind1_ordern_arg_real_ranged(order: usize, start: f64, stop: f64, step: f64) -> Vec<f64> {
    // To Do: Handle the case of 'list too long' (number of points > usize.MAX)

    let default_expr = match order {
        0   =>  {
            |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() as usize).map(|n| start + (n as f64)*step).map(|x| x.sin()/x).collect::<Vec<f64>>()
        },
        1   =>  {
            |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() 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, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() 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(), stop.classify(), step.classify()) {
        (Normal, Normal | Zero | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !stop.is_sign_negative() && 
                (stop - start).is_sign_negative() == step.is_sign_negative() {
                default_expr(start, stop, step)
            } else {
                panic!("Both 'start' and 'stop' should be non-negative, and the sign of 'step' must match the sign of 'stop - start'.");
            }
        },
        (Zero | Subnormal, Normal, Normal | Subnormal) => {
            let mut result = Vec::with_capacity((((stop-start)/step).floor() as usize) + 2);

            result.push(sph_bessel_kind1_ordern_arg_real(order, 0.0));
            if stop.is_sign_positive() && step.is_sign_positive() { 
                result.extend_from_slice(&default_expr(start+step, stop, step));
                result
            } else {
                panic!("Both 'stop' and 'step' should be non-negative");
            }
        },
        (Zero | Subnormal, Zero | Subnormal, _) => if order == 0 {
            vec![1.0]
        } else {
            vec![0.0] 
        },
        (_, _, Zero) => panic!("'step' cannot be zero"),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite"),
        _ => panic!("Improper arguments")
    }
}

/// Spherical Bessel function of the second kind, order = 0
/// j_n(x) = \sqrt{\pi/{2x}}J_{n+0.5}(x)
pub fn sph_bessel_kind2_order0_arg_real(x: f64) -> f64 {
    match x.classify() {
        Normal => -x.cos()/x,
        Zero | Subnormal => 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.
/// In the future the return type may be changed to Results to deal with errors graciously.
pub fn sph_bessel_kind2_order0_arg_real_ranged(start: f64, stop: f64, step: f64) -> Vec<f64> {
    // let default_expr = |items| items.into_iter().map(|x| x.sin()/x).collect::<Vec<f64>>();
    let default_expr = |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() as usize).map(|n| start + (n as f64)*step).map(|x| -x.cos()/x).collect::<Vec<f64>>();

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

            result.push(f64::INFINITY);
            if stop.is_sign_positive() && step.is_sign_positive() { 
                result.extend_from_slice(&default_expr(start, stop, step));
                result
            } else {
                panic!("Both 'stop' and 'step' should be non-negative");
            }
        },
        (Zero | Subnormal, Zero | Subnormal, _) => vec![f64::INFINITY],
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite"),
        _ => panic!("Improper arguments")
    }
}

/// Spherical Bessel function of the second kind, order = n
/// 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(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 => default_expr(x),
        Zero | Subnormal => f64::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, stop: f64, step: f64) -> Vec<f64> {
    // To Do: Handle the case of 'list too long' (number of points > usize.MAX)

    let default_expr = match order {
        0   =>  {
            |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() as usize).map(|n| start + (n as f64)*step).map(|x| -x.cos()/x).collect::<Vec<f64>>()
        },
        1   =>  {
            |start: f64, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() 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, stop: f64, step: f64| (0 .. ((stop-start)/step).floor() 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(), stop.classify(), step.classify()) {
        (Normal, Normal | Zero | Subnormal, Normal | Subnormal) => {
            if !start.is_sign_negative() && !stop.is_sign_negative() && 
                (stop - start).is_sign_negative() == step.is_sign_negative() {
                default_expr(start, stop, step)
            } else {
                panic!("Both 'start' and 'stop' should be non-negative, and the sign of 'step' must match the sign of 'stop - start'.");
            }
        },
        (Zero | Subnormal, Normal, Normal | Subnormal) => {
            let mut result = Vec::with_capacity((((stop-start)/step).floor() as usize) + 2);

            result.push(f64::INFINITY);
            if stop.is_sign_positive() && step.is_sign_positive() { 
                result.extend_from_slice(&default_expr(start+step, stop, step));
                result
            } else {
                panic!("Both 'stop' and 'step' should be non-negative");
            }
        },
        (Zero | Subnormal, Zero | Subnormal, _) => vec![f64::INFINITY],
        (_, _, Zero) => panic!("'step' cannot be zero"),
        (Infinite, _, _) | (_, Infinite, _) | (_, _, Infinite) => panic!("Arguments cannot be infinite"),
        _ => panic!("Improper arguments")
    }
}


#[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..100000 {
            // 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 plot_result() {
        let x_start: f64 = 0.0;
        let x_stop: f64 = 10.0;
        let x_step: f64 = 0.05;

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

        let sph_bessel2_0: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(0, x_start, x_stop, x_step);
        let sph_bessel2_1: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(1, x_start, x_stop, x_step);
        let sph_bessel2_2: Vec<f64> = sph_bessel_kind2_ordern_arg_real_ranged(2, x_start, x_stop, 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_stop]))
            .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));
    }
}