pub mod color_gradients;
pub mod configuration;
pub mod pole_figures;
use crate::configuration::{
config::Config, opt::Opt, particle_record::ParticleRecord, record::Record,
};
use crate::pole_figures::make_pole_figures::*;
use crate::pole_figures::{
crystal_axis::CrystalAxes, lambert::*, minerals::Mineral, pole_figure::PoleFigure,
};
use ndarray::{Array, Array2, Axis};
use rayon::prelude::*;
use structopt::StructOpt;
use std::fs;
use std::fs::File;
use std::io::BufRead;
use std::io::BufReader;
use std::io::Read;
use std::path::Path;
use std::path::PathBuf;
use std::time::Instant;
pub fn run() -> Result<(), Box<dyn std::error::Error>> {
let opt = Opt::from_args();
let config_file = opt.config_file;
let configuration =
load_configuration_file(config_file).expect("Could not load configuration file");
process_configuration(configuration)
}
pub fn load_configuration_file(config_file: PathBuf) -> Result<Config, toml::de::Error> {
let config_file_display = config_file.display();
let mut file = match File::open(&config_file) {
Err(why) => panic!("couldn't open {}: {}", config_file_display, why.to_string()),
Ok(file) => file,
};
let mut config_file_string = String::new();
match file.read_to_string(&mut config_file_string) {
Err(why) => panic!("couldn't read {}: {}", config_file_display, why.to_string()),
Ok(_) => (),
}
toml::from_str(&config_file_string)
}
pub fn process_configuration(config: Config) -> Result<(), Box<dyn std::error::Error>> {
let before = Instant::now();
let base_dir = config.base_dir.clone();
let experiment_dirs = config.experiment_dirs.clone();
experiment_dirs.par_iter().for_each(|experiment_dir| {
println!("Processing experiment {}", experiment_dir);
let lpo_dir = base_dir.clone() + &experiment_dir;
let statistics_file =
lpo_dir.to_owned() + &config.pole_figures.as_ref().unwrap().time_data_file;
println!("time data file:{}", statistics_file);
let file = File::open(statistics_file).unwrap();
let reader = BufReader::new(file);
let mut data: String = "".to_string();
for line in reader.lines() {
let line = line.unwrap();
let line = line.trim();
let mut line = line.replace(" ", " ");
while let Some(_) = line.find(" ") {
line = line.replace(" ", " ");
}
if line.find('#') != Some(0) {
if line.find("particle_LPO") != None {
data = data + &line + "\n";
}
}
}
let mut timestep_to_time: Vec<f64> = vec![];
let mut rdr = csv::ReaderBuilder::new()
.trim(csv::Trim::All)
.delimiter(b' ')
.comment(Some(b'#'))
.has_headers(false)
.from_reader(data.as_bytes());
for result in rdr.records() {
let record = result.unwrap().clone();
let time = record.get(1).clone();
match time {
Some(time) => timestep_to_time.push(time.to_string().parse::<f64>().unwrap()),
None => assert!(false, "Time not found"),
}
}
if config.pole_figures.is_some() {
let pole_figure_configuration = config.pole_figures.as_ref().unwrap();
let elastisity_header = pole_figure_configuration.elastisity_header;
for output_time in &pole_figure_configuration.times {
let after_time = timestep_to_time.iter().position(|x| x > &output_time);
let after_timestep = match after_time {
Some(timestep) => timestep,
None => timestep_to_time.len() - 1,
};
let mut before_timestep = after_timestep;
if after_timestep > 0 {
before_timestep = after_timestep - 1
}
let before_timestep_diff = (output_time - timestep_to_time[before_timestep]).abs();
let after_timestep_diff = (output_time - timestep_to_time[after_timestep]).abs();
let time_step = if before_timestep_diff < after_timestep_diff {
before_timestep as u64
} else {
after_timestep as u64
};
let time = timestep_to_time[time_step as usize];
println!(
"Processing time {} (requested time: {}), located in timestep : {}",
time, output_time, time_step,
);
fs::create_dir_all(
lpo_dir.to_owned() + &pole_figure_configuration.figure_output_dir,
)
.unwrap();
let file_prefix_figures = pole_figure_configuration.figure_output_dir.to_owned()
+ &pole_figure_configuration.figure_output_prefix;
let mut rank_id = 0;
let gam = 1.0;
let sphere_points = 301;
println!(
"particle ids size {}",
pole_figure_configuration.particle_ids.len()
);
for particle_id in &pole_figure_configuration.particle_ids {
println!("processing particle_id {}", particle_id);
let mut particle_olivine_a_axis_vectors = Vec::new();
let mut particle_olivine_b_axis_vectors = Vec::new();
let mut particle_olivine_c_axis_vectors = Vec::new();
let mut particle_enstatite_a_axis_vectors = Vec::new();
let mut particle_enstatite_b_axis_vectors = Vec::new();
let mut particle_enstatite_c_axis_vectors = Vec::new();
let mut file_found: bool = false;
while !file_found {
let angles_file = format!(
"{}{}-{:05}.{:04}.dat",
lpo_dir,
pole_figure_configuration.grain_data_file_prefix,
time_step,
rank_id
);
let angles_file = Path::new(&angles_file);
let mut config_mineral_string = String::new();
for mineral in pole_figure_configuration.minerals.clone() {
config_mineral_string = format!(
"{}{}",
config_mineral_string.clone(),
match mineral {
Mineral::Olivine => {
"oli_"
}
Mineral::Enstatite => {
"ens_"
}
}
)
}
let mut config_axis_string = String::new();
for axis in pole_figure_configuration.axes.clone() {
config_axis_string = format!(
"{}{}",
config_axis_string.clone(),
match axis {
CrystalAxes::AAxis => {
"A-"
}
CrystalAxes::BAxis => {
"B-"
}
CrystalAxes::CAxis => {
"C-"
}
}
)
}
config_axis_string = format!("{}Axis_", config_axis_string.clone());
let output_file = format!(
"{}{}_{}{}{}{}_g{}_sp{}_t{:05}.{:05}.png",
lpo_dir,
file_prefix_figures,
if elastisity_header {
"elastic_"
} else {
"no-elastic_"
},
config_mineral_string,
config_axis_string,
pole_figure_configuration.color_scale.to_string(),
gam,
sphere_points,
time_step,
particle_id
);
let output_file = Path::new(&output_file);
let particle_file = format!(
"{}{}-{:05}.{:04}.dat",
lpo_dir,
pole_figure_configuration.particle_data_file_prefix,
time_step,
rank_id
);
let particle_info_file = Path::new(&particle_file);
println!(" trying file name: {}", angles_file.display());
if !(fs::metadata(angles_file).is_ok()) {
println!(
"particle id {} not found for timestep {}.",
particle_id, time_step
);
break;
}
if fs::metadata(angles_file).unwrap().len() == 0 {
rank_id = rank_id + 1;
continue;
}
println!(
" found particle id {} in:{}",
particle_id,
angles_file.display()
);
let file = File::open(angles_file).unwrap();
let metadata = file.metadata().unwrap();
let mut buf_reader =
BufReader::with_capacity(metadata.len() as usize, file);
let mut decoded_data = Vec::new();
let compressed = config.compressed;
let decoded_reader = if compressed {
let mut decoder = libflate::zlib::Decoder::new(buf_reader).unwrap();
decoder.read_to_end(&mut decoded_data).unwrap();
String::from_utf8_lossy(&decoded_data)
} else {
let data = buf_reader.fill_buf().unwrap();
String::from_utf8_lossy(&data)
};
let mut rdr = csv::ReaderBuilder::new()
.has_headers(true)
.delimiter(b' ')
.from_reader(decoded_reader.as_bytes());
let mut integer = 0;
for result in rdr.deserialize() {
let record: Record = result.unwrap();
if record.id == *particle_id {
let deg_to_rad = std::f64::consts::PI / 180.;
let euler_angles = Array::from(vec![
record.mineral_0_EA_phi.unwrap() * deg_to_rad,
record.mineral_0_EA_theta.unwrap() * deg_to_rad,
record.mineral_0_EA_z.unwrap() * deg_to_rad,
]);
let rotation_matrix =
euler_angles_to_rotation_matrix(euler_angles).unwrap();
particle_olivine_a_axis_vectors
.push(rotation_matrix.row(0).to_owned());
particle_olivine_b_axis_vectors
.push(rotation_matrix.row(1).to_owned());
particle_olivine_c_axis_vectors
.push(rotation_matrix.row(2).to_owned());
let euler_angles = Array::from(vec![
record.mineral_1_EA_phi.unwrap() * deg_to_rad,
record.mineral_1_EA_theta.unwrap() * deg_to_rad,
record.mineral_1_EA_z.unwrap() * deg_to_rad,
]);
let rotation_matrix =
euler_angles_to_rotation_matrix(euler_angles).unwrap();
particle_enstatite_a_axis_vectors
.push(rotation_matrix.row(0).to_owned());
particle_enstatite_b_axis_vectors
.push(rotation_matrix.row(1).to_owned());
particle_enstatite_c_axis_vectors
.push(rotation_matrix.row(2).to_owned());
}
integer = integer + 1;
}
if particle_olivine_a_axis_vectors.len() == 0 {
rank_id = rank_id + 1;
continue;
}
file_found = true;
let mut particle_record = ParticleRecord {
id: 0,
x: 0.0,
y: 0.0,
z: Some(0.0),
olivine_deformation_type: None,
full_norm_square: None,
triclinic_norm_square_p1: None,
triclinic_norm_square_p2: None,
triclinic_norm_square_p3: None,
monoclinic_norm_square_p1: None,
monoclinic_norm_square_p2: None,
monoclinic_norm_square_p3: None,
orthohombic_norm_square_p1: None,
orthohombic_norm_square_p2: None,
orthohombic_norm_square_p3: None,
tetragonal_norm_square_p1: None,
tetragonal_norm_square_p2: None,
tetragonal_norm_square_p3: None,
hexagonal_norm_square_p1: None,
hexagonal_norm_square_p2: None,
hexagonal_norm_square_p3: None,
isotropic_norm_square: None,
};
let particle_info_file = File::open(particle_info_file).unwrap();
let buf_reader = BufReader::new(particle_info_file);
let mut rdr = csv::ReaderBuilder::new()
.has_headers(true)
.delimiter(b' ')
.from_reader(buf_reader);
for result in rdr.deserialize() {
let record: ParticleRecord = result.unwrap();
if record.id == *particle_id {
particle_record = record;
}
}
println!("end retrieve antisotropy info");
println!("create lambert equal area gridpoint");
let lambert =
create_lambert_equal_area_gridpoint(sphere_points, "upper".to_string())
.unwrap();
println!("create sphere_point_grid");
let mut sphere_point_grid =
Array2::zeros((3, sphere_points * sphere_points));
for i in 0..sphere_points {
for j in 0..sphere_points {
sphere_point_grid[[0, i * sphere_points + j]] = lambert.x[[i, j]];
sphere_point_grid[[1, i * sphere_points + j]] = lambert.y[[i, j]];
sphere_point_grid[[2, i * sphere_points + j]] = lambert.z[[i, j]];
}
}
let n_grains = particle_olivine_a_axis_vectors.len();
let mut pole_figure_grid: Vec<Vec<PoleFigure>> = vec![
vec![
PoleFigure {
crystal_axis: CrystalAxes::AAxis,
mineral: Mineral::Olivine,
counts: Array2::zeros((n_grains, 3)),
max_count: 0.0,
};
pole_figure_configuration.minerals.len()
];
pole_figure_configuration.axes.len()
];
let mut figure_horizontal_axis = 0;
for axis in pole_figure_configuration.axes.clone() {
let mut figure_vertical_axis = 0;
for mineral in pole_figure_configuration.minerals.clone() {
let mut particle_arrays = Array2::zeros((n_grains, 3));
for i in 0..n_grains {
for j in 0..3 {
particle_arrays[[i, j]] = match axis {
CrystalAxes::AAxis => match mineral {
Mineral::Olivine => {
particle_olivine_a_axis_vectors[i][j]
}
Mineral::Enstatite => {
particle_enstatite_a_axis_vectors[i][j]
}
},
CrystalAxes::BAxis => match mineral {
Mineral::Olivine => {
particle_olivine_b_axis_vectors[i][j]
}
Mineral::Enstatite => {
particle_enstatite_b_axis_vectors[i][j]
}
},
CrystalAxes::CAxis => match mineral {
Mineral::Olivine => {
particle_olivine_c_axis_vectors[i][j]
}
Mineral::Enstatite => {
particle_enstatite_c_axis_vectors[i][j]
}
},
};
}
}
let counts = gaussian_orientation_counts(
&particle_arrays,
&sphere_point_grid,
sphere_points,
)
.unwrap();
let mut max_count_value = 0.0;
for i in 0..counts.shape()[0] - 1 {
for j in 0..counts.shape()[1] - 1 {
if counts[[i, j]] > max_count_value {
max_count_value = counts[[i, j]];
}
}
}
pole_figure_grid[figure_horizontal_axis][figure_vertical_axis] =
PoleFigure {
crystal_axis: axis.clone(),
mineral: mineral.clone(),
counts: counts,
max_count: max_count_value,
};
figure_vertical_axis += 1;
}
figure_horizontal_axis += 1;
}
let mut figure_vertical_axis = 0;
for _mineral in pole_figure_configuration.minerals.clone() {
let mut max_count_value = 0.0;
let mut figure_horizontal_axis = 0;
for _axis in pole_figure_configuration.axes.clone() {
if pole_figure_grid[figure_horizontal_axis][figure_vertical_axis]
.max_count
> max_count_value
{
max_count_value = pole_figure_grid[figure_horizontal_axis]
[figure_vertical_axis]
.max_count;
}
figure_horizontal_axis += 1;
}
figure_horizontal_axis = 0;
for _axis in pole_figure_configuration.axes.clone() {
pole_figure_grid[figure_horizontal_axis][figure_vertical_axis]
.max_count = max_count_value;
figure_horizontal_axis += 1;
}
figure_vertical_axis += 1;
}
make_pole_figures(
pole_figure_configuration.small_figure,
pole_figure_configuration.no_description_text,
elastisity_header,
n_grains,
0,
&pole_figure_grid,
&lambert,
output_file,
&particle_record,
time,
gam,
&pole_figure_configuration.color_scale,
)
.unwrap();
println!(
" After make_polefigures: Elapsed time: {:.2?}",
before.elapsed()
);
}
println!("go to next id");
}
}
}
});
Ok(())
}
fn euler_angles_to_rotation_matrix(
euler_angles: Array<f64, ndarray::Dim<[usize; 1]>>,
) -> Result<Array2<f64>, Box<dyn std::error::Error>> {
let mut rotation_matrix: Array2<f64> = Array::zeros((3, 3));
rotation_matrix[[0, 0]] = euler_angles[2].cos() * euler_angles[0].cos()
- euler_angles[1].cos() * euler_angles[0].sin() * euler_angles[2].sin();
rotation_matrix[[0, 1]] = -euler_angles[2].cos() * euler_angles[0].sin()
- euler_angles[1].cos() * euler_angles[0].cos() * euler_angles[2].sin();
rotation_matrix[[0, 2]] = -euler_angles[2].sin() * euler_angles[1].sin();
rotation_matrix[[1, 0]] = euler_angles[2].sin() * euler_angles[0].cos()
+ euler_angles[1].cos() * euler_angles[0].sin() * euler_angles[2].cos();
rotation_matrix[[1, 1]] = -euler_angles[2].sin() * euler_angles[0].sin()
+ euler_angles[1].cos() * euler_angles[0].cos() * euler_angles[2].cos();
rotation_matrix[[1, 2]] = euler_angles[2].cos() * euler_angles[1].sin();
rotation_matrix[[2, 0]] = -euler_angles[1].sin() * euler_angles[0].sin();
rotation_matrix[[2, 1]] = -euler_angles[1].sin() * euler_angles[0].cos();
rotation_matrix[[2, 2]] = euler_angles[1].cos();
Ok(rotation_matrix)
}
fn gaussian_orientation_counts(
particles: &Array2<f64>,
sphere_point_grid: &Array2<f64>,
sphere_points: usize,
) -> Result<Array2<f64>, Box<dyn std::error::Error>> {
let npts = particles.shape()[0];
let k = (2. * (1. + npts as f64 / 9.)).min(100.);
let std_dev = ((npts as f64 * (k as f64 / 2. - 1.) / (k * k)) as f64).sqrt();
let mut cosalpha = particles.dot(sphere_point_grid);
cosalpha.par_mapv_inplace(f64::abs);
cosalpha = (k as f64) * (cosalpha - 1.);
cosalpha.par_mapv_inplace(f64::exp);
let counts = cosalpha.sum_axis(Axis(0));
let counts = counts.into_shape((sphere_points, sphere_points))?;
let counts = counts / (3. * std_dev);
Ok(counts)
}
#[cfg(test)]
mod tests {
use super::*;
fn euler_angles_from_rotation_matrix(
rotation_matrix: Array2<f64>,
) -> Array<f64, ndarray::Dim<[usize; 1]>> {
let mut euler_angles = Array::zeros(3);
let theta = rotation_matrix[[2, 2]].acos();
let mut phi1 = 0.0;
let phi2;
if theta != 0.0 && theta != std::f64::consts::PI {
phi1 = (rotation_matrix[[2, 0]] / -theta.sin())
.atan2(rotation_matrix[[2, 1]] / -theta.sin());
phi2 = (rotation_matrix[[0, 2]] / -theta.sin())
.atan2(rotation_matrix[[1, 2]] / theta.sin());
} else {
if theta == 0.0 {
phi2 = -phi1 - rotation_matrix[[0, 1]].atan2(rotation_matrix[[0, 0]]);
} else {
phi2 = phi1 + rotation_matrix[[0, 1]].atan2(rotation_matrix[[0, 0]]);
}
}
euler_angles[0] = phi1;
euler_angles[1] = theta;
euler_angles[2] = phi2;
euler_angles[0] = euler_angles[0]
- 2.0 * std::f64::consts::PI * (euler_angles[0] / (2.0 * std::f64::consts::PI)).floor();
euler_angles[1] = euler_angles[1]
- 2.0 * std::f64::consts::PI * (euler_angles[1] / (2.0 * std::f64::consts::PI)).floor();
euler_angles[2] = euler_angles[2]
- 2.0 * std::f64::consts::PI * (euler_angles[2] / (2.0 * std::f64::consts::PI)).floor();
euler_angles
}
#[test]
fn test_euler_angle_and_rotation_matrix_functions_part_1() {
let mut rot1 = Array2::zeros((3, 3));
rot1[[0, 0]] = 0.36;
rot1[[0, 1]] = 0.48;
rot1[[0, 2]] = -0.8;
rot1[[1, 0]] = -0.8;
rot1[[1, 1]] = 0.6;
rot1[[1, 2]] = 0.;
rot1[[2, 0]] = 0.48;
rot1[[2, 1]] = 0.64;
rot1[[2, 2]] = 0.6;
let ea1 = euler_angles_from_rotation_matrix(rot1.clone());
let rot2 = euler_angles_to_rotation_matrix(ea1).unwrap();
assert!(
(rot1[[0, 0]] - rot2[[0, 0]]).abs() < 1e-8,
"Difference between rot1[0,0] ({}) and rot2[0,0] ({}) is too large: {}.",
rot1[[0, 0]],
rot2[[0, 0]],
(rot1[[0, 0]] - rot2[[0, 0]]).abs()
);
assert!(
(rot1[[0, 1]] - rot2[[0, 1]]).abs() < 1e-8,
"Difference between rot1[0,1] ({}) and rot2[0,1] ({}) is too large: {}.",
rot1[[0, 1]],
rot2[[0, 1]],
(rot1[[0, 1]] - rot2[[0, 1]]).abs()
);
assert!(
(rot1[[0, 2]] - rot2[[0, 2]]).abs() < 1e-8,
"Difference between rot1[0,2] ({}) and rot2[0,2] ({}) is too large: {}.",
rot1[[0, 2]],
rot2[[0, 2]],
(rot1[[0, 2]] - rot2[[0, 2]]).abs()
);
assert!(
(rot1[[1, 0]] - rot2[[1, 0]]).abs() < 1e-8,
"Difference between rot1[1,0] ({}) and rot2[1,0] ({}) is too large: {}.",
rot1[[1, 0]],
rot2[[1, 0]],
(rot1[[1, 0]] - rot2[[1, 0]]).abs()
);
assert!(
(rot1[[1, 1]] - rot2[[1, 1]]).abs() < 1e-8,
"Difference between rot1[1,1] ({}) and rot2[1,1] ({}) is too large: {}.",
rot1[[1, 1]],
rot2[[1, 1]],
(rot1[[1, 1]] - rot2[[1, 1]]).abs()
);
assert!(
(rot1[[1, 2]] - rot2[[1, 2]]).abs() < 1e-8,
"Difference between rot1[1,2] ({}) and rot2[1,2] ({}) is too large: {}.",
rot1[[1, 2]],
rot2[[1, 2]],
(rot1[[1, 2]] - rot2[[1, 2]]).abs()
);
assert!(
(rot1[[2, 0]] - rot2[[2, 0]]).abs() < 1e-8,
"Difference between rot1[2,0] ({}) and rot2[2,0] ({}) is too large: {}.",
rot1[[2, 0]],
rot2[[2, 0]],
(rot1[[2, 0]] - rot2[[2, 0]]).abs()
);
assert!(
(rot1[[2, 1]] - rot2[[2, 1]]).abs() < 1e-8,
"Difference between rot1[2,1] ({}) and rot2[2,1] ({}) is too large: {}.",
rot1[[2, 1]],
rot2[[2, 1]],
(rot1[[2, 1]] - rot2[[2, 1]]).abs()
);
assert!(
(rot1[[2, 2]] - rot2[[2, 2]]).abs() < 1e-8,
"Difference between rot1[2,2] ({}) and rot2[2,2] ({}) is too large: {}.",
rot1[[2, 2]],
rot2[[2, 2]],
(rot1[[2, 2]] - rot2[[2, 2]]).abs()
);
let ea2 = euler_angles_from_rotation_matrix(rot2.clone());
let rot3 = euler_angles_to_rotation_matrix(ea2).unwrap();
assert!(
(rot1[[0, 0]] - rot3[[0, 0]]).abs() < 1e-8,
"Difference between rot1[0,0] ({}) and rot3[0,0] ({}) is too large: {}.",
rot1[[0, 0]],
rot3[[0, 0]],
(rot1[[0, 0]] - rot3[[0, 0]]).abs()
);
assert!(
(rot1[[0, 1]] - rot3[[0, 1]]).abs() < 1e-8,
"Difference between rot1[0,1] ({}) and rot3[0,1] ({}) is too large: {}.",
rot1[[0, 1]],
rot3[[0, 1]],
(rot1[[0, 1]] - rot3[[0, 1]]).abs()
);
assert!(
(rot1[[0, 2]] - rot3[[0, 2]]).abs() < 1e-8,
"Difference between rot1[0,2] ({}) and rot3[0,2] ({}) is too large: {}.",
rot1[[0, 2]],
rot3[[0, 2]],
(rot1[[0, 2]] - rot3[[0, 2]]).abs()
);
assert!(
(rot1[[1, 0]] - rot3[[1, 0]]).abs() < 1e-8,
"Difference between rot1[1,0] ({}) and rot3[1,0] ({}) is too large: {}.",
rot1[[1, 0]],
rot3[[1, 0]],
(rot1[[1, 0]] - rot3[[1, 0]]).abs()
);
assert!(
(rot1[[1, 1]] - rot3[[1, 1]]).abs() < 1e-8,
"Difference between rot1[1,1] ({}) and rot3[1,1] ({}) is too large: {}.",
rot1[[1, 1]],
rot3[[1, 1]],
(rot1[[1, 1]] - rot3[[1, 1]]).abs()
);
assert!(
(rot1[[1, 2]] - rot3[[1, 2]]).abs() < 1e-8,
"Difference between rot1[1,2] ({}) and rot3[1,2] ({}) is too large: {}.",
rot1[[1, 2]],
rot3[[1, 2]],
(rot1[[1, 2]] - rot3[[1, 2]]).abs()
);
assert!(
(rot1[[2, 0]] - rot3[[2, 0]]).abs() < 1e-8,
"Difference between rot1[2,0] ({}) and rot3[2,0] ({}) is too large: {}.",
rot1[[2, 0]],
rot3[[2, 0]],
(rot1[[2, 0]] - rot3[[2, 0]]).abs()
);
assert!(
(rot1[[2, 1]] - rot3[[2, 1]]).abs() < 1e-8,
"Difference between rot1[2,1] ({}) and rot3[2,1] ({}) is too large: {}.",
rot1[[2, 1]],
rot3[[2, 1]],
(rot1[[2, 1]] - rot3[[2, 1]]).abs()
);
assert!(
(rot1[[2, 2]] - rot3[[2, 2]]).abs() < 1e-8,
"Difference between rot1[2,2] ({}) and rot3[2,2] ({}) is too large: {}.",
rot1[[2, 2]],
rot3[[2, 2]],
(rot1[[2, 2]] - rot3[[2, 2]]).abs()
);
}
#[test]
fn test_euler_angle_and_rotation_matrix_functions_part_2() {
let mut rot1 = Array2::zeros((3, 3));
rot1[[0, 0]] = 0.36;
rot1[[0, 1]] = 0.48;
rot1[[0, 2]] = -0.8;
rot1[[1, 0]] = -0.8;
rot1[[1, 1]] = 0.6;
rot1[[1, 2]] = 0.;
rot1[[2, 0]] = 0.0;
rot1[[2, 1]] = 0.0;
rot1[[2, 2]] = 0.0;
let mut rot2_expected = Array2::zeros((3, 3));
rot2_expected[[0, 0]] = 0.;
rot2_expected[[0, 1]] = 0.;
rot2_expected[[0, 2]] = -1.;
rot2_expected[[1, 0]] = -1.;
rot2_expected[[1, 1]] = 0.0;
rot2_expected[[1, 2]] = 0.0;
rot2_expected[[2, 0]] = 0.0;
rot2_expected[[2, 1]] = 1.0;
rot2_expected[[2, 2]] = 0.0;
let ea1 = euler_angles_from_rotation_matrix(rot1.clone());
let rot2 = euler_angles_to_rotation_matrix(ea1).unwrap();
assert!(
(rot2_expected[[0, 0]] - rot2[[0, 0]]).abs() < 1e-8,
"Difference between rot2_expected[0,0] ({}) and rot2[0,0] ({}) is too large: {}.",
rot2_expected[[0, 0]],
rot2[[0, 0]],
(rot2_expected[[0, 0]] - rot2[[0, 0]]).abs()
);
assert!(
(rot2_expected[[0, 1]] - rot2[[0, 1]]).abs() < 1e-8,
"Difference between rot2_expected[0,1] ({}) and rot2[0,1] ({}) is too large: {}.",
rot2_expected[[0, 1]],
rot2[[0, 1]],
(rot2_expected[[0, 1]] - rot2[[0, 1]]).abs()
);
assert!(
(rot2_expected[[0, 2]] - rot2[[0, 2]]).abs() < 1e-8,
"Difference between rot2_expected[0,2] ({}) and rot2[0,2] ({}) is too large: {}.",
rot2_expected[[0, 2]],
rot2[[0, 2]],
(rot2_expected[[0, 2]] - rot2[[0, 2]]).abs()
);
assert!(
(rot2_expected[[1, 0]] - rot2[[1, 0]]).abs() < 1e-8,
"Difference between rot2_expected[1,0] ({}) and rot2[1,0] ({}) is too large: {}.",
rot2_expected[[1, 0]],
rot2[[1, 0]],
(rot2_expected[[1, 0]] - rot2[[1, 0]]).abs()
);
assert!(
(rot2_expected[[1, 1]] - rot2[[1, 1]]).abs() < 1e-8,
"Difference between rot2_expected[1,1] ({}) and rot2[1,1] ({}) is too large: {}.",
rot2_expected[[1, 1]],
rot2[[1, 1]],
(rot2_expected[[1, 1]] - rot2[[1, 1]]).abs()
);
assert!(
(rot2_expected[[1, 2]] - rot2[[1, 2]]).abs() < 1e-8,
"Difference between rot2_expected[1,2] ({}) and rot2[1,2] ({}) is too large: {}.",
rot2_expected[[1, 2]],
rot2[[1, 2]],
(rot2_expected[[1, 2]] - rot2[[1, 2]]).abs()
);
assert!(
(rot2_expected[[2, 0]] - rot2[[2, 0]]).abs() < 1e-8,
"Difference between rot2_expected[2,0] ({}) and rot2[2,0] ({}) is too large: {}.",
rot2_expected[[2, 0]],
rot2[[2, 0]],
(rot2_expected[[2, 0]] - rot2[[2, 0]]).abs()
);
assert!(
(rot2_expected[[2, 1]] - rot2[[2, 1]]).abs() < 1e-8,
"Difference between rot2_expected[2,1] ({}) and rot2[2,1] ({}) is too large: {}.",
rot2_expected[[2, 1]],
rot2[[2, 1]],
(rot2_expected[[2, 1]] - rot2[[2, 1]]).abs()
);
assert!(
(rot2_expected[[2, 2]] - rot2[[2, 2]]).abs() < 1e-8,
"Difference between rot2_expected[2,2] ({}) and rot2[2,2] ({}) is too large: {}.",
rot2_expected[[2, 2]],
rot2[[2, 2]],
(rot2_expected[[2, 2]] - rot2[[2, 2]]).abs()
);
let ea2 = euler_angles_from_rotation_matrix(rot2.clone());
let rot3 = euler_angles_to_rotation_matrix(ea2).unwrap();
assert!(
(rot2_expected[[0, 0]] - rot3[[0, 0]]).abs() < 1e-8,
"Difference between rot2_expected[0,0] ({}) and rot3[0,0] ({}) is too large: {}.",
rot2_expected[[0, 0]],
rot3[[0, 0]],
(rot2_expected[[0, 0]] - rot3[[0, 0]]).abs()
);
assert!(
(rot2_expected[[0, 1]] - rot3[[0, 1]]).abs() < 1e-8,
"Difference between rot2_expected[0,1] ({}) and rot3[0,1] ({}) is too large: {}.",
rot2_expected[[0, 1]],
rot3[[0, 1]],
(rot2_expected[[0, 1]] - rot3[[0, 1]]).abs()
);
assert!(
(rot2_expected[[0, 2]] - rot3[[0, 2]]).abs() < 1e-8,
"Difference between rot2_expected[0,2] ({}) and rot3[0,2] ({}) is too large: {}.",
rot2_expected[[0, 2]],
rot3[[0, 2]],
(rot2_expected[[0, 2]] - rot3[[0, 2]]).abs()
);
assert!(
(rot2_expected[[1, 0]] - rot3[[1, 0]]).abs() < 1e-8,
"Difference between rot2_expected[1,0] ({}) and rot3[1,0] ({}) is too large: {}.",
rot2_expected[[1, 0]],
rot3[[1, 0]],
(rot2_expected[[1, 0]] - rot3[[1, 0]]).abs()
);
assert!(
(rot2_expected[[1, 1]] - rot3[[1, 1]]).abs() < 1e-8,
"Difference between rot2_expected[1,1] ({}) and rot3[1,1] ({}) is too large: {}.",
rot2_expected[[1, 1]],
rot3[[1, 1]],
(rot2_expected[[1, 1]] - rot3[[1, 1]]).abs()
);
assert!(
(rot2_expected[[1, 2]] - rot3[[1, 2]]).abs() < 1e-8,
"Difference between rot2_expected[1,2] ({}) and rot3[1,2] ({}) is too large: {}.",
rot2_expected[[1, 2]],
rot3[[1, 2]],
(rot2_expected[[1, 2]] - rot3[[1, 2]]).abs()
);
assert!(
(rot2_expected[[2, 0]] - rot3[[2, 0]]).abs() < 1e-8,
"Difference between rot2_expected[2,0] ({}) and rot3[2,0] ({}) is too large: {}.",
rot2_expected[[2, 0]],
rot3[[2, 0]],
(rot2_expected[[2, 0]] - rot3[[2, 0]]).abs()
);
assert!(
(rot2_expected[[2, 1]] - rot3[[2, 1]]).abs() < 1e-8,
"Difference between rot2_expected[2,1] ({}) and rot3[2,1] ({}) is too large: {}.",
rot2_expected[[2, 1]],
rot3[[2, 1]],
(rot2_expected[[2, 1]] - rot3[[2, 1]]).abs()
);
assert!(
(rot2_expected[[2, 2]] - rot3[[2, 2]]).abs() < 1e-8,
"Difference between rot2_expected[2,2] ({}) and rot3[2,2] ({}) is too large: {}.",
rot2_expected[[2, 2]],
rot3[[2, 2]],
(rot2_expected[[2, 2]] - rot3[[2, 2]]).abs()
);
}
}