Add 'intersect_great_circle' in polygon

CHANGELOG.md

@@ -1,5 +1,15 @@
 # `cdshealpix` Change Log
 
+## 0.6.9
+
+Released 2024-05-28
+
+### Add
+
+* Add `intersect_great_circle` as needed in Aladin Lite
+
+
+
 ## 0.6.8
 
 Released 2024-05-14

Cargo.toml

@@ -1,6 +1,6 @@
 [package]
 name = "cdshealpix"
-version = "0.6.8"
+version = "0.6.9"
 authors = ["F.-X. Pineau <[email protected]>"]
 edition = "2018"
 license = "Apache-2.0 OR MIT"

README.md

@@ -24,6 +24,7 @@ It is used in:
      + [MOCli](https://github.com/cds-astro/cds-moc-rust/tree/main/crates/cli) a standalone command line tool for linux, MacOS and Windows;
      + [MOCWasm](https://github.com/cds-astro/cds-moc-rust/tree/main/crates/wasm), a WASM library to manipulate MOCs from web browsers.
  * CDS internal developments
+ * *Please help me fill in this list*
 
 
 Initially, it is a port of a part of the CDS Java library available [here](https://github.com/cds-astro/cds-healpix-java),

src/nested/mod.rs

@@ -3812,6 +3812,7 @@ impl Layer {
         .into_boxed_slice(),
       south_pole_method,
     );
+
     let bounding_cone: Cone = Cone::bounding_cone(poly.vertices());
     let mut depth_start = 0;
     let neigs: Vec<u64> = if let ContainsSouthPoleMethod::Default = south_pole_method {

src/sph_geom/mod.rs

@@ -217,6 +217,34 @@ impl Polygon {
     vertices
   }
 
+  /// Returns all the intersections between this polygon edges and the great circle of given normal `n`.
+  pub fn intersect_great_circle(&self, n: &UnitVect3) -> Vec<UnitVect3> {
+    // Ensure a < b in longitude
+    let mut vertices = vec![];
+
+    let mut left = self.vertices.last().unwrap();
+    for (right, cross_prod) in self.vertices.iter().zip(self.cross_products.iter()) {
+      // Ensures pA < pB in longitude
+      let pa = left;
+      let pb = right;
+      let papb = dot_product(pa, pb);
+
+      let i = cross_product(cross_prod, n).normalized();
+
+      // Check if the intersection i belongs to the edge,
+      // else check if the intersection of the opposite of i belongs to the edge,
+      if dot_product(pa, &i) > papb && dot_product(pb, &i) > papb {
+        vertices.push(i);
+      } else if dot_product(pa, &i.opposite()) > papb && dot_product(pb, &i.opposite()) > papb {
+        vertices.push(i.opposite());
+      }
+
+      left = right
+    }
+
+    vertices
+  }
+
   /// Returns `true` if an edge of the polygon intersects the great-circle arc defined by the
   /// two given points (we consider the arc having a length < PI).
   pub fn is_intersecting_great_circle_arc(&self, a: &Coo3D, b: &Coo3D) -> bool {
@@ -243,7 +271,6 @@ impl Polygon {
     None
   }
 
-
   /* PREVIOUS CODE, KEPT THE TIME TO COMPARE WITH special_points_finder::intersect_parallel
   /// Returns the coordinate of the intersection from an edge of the polygon
   /// with a parallel defined by a latitude (this may suffer from numerical precision at poles