use ndarray::prelude::*;
use num_complex::Complex;
pub(crate) const TRI_QUAD_PTS_3: [[f64; 3]; 3] = [
[2.0 / 3.0, 1.0 / 6.0, 1.0 / 6.0],
[1.0 / 6.0, 2.0 / 3.0, 1.0 / 6.0],
[1.0 / 6.0, 1.0 / 6.0, 2.0 / 3.0],
];
pub(crate) const TRI_QUAD_WTS_3: [f64; 3] = [1.0 / 3.0, 1.0 / 3.0, 1.0 / 3.0];
pub(crate) const TET_QUAD_PTS_4: [[f64; 4]; 4] = [
[
0.5854101966249685,
0.1381966011250105,
0.1381966011250105,
0.1381966011250105,
],
[
0.1381966011250105,
0.5854101966249685,
0.1381966011250105,
0.1381966011250105,
],
[
0.1381966011250105,
0.1381966011250105,
0.5854101966249685,
0.1381966011250105,
],
[
0.1381966011250105,
0.1381966011250105,
0.1381966011250105,
0.5854101966249685,
],
];
pub(crate) const TET_QUAD_WTS_4: [f64; 4] = [0.25, 0.25, 0.25, 0.25];
pub(crate) fn bary_interp_band(bands: &[Vec<f64>], lam: &[f64], nsta: usize) -> Vec<f64> {
let mut out = vec![0.0; nsta];
for v in 0..bands.len() {
let lv = lam[v];
if lv == 0.0 {
continue;
}
for n in 0..nsta {
out[n] += bands[v][n] * lv;
}
}
out
}
pub(crate) fn bary_interp_matrix(
mats: &[Array2<Complex<f64>>],
lam: &[f64],
) -> Array2<Complex<f64>> {
let n = mats[0].nrows();
let mut out = Array2::<Complex<f64>>::zeros((n, n));
for (mat, &w) in mats.iter().zip(lam.iter()) {
if w == 0.0 {
continue;
}
for i in 0..n {
for j in 0..n {
out[[i, j]] += mat[[i, j]] * w;
}
}
}
out
}