very basic GLON/GLAT support. No equinox is supported

src/coo_system.rs

@@ -1,4 +1,8 @@
+use std::f64::consts::PI;
+
 use fitsrs::hdu::header::{extension::image::Image, Header};
+use mapproj::LonLat;
+use mapproj::XYZ;
 
 use crate::error::Error;
 use crate::utils;
@@ -32,12 +36,16 @@ impl RadeSys {
 }
 
 pub enum CooSystem {
+    /// ICRS/J2000
     EQUATORIAL,
     GALACTIC,
     ECLIPTIC,
     HELIOECLIPTIC,
     SUPERGALACTIC,
-    CUSTOM { radesys: RadeSys, equinox: f64 },
+    CUSTOM {
+        radesys: RadeSys,
+        equinox: f64,
+    },
 }
 
 impl CooSystem {
@@ -62,4 +70,84 @@ impl CooSystem {
 
         Ok(coo_system)
     }
+
+    /// Convert a lonlat expressed in icrs coordinate system to the
+    /// self system
+    pub fn from_icrs(&self, lonlat: LonLat) -> LonLat {
+        match self {
+            CooSystem::EQUATORIAL => lonlat,
+            CooSystem::GALACTIC => {
+                let xyz = lonlat.to_xyz();
+                // ICRS_2_GAL * xyz
+                let rotated_xyz = XYZ::new(
+                    ICRS_2_GAL[0] * xyz.x() + ICRS_2_GAL[1] * xyz.y() + ICRS_2_GAL[2] * xyz.z(),
+                    ICRS_2_GAL[3] * xyz.x() + ICRS_2_GAL[4] * xyz.y() + ICRS_2_GAL[5] * xyz.z(),
+                    ICRS_2_GAL[6] * xyz.x() + ICRS_2_GAL[7] * xyz.y() + ICRS_2_GAL[8] * xyz.z(),
+                );
+
+                xyz_to_lonlat(&rotated_xyz)
+            }
+            _ => {
+                // todo
+                lonlat
+            }
+        }
+    }
+
+    /// Convert a lonlat expressed in the self system to the icrs coo system
+    pub fn to_icrs(&self, lonlat: LonLat) -> LonLat {
+        match self {
+            CooSystem::EQUATORIAL => lonlat,
+            CooSystem::GALACTIC => {
+                let xyz = lonlat.to_xyz();
+                // ICRS_2_GAL * xyz
+                let rotated_xyz = XYZ::new(
+                    GAL_2_ICRS[0] * xyz.x() + GAL_2_ICRS[1] * xyz.y() + GAL_2_ICRS[2] * xyz.z(),
+                    GAL_2_ICRS[3] * xyz.x() + GAL_2_ICRS[4] * xyz.y() + GAL_2_ICRS[5] * xyz.z(),
+                    GAL_2_ICRS[6] * xyz.x() + GAL_2_ICRS[7] * xyz.y() + GAL_2_ICRS[8] * xyz.z(),
+                );
+
+                xyz_to_lonlat(&rotated_xyz)
+            }
+            _ => {
+                // todo
+                lonlat
+            }
+        }
+    }
+}
+
+const GAL_2_ICRS: &[f64; 9] = &[
+    -0.4448296299195045,
+    -0.1980763734646737,
+    -0.873437090247923,
+    0.7469822444763707,
+    0.4559837762325372,
+    -0.4838350155267381,
+    0.4941094279435681,
+    -0.867_666_148_981_161,
+    -0.0548755604024359,
+];
+
+const ICRS_2_GAL: &[f64; 9] = &[
+    -0.4448296299195045,
+    0.7469822444763707,
+    0.4941094279435681,
+    -0.1980763734646737,
+    0.4559837762325372,
+    -0.867_666_148_981_161,
+    -0.873437090247923,
+    -0.4838350155267381,
+    -0.0548755604024359,
+];
+
+use mapproj::CustomFloat;
+fn xyz_to_lonlat(xyz: &XYZ) -> LonLat {
+    let r2 = xyz.x().pow2() + xyz.y().pow2();
+    // Latitude in [-pi/2, pi/2] (ok, since cos always positive here)
+    let lat = xyz.z().atan2(r2.sqrt());
+    // Compute the longitude in [-pi, pi]
+    let lon = xyz.y().atan2(xyz.x());
+    // Conforms to convention: Longitude in [0, 2*PI]
+    LonLat::new(if lon < 0.0 { 2.0 * PI + lon } else { lon }, lat)
 }

src/lib.rs

@@ -333,7 +333,6 @@ impl WCSProj {
 
         // 2. Identify the projection type
         let ctype1 = utils::retrieve_mandatory_parsed_keyword::<String>(header, "CTYPE1  ")?;
-        let _ = utils::retrieve_mandatory_parsed_keyword::<String>(header, "CTYPE2  ")?;
 
         let proj_name = &ctype1[5..=7];
 
@@ -420,12 +419,16 @@ impl WCSProj {
         Ok(WCSProj { proj, coo_system })
     }
 
-    /// Project a (lon, lat) 3D sphere position to get its corresponding location on the image
+    /// Project a (lon, lat) given in ICRS frame to get its corresponding location on the image
+    ///
     /// The result is given a (X, Y) tuple expressed in pixel coordinates.
     ///
-    /// # Param
-    /// * `lonlat`: the 3D sphere vertex expressed as a (lon, lat) tuple given in degrees
+    /// # Arguments
+    ///
+    /// * `lonlat`: a coo expressed as (lon, lat) tuple given in degrees and in ICRS system
     pub fn proj_lonlat(&self, lonlat: &LonLat) -> Option<ImgXY> {
+        let lonlat = &self.coo_system.from_icrs(lonlat.clone());
+
         let img_xy = match &self.proj {
             // Zenithal
             WCSCelestialProj::Azp(wcs) => wcs.lonlat2img(lonlat),
@@ -491,13 +494,15 @@ impl WCSProj {
     }
 
     /// Unproject a (X, Y) point from the image space to get its corresponding location on the sphere
-    /// The result is given a (lon, lat) tuple expressed in degrees.
     ///
-    /// # Param
+    /// The result is (lon, lat) tuple expressed in degrees in ICRS
+    ///
+    /// # Arguments
+    ///
     /// * `img_pos`: the image space point expressed as a (X, Y) tuple given en pixels
     pub fn unproj_lonlat(&self, img_pos: &ImgXY) -> Option<LonLat> {
         let img_pos = ImgXY::new(img_pos.x() + 1.0, img_pos.y() + 1.0);
-        match &self.proj {
+        let lonlat = match &self.proj {
             // Zenithal
             WCSCelestialProj::Azp(wcs) => wcs.img2lonlat(&img_pos),
             WCSCelestialProj::Szp(wcs) => wcs.img2lonlat(&img_pos),
@@ -556,7 +561,9 @@ impl WCSProj {
             WCSCelestialProj::CooSip(wcs) => wcs.img2lonlat(&img_pos),
             // Hybrid
             WCSCelestialProj::HpxSip(wcs) => wcs.img2lonlat(&img_pos),
-        }
+        };
+
+        lonlat.map(|ll| self.coo_system.to_icrs(ll))
     }
 
     /// Getter of the coordinate system