use crate::MU0_4PI;
use crate::math::{cross, distance, dot3, mag3, unit_vector};
use crate::sources::element::edge_integral;
// transform a global xyz position into a local coordinate system
fn transform(g: &[f64; 3], xhat: &[f64; 3], yhat: &[f64; 3], zhat: &[f64; 3]) -> [f64; 3] {
    [dot3(xhat, g), dot3(yhat, g), dot3(zhat, g)]
}

// Compute the unit vectors of an edge csys given 3 nodes
// Nodes should be in right hand rule order
// middle node is the 'anchor', xhat points from node 1 to node 0
// yhat points away from the center of the face
fn edge_csys(x: &[f64; 3], y: &[f64; 3], z: &[f64; 3]) -> ([f64; 3], [f64; 3], [f64; 3]) {
    // first edge defines x'-axis
    let xp_hat = unit_vector(&[x[1], y[1], z[1]], &[x[0], y[0], z[0]]);

    // second edge defines the next in plane direction to find the out of plane z
    // it is not returned because it is not guaranteed to be orthogonal to x and z
    let yp_hat_temp = unit_vector(&[x[2], y[2], z[2]], &[x[1], y[1], z[1]]);

    // unit vector normal to surface
    let mut zp_hat = [0.0; 3];
    cross(&xp_hat, &yp_hat_temp, &mut zp_hat);
    let norm = mag3(zp_hat[0], zp_hat[1], zp_hat[2]);
    zp_hat[0] /= norm;
    zp_hat[1] /= norm;
    zp_hat[2] /= norm;

    // other in-plane direction, not necessarily aligned with an edge
    let mut yp_hat = [0.0; 3];
    cross(&xp_hat, &zp_hat, &mut yp_hat);
    (xp_hat, yp_hat, zp_hat)
}

/// Compute the magnetic flux density generated by an arbitrary hexahedron in 3D space
///
/// # Arguments
/// - `x`: x coordinates of each corner node in 3D space
/// - `y`: y coordinates of each corner node in 3D space
/// - `z`: z coordinates of each corner node in 3D space
/// - `jdensity`: current density vector, assumed constant throughout the element
///
/// # Returns
///     magnetic flux density vector generated by the hexahedron
///
/// # Reference  
///
/// CRYO-06-034: "MAGNUM - A Fortran Library for the Calculation of Magnetic Configurations"
/// <https://supermagnet.sourceforge.io/notes/CRYO-06-034.pdf>
///
pub fn bfield_hexahedron(
    nx: &[f64],       // 8-length
    ny: &[f64],       // 8-length
    nz: &[f64],       // 8-length
    jdensity: &[f64], // 3-length
    target: &[f64],   // 3-length
) -> [f64; 3] {
    // Compute a total of 24 line integrals, along each edge (4 on each face)
    // Refer to the gmsh or ansys documentation for node numbering on a hex element:
    // <https://gmsh.info/doc/texinfo/gmsh.html#Node-ordering>
    // <https://ansyshelp.ansys.com/public/account/secured?returnurl=//////////Views/Secured/corp/v242/en/ans_elem/Hlp_E_SOLID236.html>
    let node_map: [[[usize; 3]; 4]; 6] = [
        [[0, 3, 2], [3, 2, 1], [2, 1, 0], [1, 0, 3]], // bottom face
        [[4, 5, 6], [5, 6, 7], [6, 7, 4], [7, 4, 5]], // top face
        [[0, 4, 7], [4, 7, 3], [7, 3, 0], [3, 0, 4]], // left face
        [[1, 2, 6], [2, 6, 5], [6, 5, 1], [5, 1, 2]], // right face
        [[0, 1, 5], [1, 5, 4], [5, 4, 0], [4, 0, 1]], // front face
        [[2, 3, 7], [3, 7, 6], [7, 6, 2], [6, 2, 3]], // back face
    ];

    let mut s = [0.0; 3];
    // outer loop over faces
    for i in 0..6 {
        let mut ni = [0.0; 3]; // normal on face 
        let mut s_face = 0.0;

        for j in 0..4 {
            // Node numbers from the map
            let n0 = node_map[i][j][0];
            let n1 = node_map[i][j][1];
            let n2 = node_map[i][j][2];

            // Local coordinate system tensor; transform the target point into the local system
            let (xhat, yhat, zhat) = edge_csys(
                &[nx[n0], nx[n1], nx[n2]],
                &[ny[n0], ny[n1], ny[n2]],
                &[nz[n0], nz[n1], nz[n2]],
            );
            ni = zhat;
            let target_local = [target[0] - nx[n1], target[1] - ny[n1], target[2] - nz[n1]];
            let p = transform(&target_local, &xhat, &yhat, &zhat);
            let xppq2 = distance(&[nx[n1], ny[n1], nz[n1]], &[nx[n0], ny[n0], nz[n0]]); // length of the edge, basically

            s_face +=
                p[1] * (edge_integral(xppq2 - p[0], p[1], p[2]) - edge_integral(-p[0], p[1], p[2]));
        }
        s[0] += ni[0] * s_face;
        s[1] += ni[1] * s_face;
        s[2] += ni[2] * s_face;
    }

    let mut b = [0.0; 3];
    cross(&[jdensity[0], jdensity[1], jdensity[2]], &s, &mut b);

    [MU0_4PI * b[0], MU0_4PI * b[1], MU0_4PI * b[2]]
}

#[cfg(test)]
mod tests {
    use super::*;

    #[test]
    fn test_hex() {
        // unit cube with origin at center
        let nx = [0.5, 0.5, -0.5, -0.5, 0.5, 0.5, -0.5, -0.5];
        let ny = [-0.5, 0.5, 0.5, -0.5, -0.5, 0.5, 0.5, -0.5];
        let nz = [-0.5, -0.5, -0.5, -0.5, 0.5, 0.5, 0.5, 0.5];

        let target = [0.5, 0.55, -0.55];
        let jdensity = [0.0, 1000.0e5, 0.0];

        let b_hex = bfield_hexahedron(&nx, &ny, &nz, &jdensity, &target);
        let mut jxr = [0.0; 3];
        cross(&jdensity, &target, &mut jxr);
        let multiplier = MU0_4PI * 1.0 / distance(&target, &[0.0; 3]).powi(3);
        let b_point = [
            multiplier * jxr[0],
            multiplier * jxr[1],
            multiplier * jxr[2],
        ];

        println!("Point: {:?}", b_point);
        println!("Hex:   {:?}", b_hex);
    }
}
