diffgeom/
coordinates.rs

1//! Module containing basic types representing coordinate systems.
2
3use super::tensors::{ContravariantIndex, CovariantIndex, Matrix, Tensor};
4use crate::typenum::consts::U2;
5use crate::typenum::uint::Unsigned;
6use crate::typenum::Pow;
7use generic_array::{ArrayLength, GenericArray};
8use std::fmt;
9use std::ops::{Index, IndexMut};
10
11/// `CoordinateSystem` marks a struct (usually a unit struct) as representing a coordinate system.
12pub trait CoordinateSystem: Sized {
13    /// An associated type representing the dimension of the coordinate system
14    type Dimension: Unsigned + ArrayLength<f64> + ArrayLength<usize>;
15
16    /// Function returning a small value for purposes of numerical differentiation.
17    /// What is considered a small value may depend on the point, hence the parameter.
18    /// Returns just 0.01 by default.
19    fn small(_: &Point<Self>) -> f64 {
20        0.01
21    }
22
23    /// Function returning the dimension
24    fn dimension() -> usize {
25        Self::Dimension::to_usize()
26    }
27}
28
29/// Struct representing a point on the manifold. The information about the coordinate system
30/// is saved in the type parameter, so that only operations on objects belonging to the same
31/// coordinate system will be allowed.
32pub struct Point<T: CoordinateSystem> {
33    /// The coordinates of the point.
34    x: GenericArray<f64, T::Dimension>,
35}
36
37impl<T> Point<T>
38where
39    T: CoordinateSystem,
40{
41    /// Creates a new point with coordinates described by the array
42    pub fn new(coords: GenericArray<f64, T::Dimension>) -> Point<T> {
43        Point { x: coords }
44    }
45
46    /// Creates a new point with coordinates passed in the slice
47    pub fn from_slice(coords: &[f64]) -> Point<T> {
48        Point {
49            x: GenericArray::clone_from_slice(coords),
50        }
51    }
52
53    /// Returns the point's coordinates as an array
54    pub fn coords_array(&self) -> &GenericArray<f64, T::Dimension> {
55        &self.x
56    }
57}
58
59impl<T> Clone for Point<T>
60where
61    T: CoordinateSystem,
62{
63    fn clone(&self) -> Point<T> {
64        Point::new(self.x.clone())
65    }
66}
67
68impl<T> Copy for Point<T>
69where
70    T: CoordinateSystem,
71    <T::Dimension as ArrayLength<f64>>::ArrayType: Copy,
72{
73}
74
75impl<T> Index<usize> for Point<T>
76where
77    T: CoordinateSystem,
78{
79    type Output = f64;
80
81    fn index(&self, idx: usize) -> &f64 {
82        &self.x[idx]
83    }
84}
85
86impl<T> IndexMut<usize> for Point<T>
87where
88    T: CoordinateSystem,
89{
90    fn index_mut(&mut self, idx: usize) -> &mut f64 {
91        &mut self.x[idx]
92    }
93}
94
95impl<T> PartialEq<Point<T>> for Point<T>
96where
97    T: CoordinateSystem,
98{
99    fn eq(&self, rhs: &Point<T>) -> bool {
100        (0..T::dimension()).all(|i| self[i] == rhs[i])
101    }
102}
103
104impl<T> Eq for Point<T> where T: CoordinateSystem {}
105
106impl<T> fmt::Debug for Point<T>
107where
108    T: CoordinateSystem,
109{
110    fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
111        write!(f, "Point{:?}", &self.x)
112    }
113}
114
115/// Trait used for conversions between different coordinate systems. Implementing `ConversionTo<T>`
116/// for a `CoordinateSystem` will allow objects in that system to be converted to the system `T`
117/// (note that `T` also has to be a `CoordinateSystem`).
118pub trait ConversionTo<T: CoordinateSystem + 'static>: CoordinateSystem
119where
120    T::Dimension: Pow<U2>,
121    <T::Dimension as Pow<U2>>::Output: ArrayLength<f64>,
122{
123    /// Function converting the coordinates of a point.
124    fn convert_point(p: &Point<Self>) -> Point<T>;
125
126    /// Function calculating a Jacobian at a point - that is, the matrix of derivatives
127    /// of the coordinate conversions.
128    ///
129    /// This will be contracted with contravariant indices in the tensor.
130    fn jacobian(p: &Point<Self>) -> Matrix<T> {
131        let d = Self::dimension();
132        let mut result = Matrix::zero(Self::convert_point(p));
133        let h = Self::small(p);
134
135        for j in 0..d {
136            let mut x = p.clone();
137            x[j] = x[j] - h;
138            let y1 = Self::convert_point(&x);
139
140            x[j] = x[j] + h * 2.0;
141            let y2 = Self::convert_point(&x);
142
143            for i in 0..d {
144                // calculate dyi/dxj
145                let index = [i, j];
146                result[&index[..]] = (y2[i] - y1[i]) / (2.0 * h);
147            }
148        }
149
150        result
151    }
152
153    /// The inverse matrix of the Jacobian at a point.
154    ///
155    /// In conversions, it will be contracted with covariant indices.
156    fn inv_jacobian(p: &Point<Self>) -> Tensor<T, (CovariantIndex, ContravariantIndex)> {
157        ConversionTo::<T>::jacobian(p).inverse().unwrap()
158    }
159}
diffgeom/coordinates.rs

diffgeom/
coordinates.rs
