From 8e9b4a90fd5472ae0026b0dc98c92a4be34b9674 Mon Sep 17 00:00:00 2001 From: Rouven Spreckels Date: Mon, 11 Mar 2024 18:22:32 +0100 Subject: [PATCH] Refactor. --- Cargo.toml | 3 +- README.md | 4 +- RELEASES.md | 8 + src/lib.rs | 597 ++++++++++++++++++++++++++-------------------------- 4 files changed, 310 insertions(+), 302 deletions(-) diff --git a/Cargo.toml b/Cargo.toml index a998ee6..5abe5e0 100644 --- a/Cargo.toml +++ b/Cargo.toml @@ -1,6 +1,6 @@ [package] name = "nalgebra-spacetime" -version = "0.4.0" +version = "0.5.0" authors = ["Rouven Spreckels "] edition = "2021" description = "Spacetime Extension for nalgebra" @@ -34,4 +34,3 @@ rustdoc-args = [ "--html-in-header", "katex.html" ] [dependencies] approx = { version = "0.5.1", default-features = false } nalgebra = { version = "0.32.4", features = ["rand"] } -num-traits = "0.2.18" diff --git a/README.md b/README.md index 653ff6d..c9cf0de 100644 --- a/README.md +++ b/README.md @@ -20,11 +20,11 @@ Spacetime Extension for [nalgebra] # Present Features - * Minkowski space as special case of `LorentzianN` space. + * Minkowski space as special case of n-dimensional `Lorentzian` space. * Raising/Lowering tensor indices: `dual()`/`r_dual()`/`c_dual()`. * Metric contraction of degree-1/degree-2 tensors: `contr()`/`scalar()`. * Spacetime `interval()` with `LightCone` depiction. - * Inertial `FrameN` of reference holding boost parameters. + * Inertial `OFrame` of reference holding boost parameters. * Lorentz boost as `new_boost()` matrix. * Direct Lorentz `boost()` to `compose()` velocities. * Wigner `rotation()` and `axis_angle()` between to-be-composed boosts. diff --git a/RELEASES.md b/RELEASES.md index 6ce9cc4..1c25b0a 100644 --- a/RELEASES.md +++ b/RELEASES.md @@ -1,3 +1,11 @@ +# Version 0.5.0 (2024-03-11) + + * Replace `N: SimdRealField + Signed + Real` with `T: RealField`. + * Rename `LorentzianN` to `Lorentzian`. + * Rename `FrameN` to `OFrame`. + * Rename `MomentumN` to `OMomentum`. + * Make `OFrame::rotation()` robust close to identity. + # Version 0.4.0 (2024-02-27) * Make Wigner rotation n-dimensional. diff --git a/src/lib.rs b/src/lib.rs index 4053201..c80ff5b 100644 --- a/src/lib.rs +++ b/src/lib.rs @@ -4,22 +4,22 @@ //! //! # Present Features //! -//! * Minkowski space as special case of [`LorentzianN`] space. +//! * Minkowski space as special case of $n$-dimensional [`Lorentzian`] space. //! * Raising/Lowering tensor indices: -//! [`LorentzianN::dual`]/[`LorentzianN::r_dual`]/[`LorentzianN::c_dual`]. +//! [`Lorentzian::dual`]/[`Lorentzian::r_dual`]/[`Lorentzian::c_dual`]. //! * Metric contraction of degree-1/degree-2 tensors: -//! [`LorentzianN::contr`]/[`LorentzianN::scalar`]. -//! * Spacetime [`LorentzianN::interval`] with [`LightCone`] depiction. -//! * Inertial [`FrameN`] of reference holding boost parameters. -//! * Lorentz boost as [`LorentzianN::new_boost`] matrix. -//! * Direct Lorentz [`LorentzianN::boost`] to [`FrameN::compose`] velocities. -//! * Wigner [`FrameN::rotation`] and [`FrameN::axis_angle`] between to-be-composed boosts. +//! [`Lorentzian::contr`]/[`Lorentzian::scalar`]. +//! * Spacetime [`Lorentzian::interval`] with [`LightCone`] depiction. +//! * Inertial [`OFrame`] of reference holding boost parameters. +//! * Lorentz boost as [`Lorentzian::new_boost`] matrix. +//! * Direct Lorentz [`Lorentzian::boost`] to [`OFrame::compose`] velocities. +//! * Wigner [`OFrame::rotation`] and [`OFrame::axis_angle`] between to-be-composed boosts. //! //! # Future Features //! //! * `Event4`/`Velocity4`/`Momentum4`/`...` equivalents of `Point4`/`...`. //! * Categorize `Rotation4`/`PureBoost4`/`...` as `Boost4`/`...`. -//! * Wigner [`FrameN::rotation`] and [`FrameN::axis_angle`] of an already-composed `Boost4`. +//! * Wigner [`OFrame::rotation`] and [`OFrame::axis_angle`] of an already-composed `Boost4`. //! * Distinguish pre/post-rotation and active/passive `Boost4` compositions. #![forbid(unsafe_code)] @@ -29,7 +29,6 @@ pub use approx; pub use nalgebra; -pub use num_traits; use approx::{abs_diff_eq, AbsDiffEq}; use nalgebra::{ @@ -41,11 +40,15 @@ use nalgebra::{ }, storage::{Owned, RawStorage, RawStorageMut, Storage, StorageMut}, DefaultAllocator, Dim, DimName, DimNameDiff, DimNameSub, Matrix, MatrixView, MatrixViewMut, - OMatrix, OVector, Scalar, SimdRealField, Unit, VectorView, + OMatrix, OVector, RealField, Scalar, Unit, VectorView, }; -use num_traits::{real::Real, sign::Signed}; use std::ops::{Add, Neg, Sub}; +#[inline] +fn neg(n: &mut T) { + *n = -n.clone(); +} + /// Extension for $n$-dimensional Lorentzian space $\R^{-,+} = \R^{1,n-1}$ with /// metric signature in spacelike sign convention. /// @@ -53,9 +56,9 @@ use std::ops::{Add, Neg, Sub}; /// /// A statically sized column-major matrix whose `R` rows and `C` columns /// coincide with degree-1/degree-2 tensor indices. -pub trait LorentzianN +pub trait Lorentzian where - N: Scalar, + T: Scalar, R: DimName, C: DimName, { @@ -88,7 +91,7 @@ where /// ``` /// use approx::assert_ulps_eq; /// use nalgebra::{Matrix4, Vector4}; - /// use nalgebra_spacetime::LorentzianN; + /// use nalgebra_spacetime::Lorentzian; /// /// let eta = Matrix4::::metric(); /// let sc = Vector4::new(-1.0, 1.0, 1.0, 1.0); @@ -148,26 +151,26 @@ where /// Equals `self.c_dual() * rhs`, the metric contraction of its *column* index with `rhs`'s /// *row* index. #[must_use] - fn contr(&self, rhs: &Matrix) -> OMatrix + fn contr(&self, rhs: &Matrix) -> OMatrix where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: AreMultipliable, - DefaultAllocator: Allocator; + DefaultAllocator: Allocator; /// Same as [`Self::contr`] but with transposed tensor indices. /// /// Equals `self.r_dual().tr_mul(rhs)`, the metric contraction of its *transposed row* index /// with `rhs`'s *row* index. #[must_use] - fn tr_contr(&self, rhs: &Matrix) -> OMatrix + fn tr_contr(&self, rhs: &Matrix) -> OMatrix where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: SameNumberOfRows, - DefaultAllocator: Allocator; + DefaultAllocator: Allocator; /// Lorentzian inner product of degree-1/degree-2 tensors. /// @@ -183,22 +186,22 @@ where /// * Lorentz scalar, invariant under Lorentz transformations, or /// * spacetime interval between two events, see [`Self::interval`]. #[must_use] - fn scalar(&self, rhs: &Matrix) -> N + fn scalar(&self, rhs: &Matrix) -> T where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: DimEq + DimEq; /// Same as [`Self::scalar`] but with transposed tensor indices. /// /// Equals `self.dual().tr_dot(rhs)`. #[must_use] - fn tr_scalar(&self, rhs: &Matrix) -> N + fn tr_scalar(&self, rhs: &Matrix) -> T where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: DimEq + DimEq; /// Lorentzian norm of timelike degree-1/degree-2 tensors. @@ -207,7 +210,7 @@ where /// /// If spacelike, returns *NaN* or panics if `N` doesn't support it. #[must_use] - fn timelike_norm(&self) -> N; + fn timelike_norm(&self) -> T; /// Lorentzian norm of spacelike degree-1/degree-2 tensors. /// @@ -215,18 +218,18 @@ where /// /// If timelike, returns *NaN* or panics if `N` doesn't support it. #[must_use] - fn spacelike_norm(&self) -> N; + fn spacelike_norm(&self) -> T; /// Spacetime interval between two events and region of `self`'s light cone. /// /// Same as [`Self::interval_fn`] but with [`AbsDiffEq::default_epsilon`] as in: /// - /// * `is_present = |time| abs_diff_eq!(time, N::zero())`, and - /// * `is_lightlike = |interval| abs_diff_eq!(interval, N::zero())`. + /// * `is_present = |time| abs_diff_eq!(time, T::zero())`, and + /// * `is_lightlike = |interval| abs_diff_eq!(interval, T::zero())`. #[must_use] - fn interval(&self, rhs: &Self) -> (N, LightCone) + fn interval(&self, rhs: &Self) -> (T, LightCone) where - N: AbsDiffEq, + T: AbsDiffEq, ShapeConstraint: DimEq; /// Spacetime interval between two events and region of `self`'s light cone. @@ -234,7 +237,7 @@ where /// Equals `(rhs - self).scalar(&(rhs - self))` where `self` is subtracted from `rhs` to depict /// `rhs` in `self`'s light cone. /// - /// Requires you to approximate when `N` equals `N::zero()` via: + /// Requires you to approximate when `N` equals `T::zero()` via: /// /// * `is_present` for the time component of `rhs - self`, and /// * `is_lightlike` for the interval. @@ -243,11 +246,11 @@ where /// /// See [`Self::interval`] for using defaults and [`approx`] for further details. #[must_use] - fn interval_fn(&self, rhs: &Self, is_present: P, is_lightlike: L) -> (N, LightCone) + fn interval_fn(&self, rhs: &Self, is_present: P, is_lightlike: L) -> (T, LightCone) where ShapeConstraint: DimEq, - P: Fn(N) -> bool, - L: Fn(N) -> bool; + P: Fn(&T) -> bool, + L: Fn(&T) -> bool; /// $ /// \gdef \uk {\hat u \cdot \vec K} @@ -288,7 +291,7 @@ where /// ``` /// use approx::assert_ulps_eq; /// use nalgebra::{Matrix4, Vector3, Vector4}; - /// use nalgebra_spacetime::{Frame4, LorentzianN}; + /// use nalgebra_spacetime::{Frame4, Lorentzian}; /// /// let event = Vector4::new_random(); /// let frame = Frame4::from_beta(Vector3::new(0.3, -0.4, 0.6)); @@ -296,11 +299,11 @@ where /// assert_ulps_eq!(boost * event, event.boost(&frame), epsilon = 1e-14); /// ``` #[must_use] - fn new_boost(frame: &FrameN) -> Self + fn new_boost(frame: &OFrame) -> Self where D: DimNameSub, ShapeConstraint: AreMultipliable + DimEq, - DefaultAllocator: Allocator>; + DefaultAllocator: Allocator>; /// Boosts this degree-1 tensor $x^\mu$ to inertial `frame` of reference along $\hat u$ with /// $\zeta$. @@ -308,7 +311,7 @@ where /// ``` /// use approx::assert_ulps_eq; /// use nalgebra::{Vector3, Vector4}; - /// use nalgebra_spacetime::{Frame4, LorentzianN}; + /// use nalgebra_spacetime::{Frame4, Lorentzian}; /// /// let muon_lifetime_at_rest = Vector4::new(2.2e-6, 0.0, 0.0, 0.0); /// let muon_frame = Frame4::from_axis_beta(Vector3::z_axis(), 0.9952); @@ -319,7 +322,7 @@ where /// /// See `boost_mut()` for further details. #[must_use] - fn boost(&self, frame: &FrameN) -> Self + fn boost(&self, frame: &OFrame) -> Self where R: DimNameSub, D: DimNameSub, @@ -327,7 +330,7 @@ where + SameNumberOfColumns + DimEq<>::Output, >::Output> + SameNumberOfRows<>::Output, >::Output>, - DefaultAllocator: Allocator>; + DefaultAllocator: Allocator>; /// $ /// \gdef \xu {(\vec x \cdot \hat u)} @@ -348,7 +351,7 @@ where /// ``` /// use approx::assert_ulps_eq; /// use nalgebra::Vector4; - /// use nalgebra_spacetime::LorentzianN; + /// use nalgebra_spacetime::Lorentzian; /// /// // Arbitrary timelike four-momentum. /// let mut momentum = Vector4::new(24.3, 5.22, 16.8, 9.35); @@ -366,7 +369,7 @@ where /// // Verify boosting four-momentum to center-of-momentum frame. /// assert_ulps_eq!(momentum, mass_at_rest, epsilon = 1e-14); /// ``` - fn boost_mut(&mut self, frame: &FrameN) + fn boost_mut(&mut self, frame: &OFrame) where R: DimNameSub, D: DimNameSub, @@ -374,43 +377,43 @@ where + SameNumberOfColumns + DimEq<>::Output, >::Output> + SameNumberOfRows<>::Output, >::Output>, - DefaultAllocator: Allocator>; + DefaultAllocator: Allocator>; /// Velocity $u^\mu$ of inertial `frame` of reference. #[must_use] - fn new_velocity(frame: &FrameN) -> OVector + fn new_velocity(frame: &OFrame) -> OVector where D: DimNameSub, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, - DefaultAllocator: Allocator> + Allocator; + DefaultAllocator: Allocator> + Allocator; /// Inertial frame of reference of this velocity $u^\mu$. #[must_use] - fn frame(&self) -> FrameN + fn frame(&self) -> OFrame where R: DimNameSub, ShapeConstraint: SameNumberOfColumns, - DefaultAllocator: Allocator>; + DefaultAllocator: Allocator>; /// From `temporal` and `spatial` spacetime split. #[must_use] - fn from_split(temporal: &N, spatial: &OVector>) -> OVector + fn from_split(temporal: &T, spatial: &OVector>) -> OVector where D: DimNameSub, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, - DefaultAllocator: Allocator> + Allocator, - >::Buffer: - StorageMut; + DefaultAllocator: Allocator> + Allocator, + >::Buffer: + StorageMut; /// Spacetime split into [`Self::temporal`] and [`Self::spatial`]. #[must_use] - fn split(&self) -> (&N, MatrixView, C, U1, R>) + fn split(&self) -> (&T, MatrixView, C, U1, R>) where R: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator, - >::Buffer: - Storage; + DefaultAllocator: Allocator, + >::Buffer: + Storage; /// Mutable spacetime split into [`Self::temporal_mut`] and /// [`Self::spatial_mut`]. @@ -418,7 +421,7 @@ where /// ``` /// use approx::assert_ulps_eq; /// use nalgebra::Vector4; - /// use nalgebra_spacetime::LorentzianN; + /// use nalgebra_spacetime::Lorentzian; /// /// let mut spacetime = Vector4::new(1.0, 2.0, 3.0, 4.0); /// let (temporal, mut spatial) = spacetime.split_mut(); @@ -428,53 +431,53 @@ where /// spatial[2] += 4.0; /// assert_ulps_eq!(spacetime, Vector4::new(2.0, 4.0, 6.0, 8.0)); /// ``` - fn split_mut(&mut self) -> (&mut N, MatrixViewMut, C>) + fn split_mut(&mut self) -> (&mut T, MatrixViewMut, C>) where R: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator, - >::Buffer: - Storage; + DefaultAllocator: Allocator, + >::Buffer: + Storage; /// Temporal component. #[must_use] - fn temporal(&self) -> &N + fn temporal(&self) -> &T where R: DimNameSub, ShapeConstraint: DimEq; /// Mutable temporal component. - fn temporal_mut(&mut self) -> &mut N + fn temporal_mut(&mut self) -> &mut T where R: DimNameSub, ShapeConstraint: DimEq; /// Spatial components. #[must_use] - fn spatial(&self) -> MatrixView, C, U1, R> + fn spatial(&self) -> MatrixView, C, U1, R> where R: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator, - >::Buffer: - Storage; + DefaultAllocator: Allocator, + >::Buffer: + Storage; /// Mutable spatial components. - fn spatial_mut(&mut self) -> MatrixViewMut, C, U1, R> + fn spatial_mut(&mut self) -> MatrixViewMut, C, U1, R> where R: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator, - >::Buffer: - StorageMut; + DefaultAllocator: Allocator, + >::Buffer: + StorageMut; } -impl LorentzianN for OMatrix +impl Lorentzian for OMatrix where - N: SimdRealField + Signed + Real, + T: RealField, R: DimName + DimNameSub, - C: DimName, - DefaultAllocator: Allocator, + C: DimName + DimNameSub, + DefaultAllocator: Allocator, { #[inline] fn metric() -> Self @@ -541,92 +544,92 @@ where } #[inline] - fn contr(&self, rhs: &Matrix) -> OMatrix + fn contr(&self, rhs: &Matrix) -> OMatrix where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: AreMultipliable, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { self.c_dual() * rhs } #[inline] - fn tr_contr(&self, rhs: &Matrix) -> OMatrix + fn tr_contr(&self, rhs: &Matrix) -> OMatrix where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: SameNumberOfRows, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { self.r_dual().tr_mul(rhs) } #[inline] - fn scalar(&self, rhs: &Matrix) -> N + fn scalar(&self, rhs: &Matrix) -> T where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: DimEq + DimEq, { self.dual().dot(rhs) } #[inline] - fn tr_scalar(&self, rhs: &Matrix) -> N + fn tr_scalar(&self, rhs: &Matrix) -> T where R2: Dim, C2: Dim, - SB: Storage, + SB: Storage, ShapeConstraint: DimEq + DimEq, { self.dual().tr_dot(rhs) } #[inline] - fn timelike_norm(&self) -> N { + fn timelike_norm(&self) -> T { self.scalar(self).neg().sqrt() } #[inline] - fn spacelike_norm(&self) -> N { + fn spacelike_norm(&self) -> T { self.scalar(self).sqrt() } #[inline] - fn interval(&self, rhs: &Self) -> (N, LightCone) + fn interval(&self, rhs: &Self) -> (T, LightCone) where - N: AbsDiffEq, + T: AbsDiffEq, ShapeConstraint: DimEq, { self.interval_fn( rhs, - |time| abs_diff_eq!(time, N::zero()), - |interval| abs_diff_eq!(interval, N::zero()), + |time| abs_diff_eq!(time, &T::zero()), + |interval| abs_diff_eq!(interval, &T::zero()), ) } - fn interval_fn(&self, rhs: &Self, is_present: P, is_lightlike: L) -> (N, LightCone) + fn interval_fn(&self, rhs: &Self, is_present: P, is_lightlike: L) -> (T, LightCone) where ShapeConstraint: DimEq, - P: Fn(N) -> bool, - L: Fn(N) -> bool, + P: Fn(&T) -> bool, + L: Fn(&T) -> bool, { - let time = self[0]; + let time = self[0].clone(); let difference = rhs - self; let interval = difference.scalar(&difference); - let light_cone = if is_lightlike(interval) { - if is_present(time) { + let light_cone = if is_lightlike(&interval) { + if is_present(&time) { LightCone::Origin } else if time.is_sign_positive() { LightCone::LightlikeFuture } else { LightCone::LightlikePast } - } else if interval.is_sign_positive() || is_present(time) { + } else if interval.is_sign_positive() || is_present(&time) { LightCone::Spacelike } else if time.is_sign_positive() { LightCone::TimelikeFuture @@ -636,32 +639,32 @@ where (interval, light_cone) } - fn new_boost(frame: &FrameN) -> Self + fn new_boost(frame: &OFrame) -> Self where D: DimNameSub, ShapeConstraint: AreMultipliable + DimEq, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { - let &FrameN { + let OFrame { zeta_cosh, zeta_sinh, - ref axis, + axis, } = frame; let mut b = Self::zeros(); for (i, u) in axis.iter().enumerate() { - b[(i + 1, 0)] = *u; - b[(0, i + 1)] = *u; + b[(i + 1, 0)] = u.clone(); + b[(0, i + 1)] = u.clone(); } let uk = b.clone_owned(); - b.gemm(zeta_cosh - N::one(), &uk, &uk, -zeta_sinh); + b.gemm(zeta_cosh.clone() - T::one(), &uk, &uk, -zeta_sinh.clone()); for i in 0..D::dim() { - b[(i, i)] += N::one(); + b[(i, i)] += T::one(); } b } #[inline] - fn boost(&self, frame: &FrameN) -> Self + fn boost(&self, frame: &OFrame) -> Self where R: DimNameSub, D: DimNameSub, @@ -669,14 +672,14 @@ where + SameNumberOfColumns + DimEq<>::Output, >::Output> + SameNumberOfRows<>::Output, >::Output>, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { let mut v = self.clone_owned(); v.boost_mut(frame); v } - fn boost_mut(&mut self, frame: &FrameN) + fn boost_mut(&mut self, frame: &OFrame) where R: DimNameSub, D: DimNameSub, @@ -684,87 +687,87 @@ where + SameNumberOfColumns + DimEq<>::Output, >::Output> + SameNumberOfRows<>::Output, >::Output>, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { - let &FrameN { + let OFrame { zeta_cosh, zeta_sinh, - ref axis, + axis, } = frame; let u = axis.as_ref(); - let a = self[0]; + let a = self[0].clone(); let (rows, _cols) = self.shape_generic(); let zu = self.rows_generic(1, rows.sub(U1)).dot(u); - self[0] = zeta_cosh * a - zeta_sinh * zu; + self[0] = zeta_cosh.clone() * a.clone() - zeta_sinh.clone() * zu.clone(); let mut z = self.rows_generic_mut(1, rows.sub(U1)); - z += u * ((zeta_cosh - N::one()) * zu - zeta_sinh * a); + z += u * ((zeta_cosh.clone() - T::one()) * zu - zeta_sinh.clone() * a); } #[inline] - fn new_velocity(frame: &FrameN) -> OVector + fn new_velocity(frame: &OFrame) -> OVector where D: DimNameSub, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, - DefaultAllocator: Allocator> + Allocator, + DefaultAllocator: Allocator> + Allocator, { frame.velocity() } #[inline] - fn frame(&self) -> FrameN + fn frame(&self) -> OFrame where R: DimNameSub, ShapeConstraint: SameNumberOfColumns, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { - FrameN::::from_velocity(self) + OFrame::::from_velocity(self) } #[inline] - fn from_split(temporal: &N, spatial: &OVector>) -> OVector + fn from_split(temporal: &T, spatial: &OVector>) -> OVector where D: DimNameSub, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, - DefaultAllocator: Allocator> + Allocator, - >::Buffer: - StorageMut, + DefaultAllocator: Allocator> + Allocator, + >::Buffer: + StorageMut, { - let mut v = OVector::::zeros(); - *v.temporal_mut() = *temporal; + let mut v = OVector::::zeros(); + *v.temporal_mut() = temporal.clone(); v.spatial_mut().copy_from(spatial); v } #[inline] - fn split(&self) -> (&N, MatrixView, C, U1, R>) + fn split(&self) -> (&T, MatrixView, C, U1, R>) where R: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator, - >::Buffer: - Storage, + DefaultAllocator: Allocator, + >::Buffer: + Storage, { (self.temporal(), self.spatial()) } #[inline] - fn split_mut(&mut self) -> (&mut N, MatrixViewMut, C>) + fn split_mut(&mut self) -> (&mut T, MatrixViewMut, C>) where R: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator, - >::Buffer: - Storage, + DefaultAllocator: Allocator, + >::Buffer: + Storage, { let (temporal, spatial) = self.as_mut_slice().split_at_mut(1); ( temporal.get_mut(0).unwrap(), - MatrixViewMut::, C>::from_slice(spatial), + MatrixViewMut::, C>::from_slice(spatial), ) } #[inline] - fn temporal(&self) -> &N + fn temporal(&self) -> &T where R: DimNameSub, ShapeConstraint: DimEq, @@ -773,7 +776,7 @@ where } #[inline] - fn temporal_mut(&mut self) -> &mut N + fn temporal_mut(&mut self) -> &mut T where R: DimNameSub, ShapeConstraint: DimEq, @@ -782,35 +785,30 @@ where } #[inline] - fn spatial(&self) -> MatrixView, C, U1, R> + fn spatial(&self) -> MatrixView, C, U1, R> where R: DimNameSub, ShapeConstraint: DimEq, - >::Buffer: - RawStorage, + >::Buffer: + RawStorage, { let (rows, _cols) = self.shape_generic(); self.rows_generic(1, rows.sub(U1)) } #[inline] - fn spatial_mut(&mut self) -> MatrixViewMut, C, U1, R> + fn spatial_mut(&mut self) -> MatrixViewMut, C, U1, R> where R: DimNameSub, ShapeConstraint: DimEq, - >::Buffer: - RawStorageMut, + >::Buffer: + RawStorageMut, { let (rows, _cols) = self.shape_generic(); self.rows_generic_mut(1, rows.sub(U1)) } } -#[inline] -fn neg(n: &mut N) { - *n = -n.clone(); -} - /// Spacetime regions regarding an event's light cone. #[derive(Debug, PartialEq, Eq, Clone, Copy)] pub enum LightCone { @@ -830,19 +828,19 @@ pub enum LightCone { /// Inertial frame of reference in $2$-dimensional Lorentzian space /// $\R^{-,+} = \R^{1,1}$. -pub type Frame2 = FrameN; +pub type Frame2 = OFrame; /// Inertial frame of reference in $3$-dimensional Lorentzian space /// $\R^{-,+} = \R^{1,2}$. -pub type Frame3 = FrameN; +pub type Frame3 = OFrame; /// Inertial frame of reference in $4$-dimensional Lorentzian space /// $\R^{-,+} = \R^{1,3}$. -pub type Frame4 = FrameN; +pub type Frame4 = OFrame; /// Inertial frame of reference in $5$-dimensional Lorentzian space /// $\R^{-,+} = \R^{1,4}$. -pub type Frame5 = FrameN; +pub type Frame5 = OFrame; /// Inertial frame of reference in $6$-dimensional Lorentzian space /// $\R^{-,+} = \R^{1,5}$. -pub type Frame6 = FrameN; +pub type Frame6 = OFrame; /// Inertial frame of reference in $n$-dimensional Lorentzian space $\R^{-,+} = \R^{1,n-1}$. /// @@ -860,39 +858,39 @@ pub type Frame6 = FrameN; /// /// Where $\gamma$ is the Lorentz factor. #[derive(Debug, PartialEq, Eq, Clone)] -pub struct FrameN +pub struct OFrame where - N: Scalar, + T: Scalar, D: DimNameSub, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { - zeta_cosh: N, - zeta_sinh: N, - axis: Unit>>, + zeta_cosh: T, + zeta_sinh: T, + axis: Unit>>, } -impl FrameN +impl OFrame where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { /// Inertial frame of reference with velocity $u^\mu$. #[must_use] - pub fn from_velocity(u: &OMatrix) -> Self + pub fn from_velocity(u: &OMatrix) -> Self where R: DimNameSub, C: Dim, ShapeConstraint: SameNumberOfRows + SameNumberOfColumns, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { let mut scaled_axis = OVector::zeros(); - let zeta_cosh = u[0]; + let zeta_cosh = u[0].clone(); let (rows, _cols) = u.shape_generic(); scaled_axis .iter_mut() .zip(u.rows_generic(1, rows.sub(U1)).iter()) - .for_each(|(scaled_axis, &u)| *scaled_axis = u); + .for_each(|(scaled_axis, u)| *scaled_axis = u.clone()); let (axis, zeta_sinh) = Unit::new_and_get(scaled_axis); Self { zeta_cosh, @@ -904,7 +902,7 @@ where /// Inertial frame of reference with rapidity $\vec \zeta$. #[must_use] #[inline] - pub fn from_zeta(scaled_axis: OVector>) -> Self { + pub fn from_zeta(scaled_axis: OVector>) -> Self { let (axis, zeta) = Unit::new_and_get(scaled_axis); Self::from_axis_zeta(axis, zeta) } @@ -912,9 +910,9 @@ where /// Inertial frame of reference along $\hat u$ with rapidity $\zeta$. #[must_use] #[inline] - pub fn from_axis_zeta(axis: Unit>>, zeta: N) -> Self { + pub fn from_axis_zeta(axis: Unit>>, zeta: T) -> Self { Self { - zeta_cosh: zeta.cosh(), + zeta_cosh: zeta.clone().cosh(), zeta_sinh: zeta.sinh(), axis, } @@ -923,7 +921,7 @@ where /// Inertial frame of reference with velocity ratio $\vec \beta$. #[must_use] #[inline] - pub fn from_beta(scaled_axis: OVector>) -> Self { + pub fn from_beta(scaled_axis: OVector>) -> Self { let (axis, beta) = Unit::new_and_get(scaled_axis); Self::from_axis_beta(axis, beta) } @@ -931,14 +929,14 @@ where /// Inertial frame of reference along $\hat u$ with velocity ratio $\beta$. #[must_use] #[inline] - pub fn from_axis_beta(axis: Unit>>, beta: N) -> Self { + pub fn from_axis_beta(axis: Unit>>, beta: T) -> Self { debug_assert!( - -N::one() < beta && beta < N::one(), + -T::one() < beta && beta < T::one(), "Velocity ratio `beta` is out of range (-1, +1)" ); - let gamma = N::one() / (N::one() - beta * beta).sqrt(); + let gamma = T::one() / (T::one() - beta.clone() * beta.clone()).sqrt(); Self { - zeta_cosh: gamma, + zeta_cosh: gamma.clone(), zeta_sinh: beta * gamma, axis, } @@ -946,53 +944,53 @@ where /// Velocity $u^\mu$. #[must_use] - pub fn velocity(&self) -> OVector + pub fn velocity(&self) -> OVector where - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { - let mut u = OVector::::zeros(); + let mut u = OVector::::zeros(); u[0] = self.gamma(); let (rows, _cols) = u.shape_generic(); u.rows_generic_mut(1, rows.sub(U1)) .iter_mut() .zip(self.axis.iter()) - .for_each(|(u, &axis)| *u = self.beta_gamma() * axis); + .for_each(|(u, axis)| *u = self.beta_gamma() * axis.clone()); u } /// Direction $\hat u$. #[must_use] #[inline] - pub fn axis(&self) -> Unit>> { + pub fn axis(&self) -> Unit>> { self.axis.clone() } /// Rapidity $\zeta$. #[must_use] #[inline] - pub fn zeta(&self) -> N { + pub fn zeta(&self) -> T { self.beta().atanh() } /// Velocity ratio $\beta$. #[must_use] #[inline] - pub fn beta(&self) -> N { + pub fn beta(&self) -> T { self.beta_gamma() / self.gamma() } /// Lorentz factor $\gamma$. #[must_use] #[inline] - pub const fn gamma(&self) -> N { - self.zeta_cosh + pub fn gamma(&self) -> T { + self.zeta_cosh.clone() } /// Product of velocity ratio $\beta$ and Lorentz factor $\gamma$. #[must_use] #[inline] - pub const fn beta_gamma(&self) -> N { - self.zeta_sinh + pub fn beta_gamma(&self) -> T { + self.zeta_sinh.clone() } /// Relativistic velocity addition `self`$\oplus$`frame`. @@ -1002,7 +1000,7 @@ where #[inline] pub fn compose(&self, frame: &Self) -> Self where - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { frame.velocity().boost(&-self.clone()).frame() } @@ -1017,14 +1015,14 @@ where /// /// See [`Self::rotation`] for further details. #[must_use] - pub fn angle(&self, frame: &Self) -> N + pub fn angle(&self, frame: &Self) -> T where - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { let (u, v) = (self, frame); - let ucv = u.compose(v); - let vcu = v.compose(u); - ucv.axis().angle(&vcu.axis()) + let ucv = u.compose(v).axis(); + let vcu = v.compose(u).axis(); + ucv.dot(&vcu).clamp(-T::one(), T::one()).acos() } /// $ @@ -1059,20 +1057,22 @@ where /// embeds the spatial rotation matrix /// /// $$ - /// \Kuv = I + \sin \epsilon (\hat v \hat u^T - \hat u \hat v^T) - /// + (\cos \epsilon - 1) (\hat u \hat u^T + \hat v \hat v^T) + /// \Kuv = I + (\hat b\hat a^T - \hat a\hat b^T) + \frac{1}{1+\hat a\cdot\hat b} \left[ + /// (\hat a\cdot\hat b)(\hat b\hat a^T+\hat a\hat b^T) - (\hat a\hat a^T + \hat b\hat b^T) + /// \right] /// $$ /// - /// rotating in the plane spanned by the orthonormalized pair $\hat u = \hat \beta_u$ and - /// $\hat v = \widehat{\hat \beta_v - (\hat \beta_u \cdot \hat \beta_v) \hat \beta_u}$. The - /// rotation angle $\epsilon$ is found according to [`Self::angle`]. + /// rotating in the plane spanned by the arc from + /// $\hat a = \widehat{\vec\beta_u\oplus\vec\beta_v}$ to + /// $\hat b = \widehat{\vec\beta_v\oplus\vec\beta_u}$ with the rotation angle + /// $\epsilon = \arccos(\hat a\cdot\hat b)$. /// /// See [`Self::axis_angle`] for $4$-dimensional specialization. /// /// ``` /// use approx::{assert_ulps_eq, assert_ulps_ne}; /// use nalgebra::{Matrix4, Vector3}; - /// use nalgebra_spacetime::{Frame4, LorentzianN}; + /// use nalgebra_spacetime::{Frame4, Lorentzian}; /// /// let u = Frame4::from_beta(Vector3::new(0.18, 0.73, 0.07)); /// let v = Frame4::from_beta(Vector3::new(0.41, 0.14, 0.25)); @@ -1093,27 +1093,27 @@ where /// assert_ulps_eq!(boost_vcu * matrix_ucv, boost_v * boost_u); /// ``` #[must_use] - pub fn rotation(&self, frame: &Self) -> (OMatrix, N) + #[allow(clippy::similar_names, clippy::many_single_char_names)] + pub fn rotation(&self, frame: &Self) -> (OMatrix, T) where - DefaultAllocator: Allocator - + Allocator - + Allocator> - + Allocator, DimNameDiff>, + DefaultAllocator: Allocator + + Allocator> + + Allocator + + Allocator, DimNameDiff>, { - let [u, v] = [self, frame]; - let ang = u.angle(v); - let (sin, cos) = ang.sin_cos(); - let u = u.axis().into_inner(); - let v = v.axis().into_inner(); - let v = (&v - &u * u.dot(&v)).normalize(); - let ut = u.transpose(); - let vt = v.transpose(); - let rot = (&v * &ut - &u * &vt) * sin + (&u * &ut + &v * &vt) * (cos - N::one()); - let mut mat = OMatrix::::identity(); + let (u, v) = (self, frame); + let a = &u.compose(v).axis().into_inner(); + let b = &v.compose(u).axis().into_inner(); + let ab = a.dot(b); + let d = T::one() + ab.clone(); + let [at, bt] = &[a.transpose(), b.transpose()]; + let [b_at, a_bt, a_at, b_bt] = &[b * at, a * bt, a * at, b * bt]; + let [k1, k2] = [b_at - a_bt, (b_at + a_bt) * ab.clone() - (a_at + b_bt)]; + let mut mat = OMatrix::::identity(); let (r, c) = mat.shape_generic(); - let mut sub = mat.generic_view_mut((1, 1), (r.sub(U1), c.sub(U1))); - sub += rot; - (mat, ang) + let mut rot = mat.generic_view_mut((1, 1), (r.sub(U1), c.sub(U1))); + rot += k1 + k2 / d; + (mat, ab.clamp(-T::one(), T::one()).acos()) } /// Wigner rotation axis $\widehat {\vec \beta_u \times \vec \beta_v}$ and angle $\epsilon$ of @@ -1136,7 +1136,7 @@ where /// ``` /// use approx::{assert_abs_diff_ne, assert_ulps_eq}; /// use nalgebra::Vector3; - /// use nalgebra_spacetime::{Frame4, LorentzianN}; + /// use nalgebra_spacetime::{Frame4, Lorentzian}; /// /// let u = Frame4::from_beta(Vector3::new(0.18, 0.73, 0.07)); /// let v = Frame4::from_beta(Vector3::new(0.41, 0.14, 0.25)); @@ -1152,12 +1152,12 @@ where /// assert_ulps_eq!(axis, ucv.cross(&vcu).normalize(), epsilon = 1e-15); /// ``` #[must_use] - pub fn axis_angle(&self, frame: &Self) -> (Unit>>, N) + pub fn axis_angle(&self, frame: &Self) -> (Unit>>, T) where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, ShapeConstraint: DimEq, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { let (u, v) = (self, frame); let (axis, sin) = Unit::new_and_get(u.axis().cross(&v.axis())); @@ -1165,44 +1165,44 @@ where let bg = u.beta_gamma() * v.beta_gamma(); let ug = u.gamma(); let vg = v.gamma(); - let cg = ug * vg + bg * uv; - let sum = N::one() + cg + ug + vg; - let prod = (N::one() + cg) * (N::one() + ug) * (N::one() + vg); + let cg = ug.clone() * vg.clone() + bg.clone() * uv; + let sum = T::one() + cg.clone() + ug.clone() + vg.clone(); + let prod = (T::one() + cg) * (T::one() + ug) * (T::one() + vg); (axis, (sum / prod * bg * sin).asin()) } } -impl From> for FrameN +impl From> for OFrame where - N: SimdRealField + Signed + Real, + T: RealField, R: DimNameSub, - C: DimName, + C: DimNameSub, ShapeConstraint: SameNumberOfColumns, - DefaultAllocator: Allocator + Allocator>, + DefaultAllocator: Allocator + Allocator>, { #[inline] - fn from(u: OMatrix) -> Self { + fn from(u: OMatrix) -> Self { u.frame() } } -impl From> for OVector +impl From> for OVector where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator> + Allocator, + DefaultAllocator: Allocator> + Allocator, { #[inline] - fn from(frame: FrameN) -> Self { + fn from(frame: OFrame) -> Self { frame.velocity() } } -impl Neg for FrameN +impl Neg for OFrame where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { type Output = Self; @@ -1213,11 +1213,11 @@ where } } -impl Add for FrameN +impl Add for OFrame where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator> + Allocator, + DefaultAllocator: Allocator> + Allocator, { type Output = Self; @@ -1227,11 +1227,11 @@ where } } -impl Sub for FrameN +impl Sub for OFrame where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator> + Allocator, + DefaultAllocator: Allocator> + Allocator, { type Output = Self; @@ -1241,12 +1241,12 @@ where } } -impl Copy for FrameN +impl Copy for OFrame where - N: SimdRealField + Signed + Real, + T: RealField + Copy, D: DimNameSub, - DefaultAllocator: Allocator>, - Owned>: Copy, + DefaultAllocator: Allocator>, + Owned>: Copy, { } @@ -1271,43 +1271,44 @@ where /// $$ /// /// With $n$-velocity $u^\mu$, Lorentz factor $\gamma$, and velocity ratio $\vec \beta$. +#[repr(transparent)] #[derive(Debug, PartialEq, Clone)] -pub struct MomentumN +pub struct OMomentum where - N: Scalar, + T: Scalar, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { - momentum: OVector, + momentum: OVector, } -impl MomentumN +impl OMomentum where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { - /// Momentum with spacetime [`LorentzianN::split`], `energy` $E$ and + /// Momentum with spacetime [`Lorentzian::split`], `energy` $E$ and /// `momentum` $\vec p$. #[must_use] #[inline] - pub fn from_split(energy: &N, momentum: &OVector>) -> Self + pub fn from_split(energy: &T, momentum: &OVector>) -> Self where - DefaultAllocator: Allocator>, - >::Buffer: - StorageMut, + DefaultAllocator: Allocator>, + >::Buffer: + StorageMut, { Self { - momentum: OVector::::from_split(energy, momentum), + momentum: OVector::::from_split(energy, momentum), } } /// Momentum $p^\mu=m u^\mu$ with rest `mass` $m$ at `velocity` $u^\mu$. #[must_use] #[inline] - pub fn from_mass_at_velocity(mass: N, velocity: OVector) -> Self + pub fn from_mass_at_velocity(mass: T, velocity: OVector) -> Self where - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { Self { momentum: velocity * mass, @@ -1319,9 +1320,9 @@ where /// Equals `frame.velocity() * mass`. #[must_use] #[inline] - pub fn from_mass_in_frame(mass: N, frame: &FrameN) -> Self + pub fn from_mass_in_frame(mass: T, frame: &OFrame) -> Self where - DefaultAllocator: Allocator>, + DefaultAllocator: Allocator>, { Self::from_mass_at_velocity(mass, frame.velocity()) } @@ -1329,8 +1330,8 @@ where /// Momentum $p^\mu$ with rest `mass` $m$ in center-of-momentum frame. #[must_use] #[inline] - pub fn from_mass_at_rest(mass: N) -> Self { - let mut momentum = OVector::::zeros(); + pub fn from_mass_at_rest(mass: T) -> Self { + let mut momentum = OVector::::zeros(); *momentum.temporal_mut() = mass; Self { momentum } } @@ -1338,66 +1339,66 @@ where /// Rest mass $m$ as timelike norm $\sqrt{-p_\mu p^\mu}$ in spacelike sign convention. #[must_use] #[inline] - pub fn mass(&self) -> N { + pub fn mass(&self) -> T { self.momentum.timelike_norm() } /// Velocity $u^\mu$ as momentum $p^\mu$ divided by rest `mass()` $m$. #[must_use] #[inline] - pub fn velocity(&self) -> OVector { + pub fn velocity(&self) -> OVector { self.momentum.clone() / self.mass() } - /// Energy $E$ as [`LorentzianN::temporal`] component. + /// Energy $E$ as [`Lorentzian::temporal`] component. #[must_use] #[inline] - pub fn energy(&self) -> &N { + pub fn energy(&self) -> &T { self.momentum.temporal() } - /// Momentum $\vec p$ as [`LorentzianN::spatial`] components. + /// Momentum $\vec p$ as [`Lorentzian::spatial`] components. #[must_use] #[inline] - pub fn momentum(&self) -> VectorView, U1, D> + pub fn momentum(&self) -> VectorView, U1, D> where - DefaultAllocator: Allocator, - >::Buffer: - Storage, + DefaultAllocator: Allocator, + >::Buffer: + Storage, { self.momentum.spatial() } } -impl From> for MomentumN +impl From> for OMomentum where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { #[inline] - fn from(momentum: OVector) -> Self { + fn from(momentum: OVector) -> Self { Self { momentum } } } -impl From> for OVector +impl From> for OVector where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { #[inline] - fn from(momentum: MomentumN) -> Self { + fn from(momentum: OMomentum) -> Self { momentum.momentum } } -impl Neg for MomentumN +impl Neg for OMomentum where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { type Output = Self; @@ -1408,11 +1409,11 @@ where } } -impl Add for MomentumN +impl Add for OMomentum where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { type Output = Self; @@ -1424,11 +1425,11 @@ where } } -impl Sub for MomentumN +impl Sub for OMomentum where - N: SimdRealField + Signed + Real, + T: RealField, D: DimNameSub, - DefaultAllocator: Allocator, + DefaultAllocator: Allocator, { type Output = Self; @@ -1440,22 +1441,22 @@ where } } -impl Copy for MomentumN +impl Copy for OMomentum where - N: SimdRealField + Signed + Real, + T: RealField + Copy, D: DimNameSub, - DefaultAllocator: Allocator, - Owned: Copy, + DefaultAllocator: Allocator, + Owned: Copy, { } /// Momentum in $2$-dimensional Lorentzian space $\R^{-,+} = \R^{1,1}$. -pub type Momentum2 = MomentumN; +pub type Momentum2 = OMomentum; /// Momentum in $3$-dimensional Lorentzian space $\R^{-,+} = \R^{1,2}$. -pub type Momentum3 = MomentumN; +pub type Momentum3 = OMomentum; /// Momentum in $4$-dimensional Lorentzian space $\R^{-,+} = \R^{1,3}$. -pub type Momentum4 = MomentumN; +pub type Momentum4 = OMomentum; /// Momentum in $5$-dimensional Lorentzian space $\R^{-,+} = \R^{1,4}$. -pub type Momentum5 = MomentumN; +pub type Momentum5 = OMomentum; /// Momentum in $6$-dimensional Lorentzian space $\R^{-,+} = \R^{1,5}$. -pub type Momentum6 = MomentumN; +pub type Momentum6 = OMomentum;