lie 0.1.2

A numerical library for working with Lie Groups and Algebras
Documentation
use ndarray as nd;

use num_complex::Complex64;

use crate::su2;

pub fn c_g_p(r: i32, u: i32) -> f64 {
    let r = r as f64;
    let u = u as f64;
    let x: f64 = (((1. + r + u) * (2. + r + u)) / ((1. + r) * (1. + 2. * r))).sqrt();
    x / 2f64.sqrt()
}

pub fn c_g_z(r: i32, u: i32) -> f64 {
    let r = r as f64;
    let u = u as f64;
    (((1. + r - u) * (1. + r + u)) / ((1. + r) * (1. + 2. * r))).sqrt()
}

pub fn c_g_n(r: i32, u: i32) -> f64 {
    let r = r as f64;
    let u = u as f64;
    f64::sqrt(
        (2. + 3. * r + f64::powf(r, 2.) - 3. * u - 2. * r * u + f64::powf(u, 2.))
            / (2. + 6. * r + 4. * f64::powf(r, 2.)),
    )
}

pub fn c_g(r: i32, u: i32, la: i32) -> f64 {
    match la {
        1 => c_g_p(r, u),
        0 => c_g_z(r, u),
        -1 => c_g_n(r, u),
        _ => panic!("Input for Clebsh Gordan doesn't make sense!"),
    }
}

pub fn q_1_u(j: f64, u: i32) -> nd::Array2<f64> {
    use su2::{s_x, s_y, s_z};
    match u {
        0 => s_z(j),
        1 => -(s_x(j) + s_y(j)),
        -1 => (s_x(j) - s_y(j)),
        _ => panic!("Bad input for u"),
    }
}

pub fn q_r_u(j: f64, r: i32, u: i32) -> nd::Array2<f64> {
    if r == 1 {
        q_1_u(j, u)
    } else {
        let n = (j * 2. + 1.) as usize;

        let mut mat: nd::Array2<f64> = nd::Array2::zeros((n, n));

        for i in (-1 + u)..=(u + 1) {
            if i.abs() < r {
                let q1u = q_1_u(j, -i + u);
                mat = mat + c_g(r - 1, i, -i + u) * q_r_u(j, r - 1, i).dot(&q1u);
            }
        }
        mat
    }
}

pub fn basis_from_spin(j: f64) -> Vec<nd::Array2<f64>> {
    let n = (j * 2. + 1.) as u32;

    let mut basis = Vec::with_capacity(2usize.pow(n) - 1);

    for r in 1..(n as i32) {
        for u in -r..=r {
            let x = q_r_u(j, r, u);
            basis.push(x);
        }
    }
    basis
}

pub fn hermitian_basis_from_spin(j: f64) -> Vec<nd::Array2<Complex64>> {
    let n = (j * 2. + 1.) as u32;
    let basis = basis_from_spin(j);

    let mut herm_basis: Vec<nd::Array2<Complex64>> = Vec::with_capacity(2usize.pow(n) - 1);

    let mut i = 0;

    for r in 1..n {
        for u in 0..=r {
            let j: u32 = i + r + u;
            let mut sign = 1;
            if u % 2 == 1 {
                sign = -1;
            }
            let idx = j - 2 * u;
            let b = (sign as f64) * &basis[idx as usize];
            let u_a = &basis[j as usize];

            let u1 = (u_a + &b).map(|&x| Complex64::new(x, 0.));

            herm_basis.push(u1 / 2.);

            if u == 0 {
                continue;
            }

            let u2 = u_a - &b;
            herm_basis.push(u2.map(|&x| Complex64::new(0., x)) / 2.);
        }
        i += 2 * r + 1;
    }
    herm_basis
}