qsym2/auxiliary/
geometry.rs

1//! Geometrical objects and manipulations.
2
3use std::collections::HashSet;
4use std::fmt;
5
6use approx;
7use derive_builder::Builder;
8use fraction;
9use itertools::{self, Itertools};
10use nalgebra::{ClosedMulAssign, Matrix3, Point3, Rotation3, Scalar, UnitVector3, Vector3};
11use num_traits::{One, ToPrimitive};
12use serde::{Deserialize, Serialize};
13
14use crate::auxiliary::atom::Atom;
15use crate::auxiliary::misc::HashableFloat;
16
17type F32 = fraction::GenericFraction<u32>;
18
19#[cfg(test)]
20#[path = "geometry_tests.rs"]
21mod geometry_tests;
22
23// ================
24// Enum definitions
25// ================
26
27/// Enumerated type to classify the type of improper rotation given an angle and axis.
28pub enum ImproperRotationKind {
29    /// The improper rotation is a rotation by the specified angle and axis followed by a
30    /// reflection in a mirror plane perpendicular to the axis.
31    MirrorPlane,
32
33    /// The improper rotation is a rotation by the specified angle and axis followed by an
34    /// inversion through the centre of inversion.
35    InversionCentre,
36}
37
38/// Mirror-plane improper rotation kind.
39pub const IMSIG: ImproperRotationKind = ImproperRotationKind::MirrorPlane;
40
41/// Inversion-centre improper rotation kind.
42pub const IMINV: ImproperRotationKind = ImproperRotationKind::InversionCentre;
43
44// =================
45// Utility functions
46// =================
47
48/// Returns the rotation angle adjusted to be in the interval $`(-\pi, +\pi]`$ and the number of
49/// $`2\pi`$-folds required to bring the original angle to that interval.
50///
51/// # Arguments
52///
53/// * `rot_ang` - A rotation angle.
54/// * `thresh` - A threshold for comparisons.
55///
56/// # Returns
57///
58/// The normalised rotation angle and the number of folds.
59#[must_use]
60pub fn normalise_rotation_angle(rot_ang: f64, thresh: f64) -> (f64, u32) {
61    let frac_1_2 = 1.0 / 2.0;
62    let fraction = rot_ang / (2.0 * std::f64::consts::PI);
63    if fraction > frac_1_2 + thresh {
64        let integer_part = fraction.trunc().to_u32().unwrap_or_else(|| {
65            panic!("Unable to convert the integer part of `{fraction}` to `u32`.")
66        });
67        let x = if fraction.fract() <= frac_1_2 + thresh {
68            integer_part
69        } else {
70            integer_part + 1
71        };
72        (rot_ang - 2.0 * std::f64::consts::PI * f64::from(x), x)
73    } else if fraction <= -frac_1_2 + thresh {
74        let integer_part = (-fraction).trunc().to_u32().unwrap_or_else(|| {
75            panic!("Unable to convert the integer part of `{fraction}` to `u32`.")
76        });
77        let x = if (-fraction).fract() < frac_1_2 - thresh {
78            integer_part
79        } else {
80            integer_part + 1
81        };
82        (rot_ang + 2.0 * std::f64::consts::PI * f64::from(x), x)
83    } else {
84        (rot_ang, 0)
85    }
86}
87
88/// Returns the rotation fraction adjusted to be in the interval $`(-1/2, +1/2]`$ and the number of
89/// $`1`$-folds required to bring the original fraction to that interval.
90///
91/// # Arguments
92///
93/// * `fraction` - A rotation fraction.
94///
95/// # Returns
96///
97/// The normalised rotation fraction and the number of folds.
98#[must_use]
99pub fn normalise_rotation_fraction(fraction: F32) -> (F32, u32) {
100    // Consider a fraction f.
101    //
102    // If f > 1/2, we seek a positive integer x such that
103    //  -1/2 < f - x <= 1/2.
104    // It turns out that x ∈ [f - 1/2, f + 1/2).
105    //
106    // If f <= -1/2, we seek a positive integer x such that
107    //  -1/2 < f + x <= 1/2.
108    // It turns out that x ∈ (-f - 1/2, -f + 1/2].
109    //
110    // If the proper rotation corresponding to f is reached from the identity
111    // via a continuous path in the parametric ball, x gives the number of times
112    // this path goes through a podal-antipodal jump, and thus whether x is even
113    // corresponds to whether this homotopy path is of class 0.
114    //
115    // See S.L. Altmann, Rotations, Quaternions, and Double Groups (Dover
116    // Publications, Inc., New York, 2005) for further information.
117    let frac_1_2 = F32::new(1u32, 2u32);
118    if fraction > frac_1_2 {
119        let integer_part = fraction.trunc();
120        let x = if fraction.fract() <= frac_1_2 {
121            integer_part
122        } else {
123            integer_part + F32::one()
124        };
125        (
126            fraction - x,
127            x.to_u32()
128                .expect("Unable to convert the 2π-turn number to `u32`."),
129        )
130    } else if fraction <= -frac_1_2 {
131        let integer_part = (-fraction).trunc();
132        let x = if (-fraction).fract() < frac_1_2 {
133            integer_part
134        } else {
135            integer_part + F32::one()
136        };
137        (
138            fraction + x,
139            x.to_u32()
140                .expect("Unable to convert the 2π-turn number to `u32`."),
141        )
142    } else {
143        (fraction, 0)
144    }
145}
146
147/// Determines the reduced fraction $`k/n`$ where $`k`$ and $`n`$ are both integers representing a
148/// proper rotation $`C_n^k`$ corresponding to a specified rotation angle.
149///
150/// # Arguments
151///
152/// * `angle` - An angle of rotation.
153/// * `thresh` - A threshold for checking if a floating point number is integral.
154/// * `max_trial_power` - Maximum power $`k`$ to try.
155///
156/// # Returns
157///
158/// An [`Option`] wrapping the required fraction.
159///
160/// # Panics
161///
162/// Panics if the deduced order $`n`$ is negative.
163#[must_use]
164pub fn get_proper_fraction(angle: f64, thresh: f64, max_trial_power: u32) -> Option<F32> {
165    let (normalised_angle, _) = normalise_rotation_angle(angle, thresh);
166    let rational_order = (2.0 * std::f64::consts::PI) / normalised_angle.abs();
167    let mut power: u32 = 1;
168    while approx::relative_ne!(
169        rational_order * (f64::from(power)),
170        (rational_order * (f64::from(power))).round(),
171        max_relative = thresh,
172        epsilon = thresh
173    ) && power < max_trial_power
174    {
175        power += 1;
176    }
177    if approx::relative_eq!(
178        rational_order * (f64::from(power)),
179        (rational_order * (f64::from(power))).round(),
180        max_relative = thresh,
181        epsilon = thresh
182    ) {
183        let orderf64 = (rational_order * (f64::from(power))).round();
184        assert!(orderf64.is_sign_positive());
185        assert!(orderf64 <= f64::from(u32::MAX));
186        #[allow(clippy::cast_sign_loss, clippy::cast_possible_truncation)]
187        let order = orderf64 as u32;
188        if normalised_angle > 0.0 {
189            Some(F32::new(power, order))
190        } else {
191            Some(F32::new_neg(power, order))
192        }
193    } else {
194        None
195    }
196}
197
198/// Computes the outer product between two three-dimensional vectors.
199///
200/// # Arguments
201///
202/// * `vec1` - The first vector, $`\mathbf{v}_1`$.
203/// * `vec2` - The second vector, $`\mathbf{v}_2`$.
204///
205/// # Returns
206///
207/// The outer product $`\mathbf{v}_1 \otimes \mathbf{v}_2`$.
208fn outer<T: Scalar + ClosedMulAssign + Copy>(vec1: &Vector3<T>, vec2: &Vector3<T>) -> Matrix3<T> {
209    let outer_product_iter: Vec<T> = vec2
210        .iter()
211        .flat_map(|&item_x| vec1.iter().map(move |&item_y| item_x * item_y))
212        .collect();
213    Matrix3::from_iterator(outer_product_iter)
214}
215
216/// Returns a $`3 \times 3`$ rotation matrix in $`\mathbb{R}^3`$ corresponding to a rotation
217/// through `angle` about `axis` raised to the power `power`.
218///
219/// # Arguments
220///
221/// * `angle` - The angle of rotation.
222/// * `axis` - The axis of rotation.
223/// * `power` - The power of rotation.
224///
225/// # Returns
226///
227/// The rotation matrix.
228#[must_use]
229pub fn proper_rotation_matrix(angle: f64, axis: &Vector3<f64>, power: i8) -> Matrix3<f64> {
230    let normalised_axis = UnitVector3::new_normalize(*axis);
231    Rotation3::from_axis_angle(&normalised_axis, (f64::from(power)) * angle).into_inner()
232}
233
234/// Returns a $`3 \times 3`$ transformation matrix in $`\mathbb{R}^3`$ corresponding to an improper
235/// rotation through `angle` about `axis` raised to the power `power`.
236///
237/// # Arguments
238///
239/// * `angle` - The angle of rotation.
240/// * `axis` - The axis of rotation.
241/// * `power` - The power of transformation.
242/// * `kind` - The convention in which the improper rotation is defined.
243///
244/// # Returns
245///
246/// The transformation matrix.
247#[must_use]
248pub fn improper_rotation_matrix(
249    angle: f64,
250    axis: &Vector3<f64>,
251    power: i8,
252    kind: &ImproperRotationKind,
253) -> Matrix3<f64> {
254    let rotmat = proper_rotation_matrix(angle, axis, power);
255    let normalised_axis = UnitVector3::new_normalize(*axis);
256    match kind {
257        ImproperRotationKind::MirrorPlane => {
258            let refmat = Matrix3::identity()
259                - 2.0 * (f64::from(power % 2)) * outer(&normalised_axis, &normalised_axis);
260            refmat * rotmat
261        }
262        ImproperRotationKind::InversionCentre => {
263            if power % 2 == 1 {
264                -rotmat
265            } else {
266                rotmat
267            }
268        }
269    }
270}
271
272/// Checks if a sequence of atoms are vertices of a regular polygon.
273///
274/// # Arguments
275///
276/// * `atoms` - A sequence of atoms to be tested.
277///
278/// # Returns
279///
280/// A flag indicating if the atoms form the vertices of a regular polygon.
281///
282/// # Panics
283///
284/// Panics if `atoms` contains fewer than three atoms.
285#[must_use]
286pub fn check_regular_polygon(atoms: &[&Atom]) -> bool {
287    assert!(
288        atoms.len() >= 3,
289        "Polygons can only be formed by three atoms or more."
290    );
291
292    let tot_m: f64 = atoms.iter().fold(0.0, |acc, atom| acc + atom.atomic_mass);
293    let com: Point3<f64> = atoms.iter().fold(Point3::origin(), |acc, atom| {
294        acc + (atom.coordinates * atom.atomic_mass - Point3::origin())
295    }) / tot_m;
296
297    let radial_dists: HashSet<(u64, i16, i8)> = atoms
298        .iter()
299        .map(|atom| {
300            (atom.coordinates - com)
301                .norm()
302                .round_factor(atom.threshold)
303                .integer_decode()
304        })
305        .collect();
306
307    // Check if all atoms are equidistant from the centre of mass
308    if radial_dists.len() == 1 {
309        let regular_angle = 2.0 * std::f64::consts::PI
310            / atoms
311                .len()
312                .to_f64()
313                .unwrap_or_else(|| panic!("Unable to convert `{}` to `f64`.", atoms.len()));
314        let thresh = atoms
315            .iter()
316            .fold(0.0_f64, |acc, atom| acc.max(atom.threshold));
317        let mut rad_vectors: Vec<Vector3<f64>> =
318            atoms.iter().map(|atom| atom.coordinates - com).collect();
319        let (vec_i, vec_j) = itertools::iproduct!(rad_vectors.iter(), rad_vectors.iter())
320            .max_by(|&(v_i1, v_j1), &(v_i2, v_j2)| {
321                v_i1.cross(v_j1)
322                    .norm()
323                    .partial_cmp(&v_i2.cross(v_j2).norm())
324                    .expect("Unable to compare the cross products of two vector pairs.")
325            })
326            .expect("Unable to find the vector pair with the largest norm cross product.");
327        let normal = UnitVector3::new_normalize(vec_i.cross(vec_j));
328        if normal.norm() < thresh {
329            return false;
330        }
331
332        let vec0 = atoms[0].coordinates - com;
333        rad_vectors.sort_by(|a, b| {
334            get_anticlockwise_angle(&vec0, a, &normal, thresh)
335                .partial_cmp(&get_anticlockwise_angle(&vec0, b, &normal, thresh))
336                .unwrap_or_else(|| {
337                    panic!(
338                        "Unable to compare anticlockwise angles of {a} and {b} relative to {vec0}."
339                    )
340                })
341        });
342        let vector_pairs: Vec<(&Vector3<f64>, &Vector3<f64>)> =
343            rad_vectors.iter().circular_tuple_windows().collect();
344        let mut angles: HashSet<(u64, i16, i8)> = vector_pairs
345            .iter()
346            .map(|(v1, v2)| {
347                get_anticlockwise_angle(v1, v2, &normal, thresh)
348                    .round_factor(thresh)
349                    .integer_decode()
350            })
351            .collect();
352        angles.insert(regular_angle.round_factor(thresh).integer_decode());
353
354        angles.len() == 1
355    } else {
356        false
357    }
358}
359
360/// Returns the anticlockwise angle $\phi$ from `vec1` to `vec2` when viewed down
361/// the `normal` vector.
362///
363/// This is only well-defined in $\mathbb{R}^3$. The range of the anticlockwise
364/// angle is $[0, 2\pi]$.
365///
366/// # Arguments
367///
368/// * `vec1` - The first vector.
369/// * `vec2` - The second vector.
370/// * `normal` - A normal unit vector defining the view.
371/// * `thresh` - Threshold for checking if either `vec1` or `vec2` is a null vector.
372///
373/// # Returns
374///
375/// The anticlockwise angle $\phi$.
376fn get_anticlockwise_angle(
377    vec1: &Vector3<f64>,
378    vec2: &Vector3<f64>,
379    normal: &UnitVector3<f64>,
380    thresh: f64,
381) -> f64 {
382    assert!(thresh >= std::f64::EPSILON);
383    assert!(vec1.norm() >= thresh);
384    assert!(vec2.norm() >= thresh);
385    let dot = vec1.dot(vec2);
386    let det = normal.into_inner().dot(&vec1.cross(vec2));
387    let mut angle = det.atan2(dot);
388    while angle < -thresh {
389        angle += 2.0 * std::f64::consts::PI;
390    }
391    angle
392}
393
394/// Geometrical transformability in three dimensions.
395pub trait Transform {
396    /// Transforms in-place the coordinates about the origin by a given
397    /// transformation.
398    ///
399    /// # Arguments
400    ///
401    /// * mat - A three-dimensional transformation matrix.
402    fn transform_mut(&mut self, mat: &Matrix3<f64>);
403
404    /// Rotates in-place the coordinates through `angle` about `axis`.
405    ///
406    /// # Arguments
407    ///
408    /// * angle - The angle of rotation.
409    /// * axis - The axis of rotation.
410    fn rotate_mut(&mut self, angle: f64, axis: &Vector3<f64>);
411
412    /// Improper-rotates in-place the coordinates through `angle` about `axis`.
413    ///
414    /// # Arguments
415    ///
416    /// * `angle` - The angle of rotation.
417    /// * `axis` - The axis of rotation.
418    /// * `kind` - The convention in which the improper rotation is defined.
419    fn improper_rotate_mut(&mut self, angle: f64, axis: &Vector3<f64>, kind: &ImproperRotationKind);
420
421    /// Translates in-place the coordinates by a specified translation vector in
422    /// three dimensions.
423    ///
424    /// # Arguments
425    ///
426    /// * `tvec` - The translation vector.
427    fn translate_mut(&mut self, tvec: &Vector3<f64>);
428
429    /// Recentres in-place to put the centre of mass at the origin.
430    fn recentre_mut(&mut self);
431
432    /// Reverses time by reversing in-place the polarity of any magnetic special atoms.
433    fn reverse_time_mut(&mut self);
434
435    /// Clones and transforms the coordinates about the origin by a given
436    /// transformation.
437    ///
438    /// # Arguments
439    ///
440    /// * `mat` - A three-dimensional transformation matrix.
441    ///
442    /// # Returns
443    ///
444    /// A transformed copy.
445    #[must_use]
446    fn transform(&self, mat: &Matrix3<f64>) -> Self;
447
448    /// Clones and rotates the coordinates through `angle` about `axis`.
449    ///
450    /// # Arguments
451    ///
452    /// * `angle` - The angle of rotation.
453    /// * `axis` - The axis of rotation.
454    ///
455    /// # Returns
456    ///
457    /// A rotated copy.
458    #[must_use]
459    fn rotate(&self, angle: f64, axis: &Vector3<f64>) -> Self;
460
461    /// Clones and improper-rotates the coordinates through `angle` about `axis`.
462    ///
463    /// # Arguments
464    ///
465    /// * `angle` - The angle of rotation.
466    /// * `axis` - The axis of rotation.
467    /// * `kind` - The convention in which the improper rotation is defined.
468    ///
469    /// # Returns
470    ///
471    /// An improper-rotated copy.
472    #[must_use]
473    fn improper_rotate(&self, angle: f64, axis: &Vector3<f64>, kind: &ImproperRotationKind)
474        -> Self;
475
476    /// Clones and translates in-place the coordinates by a specified
477    /// translation in three dimensions.
478    ///
479    /// # Arguments
480    ///
481    /// * `tvec` - The translation vector.
482    ///
483    /// # Returns
484    ///
485    /// A translated copy.
486    #[must_use]
487    fn translate(&self, tvec: &Vector3<f64>) -> Self;
488
489    /// Clones and recentres to put the centre of mass at the origin.
490    ///
491    /// # Returns
492    ///
493    /// A recentred copy.
494    #[must_use]
495    fn recentre(&self) -> Self;
496
497    /// Clones the molecule and reverses time by reversing the polarity of any magnetic special
498    /// atoms.
499    ///
500    /// # Returns
501    ///
502    /// A time-reversed copy.
503    #[must_use]
504    fn reverse_time(&self) -> Self;
505}
506
507// ===================
508// Positive Hemisphere
509// ===================
510
511// ----------------
512// ImproperOrdering
513// ----------------
514
515/// Enumerated type to handle comparisons symbolically.
516#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
517enum ImproperOrdering {
518    Greater,
519    GreaterEqual,
520    Less,
521    LessEqual,
522    Equal,
523}
524
525impl fmt::Display for ImproperOrdering {
526    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
527        match self {
528            ImproperOrdering::Greater => write!(f, ">"),
529            ImproperOrdering::GreaterEqual => write!(f, "≥"),
530            ImproperOrdering::Less => write!(f, "<"),
531            ImproperOrdering::LessEqual => write!(f, "≤"),
532            ImproperOrdering::Equal => write!(f, "="),
533        }
534    }
535}
536
537// ---------
538// Cartesian
539// ---------
540
541/***
542Coordinates
543***/
544
545/// Enumerated type to handle Cartesian coordinates symbolically.
546#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
547enum CartesianCoordinate {
548    X,
549    Y,
550    Z,
551}
552
553impl CartesianCoordinate {
554    /// Converts a Cartesian coordinate to a numerical index.
555    fn to_index(&self) -> usize {
556        match self {
557            CartesianCoordinate::X => 0,
558            CartesianCoordinate::Y => 1,
559            CartesianCoordinate::Z => 2,
560        }
561    }
562}
563
564impl fmt::Display for CartesianCoordinate {
565    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
566        match self {
567            CartesianCoordinate::X => write!(f, "x"),
568            CartesianCoordinate::Y => write!(f, "y"),
569            CartesianCoordinate::Z => write!(f, "z"),
570        }
571    }
572}
573
574/***
575Conditions
576***/
577
578/// Structure to handle inequality conditions written in terms of Cartesian coordinates.
579#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
580pub struct CartesianConditions {
581    /// The Cartesian conditions. The condititions are satisfied if all of the tuples in any of the
582    /// inner vectors are satisfied.
583    conditions: Vec<Vec<(CartesianCoordinate, ImproperOrdering, f64)>>,
584}
585
586impl CartesianConditions {
587    /// Checks if a vector satisfies the current Cartesian conditions. The condititions are
588    /// satisfied if all of the tuples in any of the inner vectors are satisfied.
589    ///
590    /// # Arguments
591    ///
592    /// * `vec` - A vector to check.
593    /// * `thresh` - A threshold for numerical comparisons.
594    ///
595    /// # Returns
596    ///
597    /// A boolean indicating if `vec` satisfies the conditions.
598    fn check(&self, vec: &Vector3<f64>, thresh: f64) -> bool {
599        self.conditions.iter().any(|condition_set| {
600            condition_set.iter().all(|(i, order, target)| match order {
601                ImproperOrdering::Greater => vec[i.to_index()] > target + thresh,
602                ImproperOrdering::GreaterEqual => vec[i.to_index()] > target - thresh,
603                ImproperOrdering::Less => vec[i.to_index()] < target - thresh,
604                ImproperOrdering::LessEqual => vec[i.to_index()] < target + thresh,
605                ImproperOrdering::Equal => approx::relative_eq!(
606                    vec[i.to_index()],
607                    target,
608                    max_relative = thresh,
609                    epsilon = thresh
610                ),
611            })
612        })
613    }
614}
615
616impl Default for CartesianConditions {
617    fn default() -> Self {
618        Self {
619            conditions: vec![
620                vec![(CartesianCoordinate::Z, ImproperOrdering::Greater, 0.0)],
621                vec![
622                    (CartesianCoordinate::Z, ImproperOrdering::Equal, 0.0),
623                    (CartesianCoordinate::X, ImproperOrdering::Greater, 0.0),
624                ],
625                vec![
626                    (CartesianCoordinate::Z, ImproperOrdering::Equal, 0.0),
627                    (CartesianCoordinate::X, ImproperOrdering::Equal, 0.0),
628                    (CartesianCoordinate::Y, ImproperOrdering::Greater, 0.0),
629                ],
630            ],
631        }
632    }
633}
634
635impl fmt::Display for CartesianConditions {
636    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
637        writeln!(f, "Cartesian conditions:")?;
638        let conditions = self
639            .conditions
640            .iter()
641            .map(|condition_set| {
642                condition_set
643                    .iter()
644                    .map(|(i, order, target)| format!("{i} {order} {target}"))
645                    .join(", ")
646            })
647            .join("\n  or\n");
648        writeln!(f, "{conditions}")?;
649        Ok(())
650    }
651}
652
653// ---------
654// Spherical
655// ---------
656
657/***
658Coordinates
659***/
660
661/// Enumerated type to handle spherical angular coordinates.
662#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
663pub enum SphericalCoordinate {
664    Theta,
665    Phi,
666}
667
668impl fmt::Display for SphericalCoordinate {
669    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
670        match self {
671            SphericalCoordinate::Theta => write!(f, "θ"),
672            SphericalCoordinate::Phi => write!(f, "φ"),
673        }
674    }
675}
676
677/***
678Conditions
679***/
680
681/// Structure to handle inequality conditions written in terms of spherical angular coordinates.
682#[derive(Debug, Clone, Builder, PartialEq, Serialize, Deserialize)]
683pub struct SphericalConditions {
684    /// The polar axis relative to which the polar angle $`\theta`$ is defined.
685    #[builder(setter(custom))]
686    z_basis: Vector3<f64>,
687
688    /// The azimuthal axis relative to which the azimuthal angle $`\phi`$ is defined.
689    #[builder(setter(custom))]
690    x_basis: Vector3<f64>,
691
692    /// The spherical angular conditions. The condititions are satisfied if all of the tuples in
693    /// any of the inner vectors are satisfied.
694    #[builder(setter(custom))]
695    conditions: Vec<Vec<(SphericalCoordinate, ImproperOrdering, f64)>>,
696}
697
698impl SphericalConditionsBuilder {
699    fn z_basis(&mut self, z_bas: Vector3<f64>) -> &mut Self {
700        self.z_basis = Some(z_bas.normalize());
701        self
702    }
703
704    fn x_basis(&mut self, x_bas: Vector3<f64>) -> &mut Self {
705        self.x_basis = Some(x_bas.normalize());
706        self
707    }
708
709    fn conditions(
710        &mut self,
711        conds: &[Vec<(SphericalCoordinate, ImproperOrdering, f64)>],
712    ) -> &mut Self {
713        self.conditions = Some(conds.to_vec());
714        self
715    }
716}
717
718impl SphericalConditions {
719    /// Returns a builder to construct [`Self`].
720    fn builder() -> SphericalConditionsBuilder {
721        SphericalConditionsBuilder::default()
722    }
723
724    /// Returns a required angular component of a vector given the current set of spherical
725    /// conditions that define the polar and azimuthal axes.
726    ///
727    /// # Arguments
728    ///
729    /// * `vec` - A vector whose components are to be retrieved.
730    /// * `coord` - A spherical angular coordinate.
731    /// * `thresh` - A threshold for checking if cosines of angles are equal to $`\pm 1`$.
732    ///
733    /// # Returns
734    ///
735    /// The required component.
736    fn get_component(&self, vec: &Vector3<f64>, coord: &SphericalCoordinate, thresh: f64) -> f64 {
737        match coord {
738            SphericalCoordinate::Theta => {
739                let cos_theta = vec.dot(&self.z_basis) / (vec.norm() * self.z_basis.norm());
740                if approx::relative_eq!(cos_theta, 1.0, epsilon = thresh, max_relative = thresh) {
741                    1.0f64.acos()
742                } else if approx::relative_eq!(
743                    cos_theta,
744                    -1.0,
745                    epsilon = thresh,
746                    max_relative = thresh
747                ) {
748                    (-1.0f64).acos()
749                } else {
750                    cos_theta.acos()
751                }
752            }
753            SphericalCoordinate::Phi => {
754                let y_vector = self.z_basis.cross(&self.x_basis).normalize();
755                let xy_vec = vec - vec.dot(&self.z_basis) / self.z_basis.norm() * self.z_basis;
756                let sgn_y = xy_vec.dot(&y_vector).signum();
757                let cos_phi = xy_vec.dot(&self.x_basis) / (xy_vec.norm() * self.x_basis.norm());
758                if approx::relative_eq!(cos_phi, 1.0, epsilon = thresh, max_relative = thresh) {
759                    sgn_y * 1.0f64.acos()
760                } else if approx::relative_eq!(
761                    cos_phi,
762                    -1.0,
763                    epsilon = thresh,
764                    max_relative = thresh
765                ) {
766                    sgn_y * (-1.0f64).acos()
767                } else {
768                    sgn_y * cos_phi.acos()
769                }
770            }
771        }
772    }
773
774    /// Checks if a vector satisfies the current spherical angular conditions. The condititions are
775    /// satisfied if all of the tuples in any of the inner vectors are satisfied.
776    ///
777    /// # Arguments
778    ///
779    /// * `vec` - A vector to check.
780    /// * `thresh` - A threshold for numerical comparisons.
781    ///
782    /// # Returns
783    ///
784    /// A boolean indicating if `vec` satisfies the conditions.
785    fn check(&self, vec: &Vector3<f64>, thresh: f64) -> bool {
786        self.conditions.iter().any(|condition_set| {
787            condition_set.iter().all(|(i, order, target)| {
788                let component = self.get_component(vec, i, thresh);
789                match order {
790                    ImproperOrdering::Greater => component > target + thresh,
791                    ImproperOrdering::GreaterEqual => component > target - thresh,
792                    ImproperOrdering::Less => component < target - thresh,
793                    ImproperOrdering::LessEqual => component < target + thresh,
794                    ImproperOrdering::Equal => approx::relative_eq!(
795                        component,
796                        target,
797                        max_relative = thresh,
798                        epsilon = thresh
799                    ),
800                }
801            })
802        })
803    }
804
805    /// Constructs a positive hemisphere where the equator consists of an odd number of equal and
806    /// disjoint arcs.
807    ///
808    /// The centre of the first arc is always at $`\phi = 0`$. Each arc is open at the
809    /// smaller-$`\phi`$ end and closed at the larger-$`\phi`$ end. It can be shown (see below)
810    /// that, as `n` is odd, no arcs can cross between $`+\pi`$ and $`-\pi`$.
811    ///
812    /// For $`n`$ odd, the centres of the most-positive and most-negative arcs are given by
813    ///
814    /// ```math
815    ///     \pm \frac{2\pi}{n} \times \frac{n - 1}{2} = \pm \pi \times \frac{n - 1}{n}.
816    /// ```
817    ///
818    /// Each arc has width $`\pi / n`$, so the most positive or most negative arc
819    /// $`\phi`$-coordinate are
820    ///
821    /// ```math
822    ///     \pm \left( \pi \times \frac{n - 1}{n} + \frac{\pi}{2n} \right)
823    ///     = \pm \pi \frac{2n - 1}{2n},
824    /// ```
825    ///
826    /// thus showing clearly that the arcs never cross from $`+\pi`$ to $`-\pi`$ and *vice versa*.
827    ///
828    ///
829    /// # Arguments
830    ///
831    /// * `z_basis` - The polar axis.
832    /// * `x_basis` - The azimuthal axis.
833    /// * `n` - An odd number specifying the number of equal and disjoint arcs belonging to the
834    /// positive hemisphere on the equator.
835    ///
836    /// # Returns
837    ///
838    /// The required spherical angular conditions.
839    fn new_disjoint_equatorial_arcs(
840        z_basis: Vector3<f64>,
841        x_basis: Vector3<f64>,
842        n: usize,
843    ) -> Self {
844        assert!(n > 0 && n.rem_euclid(2) == 1);
845        let n_f64 = n
846            .to_f64()
847            .expect("Unable to convert the number of arcs to `f64`.");
848        let half_arc = std::f64::consts::PI / (2.0 * n_f64);
849        let sep = 2.0 * std::f64::consts::PI / n_f64;
850        let half_pi = 0.5 * std::f64::consts::PI;
851
852        let mut conditions = vec![vec![
853            (
854                SphericalCoordinate::Theta,
855                ImproperOrdering::GreaterEqual,
856                0.0,
857            ),
858            (SphericalCoordinate::Theta, ImproperOrdering::Less, half_pi),
859        ]];
860
861        let phi_conditions = (0..n)
862            .map(|i| {
863                let (centre, _) = normalise_rotation_angle(i.to_f64().unwrap() * sep, f64::EPSILON);
864                let min_exc = centre - half_arc;
865                let max_inc = centre + half_arc;
866                vec![
867                    (SphericalCoordinate::Theta, ImproperOrdering::Equal, half_pi),
868                    (SphericalCoordinate::Phi, ImproperOrdering::Greater, min_exc),
869                    (
870                        SphericalCoordinate::Phi,
871                        ImproperOrdering::LessEqual,
872                        max_inc,
873                    ),
874                ]
875            })
876            .collect_vec();
877
878        conditions.extend(phi_conditions.into_iter());
879        Self::builder()
880            .z_basis(z_basis)
881            .x_basis(x_basis)
882            .conditions(&conditions)
883            .build()
884            .expect("Unable to construct a set of spherical-coordinate conditions.")
885    }
886}
887
888impl Default for SphericalConditions {
889    fn default() -> Self {
890        let half_pi = 0.5 * std::f64::consts::PI;
891        let conditions = vec![
892            vec![
893                (
894                    SphericalCoordinate::Theta,
895                    ImproperOrdering::GreaterEqual,
896                    0.0,
897                ),
898                (SphericalCoordinate::Theta, ImproperOrdering::Less, half_pi),
899            ],
900            vec![
901                (SphericalCoordinate::Theta, ImproperOrdering::Equal, half_pi),
902                (
903                    SphericalCoordinate::Phi,
904                    ImproperOrdering::Greater,
905                    -half_pi,
906                ),
907                (
908                    SphericalCoordinate::Phi,
909                    ImproperOrdering::LessEqual,
910                    half_pi,
911                ),
912            ],
913        ];
914        Self::builder()
915            .z_basis(Vector3::z())
916            .x_basis(Vector3::x())
917            .conditions(&conditions)
918            .build()
919            .expect("Unable to construct a set of spherical-coordinate conditions.")
920    }
921}
922
923impl fmt::Display for SphericalConditions {
924    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
925        writeln!(f, "Spherical conditions:")?;
926        writeln!(f, "Polar axis     (z-like): {:?}", self.z_basis)?;
927        writeln!(f, "Azimuthal axis (x-like): {:?}", self.x_basis)?;
928        let conditions = self
929            .conditions
930            .iter()
931            .map(|condition_set| {
932                condition_set
933                    .iter()
934                    .map(|(i, order, target)| format!("{i} {order} {target}"))
935                    .join(", ")
936            })
937            .join("\n  or\n");
938        writeln!(f, "{conditions}")?;
939        Ok(())
940    }
941}
942
943// ------------------
944// PositiveHemisphere
945// ------------------
946
947/// Enumerated type to handle positive hemispheres in Cartesian or spherical conditions.
948#[derive(Debug, Clone, PartialEq, Serialize, Deserialize)]
949pub enum PositiveHemisphere {
950    Cartesian(CartesianConditions),
951    Spherical(SphericalConditions),
952}
953
954impl PositiveHemisphere {
955    /// Constructs a new standard positive hemisphere in the Cartesian form.
956    pub fn new_standard_cartesian() -> Self {
957        Self::Cartesian(CartesianConditions::default())
958    }
959
960    /// Constructs a new standard positive hemisphere in the spherical form.
961    pub fn new_standard_spherical() -> Self {
962        Self::Spherical(SphericalConditions::default())
963    }
964
965    /// Constructs a new positive hemisphere in the spherical form with equal and disjoint arcs on
966    /// the equator.
967    ///
968    /// # Arguments
969    ///
970    /// `z_basis` - The polar axis.
971    /// `x_basis` - The azimuthal axis.
972    /// `n` - An odd number specifying the number of equal and disjoint arcs belonging to the
973    /// positive hemisphere on the equator.
974    ///
975    /// # Returns
976    ///
977    /// The required positive hemisphere.
978    pub fn new_spherical_disjoint_equatorial_arcs(
979        z_basis: Vector3<f64>,
980        x_basis: Vector3<f64>,
981        n: usize,
982    ) -> Self {
983        Self::Spherical(SphericalConditions::new_disjoint_equatorial_arcs(
984            z_basis, x_basis, n,
985        ))
986    }
987
988    /// Check if a rotation axis is in the current positive hemisphere.
989    ///
990    /// # Arguments
991    ///
992    /// * axis - An axis of rotation.
993    /// * thresh - Threshold for comparisons.
994    ///
995    /// # Returns
996    ///
997    /// Returns `true` if `axis` is in the positive hemisphere.
998    ///
999    /// # Panics
1000    ///
1001    /// Panics if the axis is a null vector.
1002    pub fn check_positive_pole(&self, axis: &Vector3<f64>, thresh: f64) -> bool {
1003        let normalised_axis = axis.normalize();
1004        match self {
1005            PositiveHemisphere::Cartesian(cart_conditions) => {
1006                cart_conditions.check(&normalised_axis, thresh)
1007            }
1008            PositiveHemisphere::Spherical(sph_conditions) => {
1009                sph_conditions.check(&normalised_axis, thresh)
1010            }
1011        }
1012    }
1013
1014    /// Returns the positive pole of a rotation axis with respect to the current positive
1015    /// hemisphere.
1016    ///
1017    /// # Arguments
1018    ///
1019    /// * axis - An axis of rotation.
1020    /// * thresh - Threshold for comparisons.
1021    ///
1022    /// # Returns
1023    ///
1024    /// The positive pole of `axis`.
1025    ///
1026    /// # Panics
1027    ///
1028    /// Panics if the resulting pole is a null vector.
1029    pub fn get_positive_pole(&self, axis: &Vector3<f64>, thresh: f64) -> Vector3<f64> {
1030        let normalised_axis = axis.normalize();
1031        if self.check_positive_pole(&normalised_axis, thresh) {
1032            normalised_axis
1033        } else {
1034            -normalised_axis
1035        }
1036    }
1037}
1038
1039impl Default for PositiveHemisphere {
1040    fn default() -> Self {
1041        Self::Cartesian(CartesianConditions::default())
1042    }
1043}
1044
1045impl fmt::Display for PositiveHemisphere {
1046    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
1047        match self {
1048            PositiveHemisphere::Cartesian(cart_conds) => write!(f, "{cart_conds}"),
1049            PositiveHemisphere::Spherical(sph_conds) => write!(f, "{sph_conds}"),
1050        }
1051    }
1052}
1053
1054/// Returns the standard positive pole of a rotation axis.
1055///
1056/// The definition of standard positive poles can be found in S.L. Altmann, Rotations,
1057/// Quaternions, and Double Groups (Dover Publications, Inc., New York, 2005) (Chapter 9).
1058///
1059/// # Arguments
1060///
1061/// * axis - An axis of rotation (proper or improper).
1062/// * thresh - Threshold for comparisons.
1063///
1064/// # Returns
1065///
1066/// The positive pole of `axis`.
1067///
1068/// # Panics
1069///
1070/// Panics if the resulting pole is a null vector.
1071#[must_use]
1072pub fn get_standard_positive_pole(axis: &Vector3<f64>, thresh: f64) -> Vector3<f64> {
1073    let poshem = PositiveHemisphere::new_standard_cartesian();
1074    poshem.get_positive_pole(axis, thresh)
1075}
1076
1077/// Check if a rotation axis is in the standard positive hemisphere.
1078///
1079/// The definition of the standard positive hemisphere can be found in S.L. Altmann, Rotations,
1080/// Quaternions, and Double Groups (Dover Publications, Inc., New York, 2005) (Chapter 9).
1081///
1082/// # Arguments
1083///
1084/// * axis - An axis of rotation.
1085/// * thresh - Threshold for comparisons.
1086///
1087/// # Returns
1088///
1089/// Returns `true` if `axis` is in the positive hemisphere.
1090///
1091/// # Panics
1092///
1093/// Panics if the axis is a null vector.
1094#[must_use]
1095pub fn check_standard_positive_pole(axis: &Vector3<f64>, thresh: f64) -> bool {
1096    let poshem = PositiveHemisphere::new_standard_cartesian();
1097    poshem.check_positive_pole(axis, thresh)
1098}