hoomd_vector/quaternion.rs
1// Copyright (c) 2024-2026 The Regents of the University of Michigan.
2// Part of hoomd-rs, released under the BSD 3-Clause License.
3
4//! Implement [`Quaternion`] and related types.
5use serde::{Deserialize, Serialize};
6use std::{
7 fmt,
8 ops::{Add, AddAssign, Div, DivAssign, Mul, MulAssign, Sub, SubAssign},
9};
10
11use approxim::approx_derive::RelativeEq;
12use rand::{
13 Rng, RngExt,
14 distr::{Distribution, StandardUniform},
15};
16use rand_distr::StandardNormal;
17
18use crate::{Cartesian, Cross, Error, InnerProduct, Rotate, Rotation, RotationMatrix, Unit};
19
20/// Extended complex number.
21///
22/// A quaternion has a real value and three complex values, represented by scalar and 3-vector
23/// respectively:
24/// ```math
25/// \mathbf{q} = (s, \vec{v})
26/// ```
27///
28/// Looking for the quaternion representation of 3D rotations? See [`Versor`].
29///
30/// ## Constructing quaternions
31///
32/// Create a quaternion with an array of coordinates (`[scalar, vector_0, vector_1, vector_2]`).
33/// ```
34/// use hoomd_vector::Quaternion;
35///
36/// let q = Quaternion::from([1.0, 2.0, 3.0, 4.0]);
37/// assert_eq!(q.scalar, 1.0);
38/// assert_eq!(q.vector, [2.0, 3.0, 4.0].into());
39/// ```
40///
41/// ## Quaternion properties
42///
43/// Compute a quaternion's norm:
44/// ```
45/// use hoomd_vector::Quaternion;
46///
47/// let q = Quaternion::from([3.0, 0.0, 4.0, 0.0]);
48/// let norm = q.norm();
49/// assert_eq!(norm, 5.0);
50/// ```
51///
52/// Form the conjugate:
53/// ```
54/// use hoomd_vector::Quaternion;
55///
56/// let q = Quaternion::from([1.0, 2.0, 3.0, 4.0]);
57/// let q_star = q.conjugate();
58/// assert_eq!(q_star, [1.0, -2.0, -3.0, -4.0].into());
59/// ```
60///
61/// ## Operating on quaternions
62///
63/// All operation examples use the following two quaternions:
64/// ```
65/// use hoomd_vector::Quaternion;
66///
67/// let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
68/// let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
69/// ```
70///
71/// Addition:
72///
73/// ```
74/// # use hoomd_vector::Quaternion;
75/// # let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
76/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
77/// let c = a + b;
78/// assert_eq!(c, [-1.0, 4.0, 10.0, -3.0].into());
79/// ```
80///
81/// ```
82/// # use hoomd_vector::Quaternion;
83/// # let mut a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
84/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
85/// a += b;
86/// assert_eq!(a, [-1.0, 4.0, 10.0, -3.0].into());
87/// ```
88///
89/// Subtraction:
90///
91/// ```
92/// # use hoomd_vector::Quaternion;
93/// # let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
94/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
95/// let c = a - b;
96/// assert_eq!(c, [3.0, -8.0, 2.0, -5.0].into());
97/// ```
98///
99/// ```
100/// # use hoomd_vector::Quaternion;
101/// # let mut a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
102/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
103/// a -= b;
104/// assert_eq!(a, [3.0, -8.0, 2.0, -5.0].into());
105/// ```
106///
107/// Multiplication by a scalar:
108///
109/// ```
110/// # use hoomd_vector::Quaternion;
111/// # let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
112/// let c = a * 2.0;
113/// assert_eq!(c, [2.0, -4.0, 12.0, -8.0].into());
114/// ```
115///
116/// ```
117/// # use hoomd_vector::Quaternion;
118/// # let mut a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
119/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
120/// a *= 2.0;
121/// assert_eq!(a, [2.0, -4.0, 12.0, -8.0].into());
122/// ```
123///
124/// Division by a scalar:
125///
126/// ```
127/// # use hoomd_vector::Quaternion;
128/// # let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
129/// let c = a / 2.0;
130/// assert_eq!(c, [0.5, -1.0, 3.0, -2.0].into());
131/// ```
132///
133/// ```
134/// # use hoomd_vector::Quaternion;
135/// # let mut a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
136/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
137/// a /= 2.0;
138/// assert_eq!(a, [0.5, -1.0, 3.0, -2.0].into());
139/// ```
140///
141/// Quaternion multiplication:
142///
143/// ```
144/// # use hoomd_vector::Quaternion;
145/// # let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
146/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
147/// let c = a * b;
148/// assert_eq!(c, [-10.0, 32.0, -30.0, -35.0].into());
149/// ```
150///
151/// ```
152/// # use hoomd_vector::Quaternion;
153/// # let mut a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
154/// # let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
155/// a *= b;
156/// assert_eq!(a, [-10.0, 32.0, -30.0, -35.0].into());
157/// ```
158#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
159pub struct Quaternion {
160 /// Scalar component
161 pub scalar: f64,
162
163 /// Vector component
164 pub vector: Cartesian<3>,
165}
166
167impl Quaternion {
168 /// Construct a pure quaternion: $` (0, \vec{v}) `$
169 ///
170 /// # Example
171 ///
172 /// ```
173 /// use hoomd_vector::{Cartesian, Quaternion};
174 ///
175 /// let q = Quaternion::pure(Cartesian::from([1.0, -2.0, 4.0]));
176 ///
177 /// assert_eq!(q.scalar, 0.0);
178 /// assert_eq!(q.vector, [1.0, -2.0, 4.0].into());
179 /// ```
180 #[inline]
181 pub fn pure(vector: Cartesian<3>) -> Self {
182 Self {
183 scalar: 0.0,
184 vector,
185 }
186 }
187
188 /// The norm of the quaternion, squared.
189 /// ```math
190 /// |\mathbf{q}|^2
191 /// ```
192 ///
193 /// # Example
194 /// ```
195 /// use hoomd_vector::Quaternion;
196 ///
197 /// let q = Quaternion::from([3.0, 0.0, 4.0, 0.0]);
198 /// let norm_squared = q.norm_squared();
199 /// assert_eq!(norm_squared, 25.0);
200 /// ```
201 #[inline]
202 #[must_use]
203 pub fn norm_squared(&self) -> f64 {
204 self.scalar * self.scalar + self.vector.dot(&self.vector)
205 }
206
207 /// The norm of the quaternion.
208 /// ```math
209 /// |\mathbf{q}|
210 /// ```
211 ///
212 /// # Example
213 /// ```
214 /// use hoomd_vector::Quaternion;
215 ///
216 /// let q = Quaternion::from([3.0, 0.0, 4.0, 0.0]);
217 /// let norm = q.norm();
218 /// assert_eq!(norm, 5.0);
219 /// ```
220 #[inline]
221 #[must_use]
222 pub fn norm(&self) -> f64 {
223 self.norm_squared().sqrt()
224 }
225
226 /// Construct the conjugate of this quaternion.
227 /// ```math
228 /// \mathbf{q}^* = (s, -\vec{v})
229 /// ```
230 ///
231 /// # Example
232 /// ```
233 /// use hoomd_vector::Quaternion;
234 ///
235 /// let q = Quaternion::from([1.0, 2.0, 3.0, 4.0]);
236 /// let q_star = q.conjugate();
237 /// assert_eq!(q_star, [1.0, -2.0, -3.0, -4.0].into());
238 /// ```
239 #[inline]
240 #[must_use]
241 pub fn conjugate(self) -> Self {
242 Self {
243 scalar: self.scalar,
244 vector: -self.vector,
245 }
246 }
247
248 /// Create a [`Versor`] by normalizing the given quaternion.
249 ///
250 /// ```math
251 /// \mathbf{v} = \frac{\mathbf{q}}{|\mathbf{q}|}
252 /// ```
253 ///
254 /// # Example
255 ///
256 /// ```
257 /// use hoomd_vector::{Quaternion, Versor};
258 ///
259 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
260 /// let q = Quaternion::from([3.0, 0.0, 0.0, 4.0]);
261 /// let v = q.to_versor()?;
262 /// assert_eq!(*v.get(), [3.0 / 5.0, 0.0, 0.0, 4.0 / 5.0].into());
263 /// # Ok(())
264 /// # }
265 /// ```
266 ///
267 /// # Errors
268 ///
269 /// [`Error::InvalidQuaternionMagnitude`] when `self` is the 0 quaternion.
270 #[inline]
271 pub fn to_versor(self) -> Result<Versor, Error> {
272 let mag = self.norm();
273 if mag == 0.0 {
274 Err(Error::InvalidQuaternionMagnitude)
275 } else {
276 Ok(Versor(self / mag))
277 }
278 }
279
280 /// Create a [`Versor`] by normalizing the given quaternion.
281 ///
282 /// ```math
283 /// \mathbf{v} = \frac{\mathbf{q}}{|\mathbf{q}|}
284 /// ```
285 ///
286 /// # Example
287 ///
288 /// ```
289 /// use hoomd_vector::{Quaternion, Versor};
290 ///
291 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
292 /// let q = Quaternion::from([3.0, 0.0, 0.0, 4.0]);
293 /// let v = q.to_versor_unchecked();
294 /// assert_eq!(*v.get(), [3.0 / 5.0, 0.0, 0.0, 4.0 / 5.0].into());
295 /// # Ok(())
296 /// # }
297 /// ```
298 ///
299 /// # Panics
300 ///
301 /// Divide by 0 when `self` is the 0 quaternion.
302 #[inline]
303 #[must_use]
304 pub fn to_versor_unchecked(self) -> Versor {
305 Versor(self / self.norm())
306 }
307}
308
309impl From<[f64; 4]> for Quaternion {
310 /// Construct a [`Quaternion`] from 4 values.
311 ///
312 /// The first value is the real part. The 2nd through 4th are the complex vector part:
313 /// `[scalar, vector_0, vector_1, vector_2]`.
314 ///
315 /// # Example
316 /// ```
317 /// use hoomd_vector::Quaternion;
318 ///
319 /// let q = Quaternion::from([1.0, 2.0, 3.0, 4.0]);
320 /// assert_eq!(q.scalar, 1.0);
321 /// assert_eq!(q.vector, [2.0, 3.0, 4.0].into());
322 /// ```
323 #[inline]
324 fn from(value: [f64; 4]) -> Self {
325 Self {
326 scalar: value[0],
327 vector: [value[1], value[2], value[3]].into(),
328 }
329 }
330}
331
332impl fmt::Display for Quaternion {
333 /// Format a [`Quaternion`] as `[{s}, [{v[0]}, {v[1]}, {v[2]}]]`.
334 #[inline]
335 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
336 write!(f, "[{}, {}]", self.scalar, self.vector)
337 }
338}
339
340impl Add for Quaternion {
341 type Output = Self;
342
343 #[inline]
344 fn add(self, rhs: Self) -> Self {
345 Self {
346 scalar: self.scalar + rhs.scalar,
347 vector: self.vector + rhs.vector,
348 }
349 }
350}
351
352impl AddAssign for Quaternion {
353 #[inline]
354 fn add_assign(&mut self, rhs: Self) {
355 self.scalar += rhs.scalar;
356 self.vector += rhs.vector;
357 }
358}
359
360impl Div<f64> for Quaternion {
361 type Output = Self;
362
363 #[inline]
364 fn div(self, rhs: f64) -> Self {
365 Self {
366 scalar: self.scalar / rhs,
367 vector: self.vector / rhs,
368 }
369 }
370}
371
372impl DivAssign<f64> for Quaternion {
373 #[inline]
374 fn div_assign(&mut self, rhs: f64) {
375 self.scalar /= rhs;
376 self.vector /= rhs;
377 }
378}
379
380impl Mul<f64> for Quaternion {
381 type Output = Self;
382
383 #[inline]
384 fn mul(self, rhs: f64) -> Self {
385 Self {
386 scalar: self.scalar * rhs,
387 vector: self.vector * rhs,
388 }
389 }
390}
391
392impl MulAssign<f64> for Quaternion {
393 #[inline]
394 fn mul_assign(&mut self, rhs: f64) {
395 self.scalar *= rhs;
396 self.vector *= rhs;
397 }
398}
399
400impl Mul<Quaternion> for Quaternion {
401 type Output = Self;
402
403 #[inline]
404 fn mul(self, rhs: Quaternion) -> Self {
405 Self {
406 scalar: (self.scalar * rhs.scalar - self.vector.dot(&rhs.vector)),
407 vector: (rhs.vector * self.scalar
408 + self.vector * rhs.scalar
409 + self.vector.cross(&rhs.vector)),
410 }
411 }
412}
413
414impl MulAssign<Quaternion> for Quaternion {
415 #[inline]
416 fn mul_assign(&mut self, rhs: Quaternion) {
417 let result = *self * rhs;
418 self.scalar = result.scalar;
419 self.vector = result.vector;
420 }
421}
422
423impl Sub for Quaternion {
424 type Output = Self;
425
426 #[inline]
427 fn sub(self, rhs: Self) -> Self {
428 Self {
429 scalar: self.scalar - rhs.scalar,
430 vector: self.vector - rhs.vector,
431 }
432 }
433}
434
435impl SubAssign for Quaternion {
436 #[inline]
437 fn sub_assign(&mut self, rhs: Self) {
438 self.scalar -= rhs.scalar;
439 self.vector -= rhs.vector;
440 }
441}
442
443/// A unit [`Quaternion`] that represents a 3D rotation.
444///
445/// [`Versor`] represents a 3D rotation with a **unit quaternion**. Rotation follows the Hamilton
446/// convention.
447///
448/// ## Constructing a [`Versor`]:
449///
450/// The default [`Versor`] is the identity:
451///
452/// ```
453/// use hoomd_vector::Versor;
454///
455/// let v = Versor::default();
456/// assert_eq!(*v.get(), [1.0, 0.0, 0.0, 0.0].into());
457/// ```
458///
459/// Create a [`Versor`] that rotates by an angle about an axis:
460/// ```
461/// use hoomd_vector::Versor;
462/// use std::f64::consts::PI;
463///
464/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
465/// let v = Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, PI / 2.0);
466/// assert_eq!(
467/// *v.get(),
468/// [(PI / 4.0).cos(), 0.0, (PI / 4.0).sin(), 0.0].into()
469/// );
470/// # Ok(())
471/// # }
472/// ```
473///
474/// Create a [`Versor`] by normalizing a [`Quaternion`]:
475/// ```
476/// use hoomd_vector::{Quaternion, Versor};
477///
478/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
479/// let q = Quaternion::from([3.0, 0.0, 0.0, 4.0]);
480/// let v = q.to_versor()?;
481/// assert_eq!(*v.get(), [3.0 / 5.0, 0.0, 0.0, 4.0 / 5.0].into());
482/// # Ok(())
483/// # }
484/// ```
485///
486/// Create a random [`Versor`]:
487/// ```
488/// use hoomd_vector::Versor;
489/// use rand::{RngExt, SeedableRng, rngs::StdRng};
490///
491/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
492/// let mut rng = StdRng::seed_from_u64(1);
493/// let v: Versor = rng.random();
494/// # Ok(())
495/// # }
496/// ```
497///
498/// ## Operations using [`Versor`]
499///
500/// Rotate a [`Cartesian<3>`] by a [`Versor`]:
501/// ```
502/// use approxim::assert_relative_eq;
503/// use hoomd_vector::{Cartesian, Rotate, Rotation, Versor};
504/// use std::f64::consts::PI;
505///
506/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
507/// let a = Cartesian::from([-1.0, 0.0, 0.0]);
508/// let v = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0);
509/// let b = v.rotate(&a);
510/// assert_relative_eq!(b, [0.0, -1.0, 0.0].into());
511/// # Ok(())
512/// # }
513/// ```
514///
515/// Combine two rotations together:
516/// ```
517/// use hoomd_vector::{Rotation, Versor};
518/// use std::f64::consts::PI;
519///
520/// # fn main() -> Result<(), Box<dyn std::error::Error>> {
521/// let a = Versor::from_axis_angle([1.0, 0.0, 1.0].try_into()?, PI / 2.0);
522/// let b = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 4.0);
523/// let c = a.combine(&b);
524/// # Ok(())
525/// # }
526/// ```
527#[derive(Clone, Copy, Debug, PartialEq, RelativeEq, Serialize, Deserialize)]
528pub struct Versor(Quaternion);
529
530impl Versor {
531 /// Take the dot product of the Versor as an element of $`\mathbb{R}^4`$.
532 #[inline]
533 fn dot_as_cartesian(&self, other: &Self) -> f64 {
534 self.get().scalar * other.get().scalar + self.get().vector.dot(&other.get().vector)
535 }
536 /// Create a [`Versor`] that rotates by an angle (in radians)
537 /// counterclockwise about an axis.
538 ///
539 /// # Example
540 ///
541 /// ```
542 /// use hoomd_vector::Versor;
543 /// use std::f64::consts::PI;
544 ///
545 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
546 /// let v = Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, PI / 2.0);
547 /// assert_eq!(
548 /// *v.get(),
549 /// [(PI / 4.0).cos(), 0.0, (PI / 4.0).sin(), 0.0].into()
550 /// );
551 /// # Ok(())
552 /// # }
553 /// ```
554 #[inline]
555 #[must_use]
556 pub fn from_axis_angle(axis: Unit<Cartesian<3>>, angle: f64) -> Self {
557 let Unit(axis_vector) = axis;
558
559 Versor(Quaternion {
560 scalar: (angle / 2.0).cos(),
561 vector: axis_vector * (angle / 2.0).sin(),
562 })
563 }
564
565 /// Normalize the versor.
566 ///
567 /// Nominally, all [`Versor`] instances retain a unit norm. Due to limited
568 /// floating point precision, this assumption may not hold after repeated
569 /// operations. Normalize versors when needed to correct this issue.
570 ///
571 /// # Example
572 ///
573 /// ```
574 /// use hoomd_vector::Versor;
575 /// use std::f64::consts::PI;
576 ///
577 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
578 /// let a = Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, PI / 2.0);
579 /// let b = a.normalized();
580 /// # Ok(())
581 /// # }
582 /// ```
583 #[inline]
584 #[must_use]
585 pub fn normalized(self) -> Self {
586 let Versor(q) = self;
587 let f = 1.0 / q.norm();
588 Self(Quaternion {
589 scalar: q.scalar * f,
590 vector: q.vector * f,
591 })
592 }
593
594 /// Get the unit quaternion.
595 #[inline]
596 #[must_use]
597 pub fn get(&self) -> &Quaternion {
598 &self.0
599 }
600
601 /// A metric quantifying the angle (in radians) of the spherical arc separating two Versors.
602 ///
603 /// $`d : \mathbb{H} \times \mathbb{H} \to \mathbb{R}^+, \quad d(q_0, q_1) = \arccos(|q_0 \cdot q_1|)`$
604 ///
605 /// This value always lies in the range $`[0, \pi]`$, and is symmetric: while there
606 /// are multiple arcs separating a pair of quaternions, this metric always chooses
607 /// the shortest.
608 #[inline]
609 #[must_use]
610 pub fn arc_distance(&self, other: &Self) -> f64 {
611 self.dot_as_cartesian(other).acos()
612 }
613 /// A fast metric on Versors representing elements of SO(3).
614 ///
615 /// $`d : \mathbb{H} \times \mathbb{H} \to \mathbb{R}^+, \quad d(q_0, q_1) = 1 - |q_0 \cdot q_1 |`$
616 ///
617 /// This has less geometric meaning than the [`arc_distance`](Versor::arc_distance) metric. However, it
618 /// is much faster while still obeying the triangle inequality and the axiom
619 /// $`d(q_0, q_1) = d(q_1, q_0)`$. This metric always lies in the range
620 /// $`[0, 1]`$, and is symmetric such that $`d(q, q)`$ = $`d(q, -q)`$.
621 #[inline]
622 #[must_use]
623 pub fn half_euclidean_norm_squared(&self, other: &Self) -> f64 {
624 1.0 - self.dot_as_cartesian(other)
625 }
626}
627
628impl From<Versor> for RotationMatrix<3> {
629 /// Construct a rotation matrix equivalent to this versor's rotation.
630 ///
631 /// When rotating many vectors by the same [`Versor`], improve performance
632 /// by converting to a matrix first and applying that matrix to the vectors.
633 ///
634 /// # Example
635 /// ```
636 /// use approxim::assert_relative_eq;
637 /// use hoomd_vector::{Cartesian, Rotate, RotationMatrix, Versor};
638 /// use std::f64::consts::PI;
639 ///
640 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
641 /// let a = Cartesian::from([-1.0, 0.0, 0.0]);
642 /// let v = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0);
643 ///
644 /// let matrix = RotationMatrix::from(v);
645 /// let b = matrix.rotate(&a);
646 /// assert_relative_eq!(b, [0.0, -1.0, 0.0].into());
647 /// # Ok(())
648 /// # }
649 /// ```
650 #[inline]
651 fn from(versor: Versor) -> RotationMatrix<3> {
652 let Versor(quaternion) = versor;
653 let a = quaternion.scalar;
654 let b = quaternion.vector[0];
655 let c = quaternion.vector[1];
656 let d = quaternion.vector[2];
657
658 RotationMatrix {
659 rows: [
660 [
661 a * a + b * b - c * c - d * d,
662 2.0 * b * c - 2.0 * a * d,
663 2.0 * b * d + 2.0 * a * c,
664 ]
665 .into(),
666 [
667 2.0 * b * c + 2.0 * a * d,
668 a * a - b * b + c * c - d * d,
669 2.0 * c * d - 2.0 * a * b,
670 ]
671 .into(),
672 [
673 2.0 * b * d - 2.0 * a * c,
674 2.0 * c * d + 2.0 * a * b,
675 a * a - b * b - c * c + d * d,
676 ]
677 .into(),
678 ],
679 }
680 }
681}
682
683impl Default for Versor {
684 /// Create an identity rotation.
685 ///
686 /// # Example
687 /// ```
688 /// use hoomd_vector::Versor;
689 ///
690 /// let v = Versor::default();
691 /// ```
692 #[inline]
693 fn default() -> Self {
694 Self(Quaternion {
695 scalar: 1.0,
696 vector: [0.0, 0.0, 0.0].into(),
697 })
698 }
699}
700
701impl Rotate<Cartesian<3>> for Versor {
702 type Matrix = RotationMatrix<3>;
703
704 /// Rotate a [`Cartesian<3>`] by a [`Versor`]
705 ///
706 /// ```math
707 /// \mathbf{q} \vec{a} \mathbf{q}^*
708 /// ```
709 ///
710 /// # Example
711 ///
712 /// ```
713 /// use approxim::assert_relative_eq;
714 /// use hoomd_vector::{Cartesian, Rotate, Rotation, Versor};
715 /// use std::f64::consts::PI;
716 ///
717 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
718 /// let a = Cartesian::from([-1.0, 0.0, 0.0]);
719 /// let v = Versor::from_axis_angle([0.0, 0.0, 1.0].try_into()?, PI / 2.0);
720 /// let b = v.rotate(&a);
721 /// assert_relative_eq!(b, [0.0, -1.0, 0.0].into());
722 /// # Ok(())
723 /// # }
724 /// ```
725 #[inline]
726 fn rotate(&self, vector: &Cartesian<3>) -> Cartesian<3> {
727 let Versor(quaternion) = self;
728
729 *vector
730 * (quaternion.scalar * quaternion.scalar - quaternion.vector.dot(&quaternion.vector))
731 + quaternion.vector.cross(vector) * (2.0 * quaternion.scalar)
732 + quaternion.vector * (2.0 * quaternion.vector.dot(vector))
733 }
734}
735
736impl Rotation for Versor {
737 /// Combine two rotations.
738 ///
739 /// The resulting versor is obtained by quaternion multiplication.
740 /// ```math
741 /// \mathbf{q}_{ab} = \mathbf{q}_a \mathbf{q}_b
742 /// ```
743 ///
744 /// # Example
745 ///
746 /// ```
747 /// use hoomd_vector::{Rotation, Versor};
748 ///
749 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
750 /// let q_a = Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, 1.5);
751 /// let q_b = Versor::from_axis_angle([1.0, 0.0, 0.0].try_into()?, 0.125);
752 /// let q_ab = q_a.combine(&q_b);
753 /// # Ok(())
754 /// # }
755 /// ```
756 #[inline]
757 fn combine(&self, other: &Self) -> Self {
758 let Versor(a) = self;
759 let Versor(b) = other;
760
761 Versor(a.mul(*b))
762 }
763
764 /// Create the identity [`Versor`]: [1, [0, 0, 0]]
765 ///
766 /// # Example
767 ///
768 /// ```
769 /// use hoomd_vector::{Rotation, Versor};
770 ///
771 /// let identity = Versor::identity();
772 /// ```
773 #[inline]
774 fn identity() -> Self {
775 Self::default()
776 }
777
778 /// Create a [`Versor`] that performs the inverse rotation of the given versor.
779 ///
780 /// ```math
781 /// \mathbf{q}^*
782 /// ```
783 ///
784 /// # Example
785 ///
786 /// ```
787 /// use hoomd_vector::{Rotation, Versor};
788 ///
789 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
790 /// let v = Versor::from_axis_angle([0.0, 1.0, 0.0].try_into()?, 1.5);
791 /// let v_star = v.inverted();
792 /// # Ok(())
793 /// # }
794 /// ```
795 #[inline]
796 fn inverted(self) -> Self {
797 let Versor(quaternion) = self;
798
799 Versor(quaternion.conjugate())
800 }
801}
802
803impl fmt::Display for Versor {
804 /// Format a [`Versor`] as `[{s}, [{v[0]}, {v[1]}, {v[2]}]]`.
805 #[inline]
806 fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
807 write!(f, "{}", self.0)
808 }
809}
810
811impl Distribution<Versor> for StandardUniform {
812 /// Sample a random [`Versor`] from the uniform distribution over all rotations.
813 ///
814 /// # Example
815 ///
816 /// ```
817 /// use hoomd_vector::Versor;
818 /// use rand::{RngExt, SeedableRng, rngs::StdRng};
819 ///
820 /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
821 /// let mut rng = StdRng::seed_from_u64(1);
822 /// let v: Versor = rng.random();
823 /// # Ok(())
824 /// # }
825 /// ```
826 #[inline]
827 fn sample<R: Rng + ?Sized>(&self, rng: &mut R) -> Versor {
828 // See method 19 from: https://extremelearning.com.au/how-to-generate-uniformly-random-points-on-n-spheres-and-n-balls/
829 let scalar = rng.sample::<f64, _>(StandardNormal);
830 let vector = Cartesian::<3>::from(std::array::from_fn(|_| rng.sample(StandardNormal)));
831
832 let norm = (vector.norm_squared() + (scalar * scalar)).sqrt();
833
834 Versor(Quaternion {
835 scalar: scalar / norm,
836 vector: vector / norm,
837 })
838 }
839}
840
841#[cfg(test)]
842mod tests {
843 use super::*;
844 use approxim::{assert_abs_diff_eq, assert_relative_eq};
845 use rand::{SeedableRng, rngs::StdRng};
846 use rstest::*;
847 use std::f64::consts::PI;
848
849 mod quaternion {
850 use super::*;
851
852 #[test]
853 fn from_array() {
854 let q = Quaternion::from([2.0, -3.0, 4.0, 7.0]);
855 assert!(q.scalar == 2.0);
856 assert!(q.vector == [-3.0, 4.0, 7.0].into());
857 }
858
859 #[test]
860 fn norm() {
861 let q = Quaternion::from([1.0, 4.0, -3.0, -2.0]);
862 assert_eq!(q.norm_squared(), 30.0);
863 assert_eq!(q.norm(), 30.0_f64.sqrt());
864 }
865
866 #[test]
867 fn conjugate() {
868 let q1 = Quaternion::from([1.0, -2.0, 4.0, -0.5]);
869 let q2 = q1.conjugate();
870 assert_eq!(q2, [1.0, 2.0, -4.0, 0.5].into());
871 assert_relative_eq!(q2 * q1, [q2.norm() * q1.norm(), 0.0, 0.0, 0.0].into());
872 }
873
874 #[test]
875 fn to_versor() {
876 let q = Quaternion::from([5.0, 3.0, -1.0, 1.0]);
877
878 assert_relative_eq!(
879 q.to_versor()
880 .expect("hard-coded quatnernion should be non zero"),
881 Versor(Quaternion {
882 scalar: 5.0 / 6.0,
883 vector: [3.0 / 6.0, -1.0 / 6.0, 1.0 / 6.0].into()
884 })
885 );
886
887 assert_relative_eq!(
888 q.to_versor_unchecked(),
889 Versor(Quaternion {
890 scalar: 5.0 / 6.0,
891 vector: [3.0 / 6.0, -1.0 / 6.0, 1.0 / 6.0].into()
892 })
893 );
894
895 let zero = Quaternion::from([0.0, 0.0, 0.0, 0.0]);
896 assert!(matches!(
897 zero.to_versor(),
898 Err(Error::InvalidQuaternionMagnitude)
899 ));
900 }
901
902 #[test]
903 fn ops() {
904 let a = Quaternion::from([1.0, -2.0, 6.0, -4.0]);
905 let b = Quaternion::from([-2.0, 6.0, 4.0, 1.0]);
906
907 // +, +=
908 assert_eq!(a + b, [-1.0, 4.0, 10.0, -3.0].into());
909 let mut c = a;
910 c += b;
911 assert_eq!(c, [-1.0, 4.0, 10.0, -3.0].into());
912
913 // -, -=
914 assert_eq!(a - b, [3.0, -8.0, 2.0, -5.0].into());
915 let mut c = a;
916 c -= b;
917 assert_eq!(c, [3.0, -8.0, 2.0, -5.0].into());
918
919 // Scalar * and /
920 assert_eq!(a * 2.0, [2.0, -4.0, 12.0, -8.0].into());
921 let mut c = a;
922 c *= 2.0;
923 assert_eq!(c, [2.0, -4.0, 12.0, -8.0].into());
924
925 assert_eq!(a / 2.0, [0.5, -1.0, 3.0, -2.0].into());
926 let mut c = a;
927 c /= 2.0;
928 assert_eq!(c, [0.5, -1.0, 3.0, -2.0].into());
929
930 // Quaternion multiplication
931 assert_eq!(a * b, [-10.0, 32.0, -30.0, -35.0].into());
932 let mut c = a;
933 c *= b;
934 assert_eq!(c, [-10.0, 32.0, -30.0, -35.0].into());
935 }
936
937 #[test]
938 fn display() {
939 let q = Quaternion {
940 scalar: 0.5,
941 vector: [0.125, -0.875, 2.125].into(),
942 };
943 let s = format!("{q}");
944 assert_eq!(s, "[0.5, [0.125, -0.875, 2.125]]");
945 }
946 }
947
948 mod versor {
949 use super::*;
950 #[test]
951 fn default() {
952 let a = Versor::default();
953 assert!(a.get() == &[1.0, 0.0, 0.0, 0.0].into());
954 }
955
956 #[test]
957 fn identity() {
958 let a = Versor::identity();
959 assert!(a.get() == &[1.0, 0.0, 0.0, 0.0].into());
960 }
961
962 #[rstest(
963 theta => [0.0, PI / 2.0, 1e-12 * PI, -3.0, 12345.6],
964 axis => [[1.0, 0.0, 0.0].try_into().expect("hard-coded vector should have non-zero length"), [1.0, -1.0, 1.0].try_into().expect("hard-coded vector should have non-zero length")],
965 )]
966 fn from_axis_angle(theta: f64, axis: Unit<Cartesian<3>>) {
967 let Unit(axis_vector) = axis;
968
969 let Versor(q) = Versor::from_axis_angle(axis, theta);
970 assert_relative_eq!(q.scalar, (theta / 2.0).cos());
971 assert_relative_eq!(q.vector, axis_vector * (theta / 2.0).sin());
972 }
973
974 #[rstest(
975 theta_1 => [0.0, PI / 2.0, -3.0],
976 theta_2 => [-0.0, -PI / 3.0, PI, 2.0 * PI]
977 )]
978 fn combine_same_axis(theta_1: f64, theta_2: f64) {
979 let axis = [1.0, 0.0, 0.0]
980 .try_into()
981 .expect("hard-coded vector should have non-zero length");
982 let Unit(axis_vector) = axis;
983
984 let a = Versor::from_axis_angle(axis, theta_1);
985 let b = Versor::from_axis_angle(axis, theta_2);
986 let c = a.combine(&b);
987
988 let theta = theta_1 + theta_2;
989 let Versor(q) = c;
990 assert_relative_eq!(q.scalar, (theta / 2.0).cos());
991 assert_relative_eq!(q.vector, axis_vector * (theta / 2.0).sin());
992 }
993
994 fn validate_rotations<R: Rotate<Cartesian<3>>>(z_pi_2: &R, y_pi_4: &R) {
995 assert_relative_eq!(
996 z_pi_2.rotate(&[0.0, 0.0, 1.0].into()),
997 [0.0, 0.0, 1.0].into()
998 );
999 assert_relative_eq!(
1000 z_pi_2.rotate(&[1.0, 0.0, 4.25].into()),
1001 [0.0, 1.0, 4.25].into()
1002 );
1003 assert_relative_eq!(
1004 z_pi_2.rotate(&[0.0, 1.0, -8.75].into()),
1005 [-1.0, 0.0, -8.75].into()
1006 );
1007
1008 let sqrt_2_2 = 2.0_f64.sqrt() / 2.0;
1009 assert_relative_eq!(
1010 y_pi_4.rotate(&[0.0, -10.0, 0.0].into()),
1011 [0.0, -10.0, 0.0].into()
1012 );
1013 assert_relative_eq!(
1014 y_pi_4.rotate(&[1.0, -15.0, 0.0].into()),
1015 [sqrt_2_2, -15.0, -sqrt_2_2].into()
1016 );
1017 assert_relative_eq!(
1018 y_pi_4.rotate(&[sqrt_2_2, -15.0, -sqrt_2_2].into()),
1019 [0.0, -15.0, -1.0].into()
1020 );
1021 }
1022
1023 #[test]
1024 fn rotate() {
1025 let z_pi_2 = Versor::from_axis_angle(
1026 [0.0, 0.0, 1.0]
1027 .try_into()
1028 .expect("hard-coded vector should have non-zero length"),
1029 PI / 2.0,
1030 );
1031 let y_pi_4 = Versor::from_axis_angle(
1032 [0.0, 1.0, 0.0]
1033 .try_into()
1034 .expect("hard-coded vector should have non-zero length"),
1035 PI / 4.0,
1036 );
1037
1038 validate_rotations(&z_pi_2, &y_pi_4);
1039 }
1040
1041 #[test]
1042 fn precompute() {
1043 let z_pi_2 = RotationMatrix::from(Versor::from_axis_angle(
1044 [0.0, 0.0, 1.0]
1045 .try_into()
1046 .expect("hard-coded vector should have non-zero length"),
1047 PI / 2.0,
1048 ));
1049 let y_pi_4 = RotationMatrix::from(Versor::from_axis_angle(
1050 [0.0, 1.0, 0.0]
1051 .try_into()
1052 .expect("hard-coded vector should have non-zero length"),
1053 PI / 4.0,
1054 ));
1055
1056 validate_rotations(&z_pi_2, &y_pi_4);
1057 }
1058
1059 #[test]
1060 fn combine_different_axis() {
1061 let a = Versor::from_axis_angle(
1062 [1.0, 0.0, 0.0]
1063 .try_into()
1064 .expect("hard-coded vector should have non-zero length"),
1065 PI / 4.0,
1066 );
1067 let b = Versor::from_axis_angle(
1068 [0.0, 0.0, 1.0]
1069 .try_into()
1070 .expect("hard-coded vector should have non-zero length"),
1071 PI / 2.0,
1072 );
1073
1074 let q = a.combine(&b);
1075 let v = q.rotate(&[1.0, 0.0, 0.0].into());
1076 assert_relative_eq!(v, [0.0, 2.0_f64.sqrt() / 2.0, 2.0_f64.sqrt() / 2.0].into());
1077 }
1078
1079 #[rstest(theta => [0.0, 1.0, 2.125])]
1080 fn inverted(theta: f64) {
1081 let q1 = Versor::from_axis_angle(
1082 [1.0, 0.5, -2.0]
1083 .try_into()
1084 .expect("hard-coded vector should have non-zero length"),
1085 theta,
1086 );
1087 let q2 = q1.inverted();
1088 assert_relative_eq!(q1.combine(&q2), Versor::identity());
1089 }
1090
1091 #[test]
1092 fn display() {
1093 let v = Versor(Quaternion {
1094 scalar: 0.5,
1095 vector: [0.125, -0.875, 2.125].into(),
1096 });
1097 let s = format!("{v}");
1098 assert_eq!(s, "[0.5, [0.125, -0.875, 2.125]]");
1099 }
1100
1101 #[test]
1102 fn normalized() {
1103 let v = Versor(Quaternion {
1104 scalar: 5.0,
1105 vector: [3.0, -1.0, 1.0].into(),
1106 });
1107 assert_relative_eq!(
1108 v.normalized(),
1109 Versor(Quaternion {
1110 scalar: 5.0 / 6.0,
1111 vector: [3.0 / 6.0, -1.0 / 6.0, 1.0 / 6.0].into()
1112 })
1113 );
1114 }
1115
1116 #[test]
1117 fn random() {
1118 const CHECK_VECTORS: [Cartesian<3>; 3] = [
1119 Cartesian {
1120 coordinates: [1.0, 0.0, 0.0],
1121 },
1122 Cartesian {
1123 coordinates: [0.0, 1.0, 0.0],
1124 },
1125 Cartesian {
1126 coordinates: [1.0, 0.0, 1.0],
1127 },
1128 ];
1129
1130 // Perform basic checks on random versors.
1131 // 1) Ensure that each randomly generated versor is unit.
1132 // 2) Check that the result of rotating a reference vector by random versors does not
1133 // point in any special direction. The average dot product should be close to 0.
1134 let samples: u32 = 20_000;
1135
1136 let reference = Cartesian::from([1.0, 0.0, 0.0]);
1137 let mut dot_sums = [0.0; CHECK_VECTORS.len()];
1138
1139 let mut rng = StdRng::seed_from_u64(1);
1140
1141 for _ in 0..samples {
1142 let q: Versor = rng.random();
1143 assert_relative_eq!(q.get().norm_squared(), 1.0, max_relative = 1e-15);
1144
1145 let v = q.rotate(&reference);
1146 for i in 0..CHECK_VECTORS.len() {
1147 dot_sums[i] += v.dot(&CHECK_VECTORS[i]);
1148 }
1149 }
1150
1151 for dot_sum in dot_sums {
1152 assert_abs_diff_eq!(dot_sum / f64::from(samples), 0.0, epsilon = 0.01);
1153 }
1154 }
1155 }
1156}