hoomd_geometry/shape/
eighteight.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 [`EightEight`]
5
6use serde::{Deserialize, Serialize};
7
8use crate::IsPointInside;
9use hoomd_manifold::{Hyperbolic, Minkowski};
10use hoomd_vector::Metric;
11use std::f64::consts::PI;
12
13/// A regular octagon in two-dimensional hyperbolic space.
14///
15/// [`EightEight`] implements a single regular octagon in the {8,8} tiling of
16/// two-dimensional hyperbolic space. The scaling of the octagon is set such
17/// that each of the angles is $` \frac{2\pi}{8} `$ so that eight equivalent
18/// octagons meet at each vertex.
19#[derive(Clone, Debug, PartialEq, Serialize, Deserialize)]
20pub struct EightEight {}
21
22impl IsPointInside<Hyperbolic<3>> for EightEight {
23    /// Checks if a given Hyperbolic point is inside [`EightEight`].
24    ///
25    /// # Example
26    /// ```
27    /// use hoomd_geometry::{IsPointInside, shape::EightEight};
28    /// use hoomd_manifold::Hyperbolic;
29    /// use std::f64::consts::PI;
30    ///
31    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
32    /// let eight_eight = EightEight {};
33    ///
34    /// let point = Hyperbolic::<3>::from_polar_coordinates(1.0, PI / 8.0);
35    /// assert!(eight_eight.is_point_inside(&point));
36    /// # Ok(())
37    /// # }
38    /// ```
39    #[inline]
40    fn is_point_inside(&self, point: &Hyperbolic<3>) -> bool {
41        EightEight::distance_to_boundary(point) >= 0.0
42    }
43}
44
45impl EightEight {
46    /// Computes the shortest distance between a given point and the boundary
47    /// of `EightEight`.
48    ///
49    /// The shortest distance is computed by finding the arclength of the geodesic
50    /// which passes through the query point and intersects the boundary at a
51    /// right angle.
52    ///
53    /// # Example
54    /// ```
55    /// use approxim::assert_relative_eq;
56    /// use hoomd_geometry::shape::EightEight;
57    /// use hoomd_manifold::{Hyperbolic, Minkowski};
58    /// use std::f64::consts::PI;
59    ///
60    /// # fn main() -> Result<(), Box<dyn std::error::Error>> {
61    /// let v: f64 = EightEight::CUSP_TO_EDGE - 0.4;
62    /// let theta: f64 = PI / 8.0;
63    /// let x = Hyperbolic::from_minkowski_coordinates(
64    ///     [
65    ///         (v.sinh()) * (theta.cos()),
66    ///         (v.sinh()) * (theta.sin()),
67    ///         (v.cosh()),
68    ///     ]
69    ///     .into(),
70    /// );
71    /// assert_relative_eq!(
72    ///     EightEight::distance_to_boundary(&x),
73    ///     0.4,
74    ///     epsilon = 1e-12
75    /// );
76    /// # Ok(())
77    /// # }
78    /// ```
79    #[inline]
80    #[must_use]
81    pub fn distance_to_boundary(point: &Hyperbolic<3>) -> f64 {
82        let theta =
83            (point.coordinates()[1].atan2(point.coordinates()[0])).rem_euclid(PI / 4.0) - PI / 8.0;
84        let boost = (point.coordinates()[2]).acosh();
85        let (b_sinh, b_cosh) = (boost.sinh(), boost.cosh());
86        let xi = Self::CUSP_TO_EDGE;
87        let (xi_sinh, xi_cosh) = (xi.sinh(), xi.cosh());
88        // boost into frame where edge is the vertical diameter
89        let edge_as_diameter: Hyperbolic<3> =
90            Hyperbolic::<3>::from_minkowski_coordinates(Minkowski::from([
91                xi_cosh * b_sinh * (theta.cos()) - xi_sinh * b_cosh,
92                b_sinh * (theta.sin()),
93                -xi_sinh * b_sinh * (theta.cos()) + xi_cosh * b_cosh,
94            ]));
95        let flipped = Hyperbolic::<3>::from_minkowski_coordinates(Minkowski::from([
96            -edge_as_diameter.coordinates()[0],
97            edge_as_diameter.coordinates()[1],
98            edge_as_diameter.coordinates()[2],
99        ]));
100        let sign = -(edge_as_diameter.coordinates()[0]).signum();
101        sign * (edge_as_diameter.distance(&flipped)) / 2.0
102    }
103    /// Points on the boundary of the fundamental domain
104    #[inline]
105    #[must_use]
106    pub fn boundary_points(number_of_points: usize) -> Vec<(f64, f64)> {
107        let mut coords = Vec::<(f64, f64)>::new();
108        for n in 0..number_of_points {
109            let angle = (n as f64) * 2.0 * PI / (number_of_points as f64);
110            let tile_size = EightEight::EIGHTEIGHT;
111            let eta =
112                (tile_size.tanh() / (angle.cos() - angle.sin() * (1.0 - (2.0_f64).sqrt()))).atanh();
113            let x = eta.sinh() / (1.0 + eta.cosh());
114            for k in 0..8 {
115                coords.push((
116                    x * (angle + f64::from(k) * PI / 4.0).cos(),
117                    x * (angle + f64::from(k) * PI / 4.0).sin(),
118                ));
119            }
120        }
121        coords
122    }
123    /// Apply a lattice transformation to a point.
124    #[inline]
125    #[must_use]
126    pub fn gamma(eta: f64, theta: f64, point: &[f64; 3]) -> [f64; 3] {
127        let (eta_sinh_squared, two_eta_sinh, theta_sin, theta_cos) = (
128            (eta.sinh()).powi(2),
129            (2.0 * eta).sinh(),
130            theta.sin(),
131            theta.cos(),
132        );
133        [
134            (2.0 * (eta_sinh_squared) * ((theta_cos).powi(2)) + 1.0) * point[0]
135                + ((2.0 * theta).sin()) * (eta_sinh_squared) * point[1]
136                + (two_eta_sinh) * (theta_cos) * point[2],
137            ((2.0 * theta).sin()) * (eta_sinh_squared) * point[0]
138                + (2.0 * (eta_sinh_squared) * ((theta_sin).powi(2)) + 1.0) * point[1]
139                + (two_eta_sinh) * (theta_sin) * point[2],
140            (two_eta_sinh) * (theta_cos) * point[0]
141                + (two_eta_sinh) * (theta_sin) * point[1]
142                + ((2.0 * eta).cosh()) * point[2],
143        ]
144    }
145    /// Calculate the change in angle in the tangent bundle associated with a lattice
146    /// transformation.
147    #[inline]
148    #[must_use]
149    pub fn reorient(theta: f64, point: &[f64; 3]) -> f64 {
150        let (q_u, q_v) = (point[0] / (1.0 + point[2]), point[1] / (1.0 + point[2]));
151        let alpha = (1.0 + (PI / 4.0).cos()).sqrt();
152        let beta = (theta.cos()) * ((2.0 * ((PI / 4.0).cos())).sqrt());
153        let gamma = (theta.sin()) * ((2.0 * ((PI / 4.0).cos())).sqrt());
154        // let prefactor = 1.0/((b*v-c*u).powi(2)+(a+b*u+c*v).powi(2)).powi(2);
155        let p_x = alpha + beta * q_u + gamma * q_v;
156        let p_y = beta * q_v - gamma * q_u;
157        -2.0 * (p_y.atan2(p_x))
158    }
159    /// Cusp-to-vertex distance for the {8,8} tiling for Gauss curvature K = -1.
160    pub const EIGHTEIGHT: f64 = 2.448_452_447_678_076;
161    /// Cusp-to-middle-of-edge distance for the {8,8} tiling for Gauss curvature K = -1.
162    pub const CUSP_TO_EDGE: f64 = 1.528_570_919_480_998;
163    /// Length of one of the sides of the {8,8} tiling for Gauss curvature K = -1.
164    pub const EDGE_LENGTH: f64 = 3.057_141_838_961_997;
165}
166
167#[cfg(test)]
168mod tests {
169    use super::*;
170    use approxim::assert_relative_eq;
171    use hoomd_manifold::{Hyperbolic, HyperbolicDisk, Minkowski};
172    use rand::{SeedableRng, distr::Distribution, rngs::StdRng};
173    use std::ops::Not;
174
175    #[test]
176    fn boundary_distance() {
177        // Distance to the edge of the {8,8} fundamental domain
178        let e = Hyperbolic::<3>::from_polar_coordinates(1.0, 0.1);
179        let e_edge_distance = EightEight::distance_to_boundary(&e);
180        let e_edge_distance_numeric = 0.631_401_734_734_821;
181        assert_relative_eq!(e_edge_distance, e_edge_distance_numeric, epsilon = 1e-12);
182
183        let f = Hyperbolic::<3>::from_polar_coordinates(0.6, 0.2 + PI / 4.0);
184        let f_edge_distance = EightEight::distance_to_boundary(&f);
185        let f_edge_distance_numeric = 0.947_879_122_461_848;
186        assert_relative_eq!(f_edge_distance, f_edge_distance_numeric, epsilon = 1e-12);
187    }
188
189    #[test]
190    fn inside_is_inside() {
191        let eight_eight = EightEight {};
192        let r = 1.528_570_919_480_998;
193        let mut rng = StdRng::seed_from_u64(239);
194        let disk = HyperbolicDisk {
195            disk_radius: r.try_into().expect("hard-coded positive number"),
196            point: Hyperbolic::<3>::from_minkowski_coordinates(Minkowski::from([0.0, 0.0, 1.0])),
197        };
198        let random_point: Hyperbolic<3> = disk.sample(&mut rng);
199        assert!(eight_eight.is_point_inside(&random_point));
200
201        let point_1 = Hyperbolic::<3>::from_polar_coordinates(1.52, PI / 8.0);
202        assert!(eight_eight.is_point_inside(&point_1));
203
204        let point_2 = Hyperbolic::<3>::from_polar_coordinates(2.44, PI / 4.0);
205        assert!(eight_eight.is_point_inside(&point_2));
206    }
207
208    #[test]
209    fn outside_is_outside() {
210        let eight_eight = EightEight {};
211        let point_1 = Hyperbolic::<3>::from_polar_coordinates(1.54, PI / 8.0);
212        assert!((eight_eight.is_point_inside(&point_1)).not());
213
214        let point_2 = Hyperbolic::<3>::from_polar_coordinates(2.45, PI / 4.0);
215        assert!((eight_eight.is_point_inside(&point_2)).not());
216    }
217}
hoomd_geometry/shape/eighteight.rs

hoomd_geometry/shape/
eighteight.rs
