1use std::fmt;
4
5use approx;
6use log;
7use nalgebra as na;
8use serde::{Deserialize, Serialize};
9
10#[cfg(test)]
11#[path = "rotsym_tests.rs"]
12mod rotsym_tests;
13
14#[derive(Clone, Debug, Serialize, Deserialize)]
17pub enum RotationalSymmetry {
18 Spherical,
20 OblatePlanar,
23 OblateNonPlanar,
26 ProlateLinear,
28 ProlateNonLinear,
31 AsymmetricPlanar,
34 AsymmetricNonPlanar,
37}
38
39impl fmt::Display for RotationalSymmetry {
40 fn fmt(&self, f: &mut fmt::Formatter<'_>) -> fmt::Result {
41 match self {
42 RotationalSymmetry::Spherical => write!(f, "Spherical"),
43 RotationalSymmetry::OblatePlanar => write!(f, "Oblate, planar"),
44 RotationalSymmetry::OblateNonPlanar => write!(f, "Oblate, non-planar"),
45 RotationalSymmetry::ProlateLinear => write!(f, "Prolate, linear"),
46 RotationalSymmetry::ProlateNonLinear => write!(f, "Prolate, non-linear"),
47 RotationalSymmetry::AsymmetricPlanar => write!(f, "Asymmetric, planar"),
48 RotationalSymmetry::AsymmetricNonPlanar => write!(f, "Asymmetric, non-planar"),
49 }
50 }
51}
52
53#[must_use]
68pub fn calc_rotational_symmetry(
69 inertia_tensor: &na::Matrix3<f64>,
70 thresh: f64,
71) -> RotationalSymmetry {
72 let moi_mat = inertia_tensor.symmetric_eigenvalues();
73 let mut moi: Vec<&f64> = moi_mat.iter().collect();
74 moi.sort_by(|a, b| {
75 a.partial_cmp(b)
76 .unwrap_or_else(|| panic!("`{a}` and `{b}` cannot be compared."))
77 });
78 log::debug!("Moments of inertia:");
79 for component in &moi {
80 log::debug!(" {component:+.14}");
81 }
82 if approx::relative_eq!(*moi[0], *moi[1], epsilon = thresh, max_relative = thresh) {
83 if approx::relative_eq!(*moi[1], *moi[2], epsilon = thresh, max_relative = thresh) {
84 return RotationalSymmetry::Spherical;
85 }
86 if approx::relative_eq!(
87 *moi[2],
88 *moi[0] + *moi[1],
89 epsilon = thresh,
90 max_relative = thresh
91 ) {
92 return RotationalSymmetry::OblatePlanar;
93 }
94 return RotationalSymmetry::OblateNonPlanar;
95 }
96 if approx::relative_eq!(*moi[1], *moi[2], epsilon = thresh, max_relative = thresh) {
97 if moi[0].abs() < thresh {
98 return RotationalSymmetry::ProlateLinear;
99 }
100 return RotationalSymmetry::ProlateNonLinear;
101 }
102 if approx::relative_eq!(
103 *moi[2],
104 *moi[0] + *moi[1],
105 epsilon = thresh,
106 max_relative = thresh
107 ) {
108 return RotationalSymmetry::AsymmetricPlanar;
109 }
110 RotationalSymmetry::AsymmetricNonPlanar
111}