use nalgebra::{*, allocator::Allocator};
use std::f64::consts::{PI, E};

/* dynamic matrices */
pub fn rand_eigval_series(dim: usize, time_res: usize) -> Vec<OVector<Complex<f64>, Dyn>> {
    // initialize the random matrix
    let mut rand_mat = DMatrix::<f64>::from_fn(dim, dim, |j, k| {
        let n = j*dim + k;
        E*((n*n) as f64) % 2.0 - 1.0
    }) * (3.0 / (dim as f64)).sqrt();
    
    // initialize the rotation step
    let mut rot_step = DMatrix::<f64>::identity(dim, dim);
    let max_freq = 4;
    for n in (0..dim).step_by(2) {
        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigval_series = Vec::<OVector<Complex<f64>, Dyn>>::with_capacity(time_res);
    eigval_series.push(rand_mat.complex_eigenvalues());
    for _ in 1..time_res {
        rand_mat = &rot_step * rand_mat;
        eigval_series.push(rand_mat.complex_eigenvalues());
    }
    eigval_series
}

/* dynamic single float matrices */
/*pub fn rand_eigval_series(dim: usize, time_res: usize) -> Vec<OVector<Complex<f32>, Dyn>> {
    // initialize the random matrix
    let mut rand_mat = DMatrix::<f32>::from_fn(dim, dim, |j, k| {
        let n = j*dim + k;
        (E as f32)*((n*n) as f32) % 2.0_f32 - 1.0_f32
    }) * (3.0_f32 / (dim as f32)).sqrt();
    
    // initialize the rotation step
    let mut rot_step = DMatrix::<f32>::identity(dim, dim);
    let max_freq = 4;
    for n in (0..dim).step_by(2) {
        let ang = (PI as f32) * ((n % max_freq) as f32) / (time_res as f32);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigval_series = Vec::<OVector<Complex<f32>, Dyn>>::with_capacity(time_res);
    eigval_series.push(rand_mat.complex_eigenvalues());
    for _ in 1..time_res {
        rand_mat = &rot_step * rand_mat;
        eigval_series.push(rand_mat.complex_eigenvalues());
    }
    eigval_series
}*/

/* static matrices. should only be used when the dimension is really small */
/*pub fn rand_eigval_series<N>(time_res: usize) -> Vec<OVector<Complex<f64>, N>>
    where
        N: ToTypenum + DimName + DimSub<U1>,
        DefaultAllocator:
            Allocator<N> +
            Allocator<N, N> +
            Allocator<<N as DimSub<U1>>::Output> +
            Allocator<N, <N as DimSub<U1>>::Output>
{
    // initialize the random matrix
    let dim = N::try_to_usize().unwrap();
    let mut rand_mat = OMatrix::<f64, N, N>::from_fn(|j, k| {
        let n = j*dim + k;
        E*((n*n) as f64) % 2.0 - 1.0
    }) * (3.0 / (dim as f64)).sqrt();
    
    // initialize the rotation step
    let mut rot_step = OMatrix::<f64, N, N>::identity();
    let max_freq = 4;
    for n in (0..dim).step_by(2) {
        let ang = PI * ((n % max_freq) as f64) / (time_res as f64);
        let ang_cos = ang.cos();
        let ang_sin = ang.sin();
        rot_step[(n, n)] = ang_cos;
        rot_step[(n+1, n)] = ang_sin;
        rot_step[(n, n+1)] = -ang_sin;
        rot_step[(n+1, n+1)] = ang_cos;
    }
    
    // find the eigenvalues
    let mut eigval_series = Vec::<OVector<Complex<f64>, N>>::with_capacity(time_res);
    eigval_series.push(rand_mat.complex_eigenvalues());
    for _ in 1..time_res {
        rand_mat = &rot_step * rand_mat;
        eigval_series.push(rand_mat.complex_eigenvalues());
    }
    eigval_series
}*/