Add RA/DEC + rates to SimultaneousStates (#144)

Add RA/DEC + rates to SimultaneousStates

CHANGELOG.md

@@ -12,6 +12,8 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
 
 - Added support for loading orbit information from JPL Horizons
 - Added more complete testing for light delay computations in the various FOV checks.
+- Added quality of life function to `SimultaneousState` for computing the RA/Dec along
+  with the projected on sky rates of motion.
 
 ### Changed
 
@@ -28,6 +30,12 @@ and this project adheres to [Semantic Versioning](https://semver.org/spec/v2.0.0
   it was returning the position/velocity at the observed position, but the time was
   the instantaneous time.
 
+### Removed
+
+- Removed several polar coordinate transformations in the rust backend which were all
+  equivalent variations of latitude/longitude conversions. Nothing in the Python was
+  impacted.
+
 
 ## [v1.0.2]
 

src/kete/rust/simult_states.rs

@@ -14,7 +14,7 @@ use crate::{fovs::AllowedFOV, frame::PyFrames, state::PyState};
 /// The main value in this is that also includes an optional Field of View.
 /// If the FOV is provided, it is implied that the states which are present
 /// in this file were objects seen by the FOV.
-/// 
+///
 /// In the case where the FOV is provided, it is expected that the states
 /// positions will include light delay, so an object which is ~1au away from
 /// the FOV observer will have a JD which is offset by about 8 minutes.
@@ -158,6 +158,32 @@ impl PySimultaneousStates {
         Ok(vecs)
     }
 
+    /// If a FOV is present, calculate the RA/Decs and their rates for all states in this object.
+    /// This will automatically convert all frames to Equatorial.
+    ///
+    /// 4 numbers are returned for each object, [RA, DEC, RA', DEC'], where rates are provided in
+    /// degrees/day.
+    ///
+    /// The returned RA' rate is scaled by cos(dec) so that it is equivalent to a
+    /// linear projection onto the observing plane.
+    ///
+    #[getter]
+    pub fn ra_dec_with_rates(&self) -> PyResult<Vec<[f64; 4]>> {
+        Ok(self
+            .0
+            .ra_dec_with_rates()?
+            .into_iter()
+            .map(|[ra, dec, dra, ddec]| {
+                [
+                    ra.to_degrees(),
+                    dec.to_degrees(),
+                    dra.to_degrees(),
+                    ddec.to_degrees(),
+                ]
+            })
+            .collect())
+    }
+
     fn __repr__(&self) -> String {
         let n_states = self.0.states.len();
         let fov_str = match self.fov() {

src/kete/rust/vector.rs

@@ -185,8 +185,7 @@ impl Vector {
         if r < 1e-8 {
             return 0.0;
         }
-        ((FRAC_PI_2 - ((data.z / r).clamp(-1.0, 1.0)).acos()).to_degrees() + 180.0)
-            .rem_euclid(360.0)
+        ((FRAC_PI_2 - (data.z / r).clamp(-1.0, 1.0).acos()).to_degrees() + 180.0).rem_euclid(360.0)
             - 180.0
     }
 

src/kete_core/src/frames/definitions.rs

@@ -5,7 +5,7 @@
 use lazy_static::lazy_static;
 use nalgebra::{Matrix3, Rotation3, UnitVector3, Vector3};
 use serde::{Deserialize, Serialize};
-use std::f64::consts::{PI, TAU};
+use std::f64::consts::{FRAC_PI_2, PI, TAU};
 use std::fmt::{self, Debug, Display};
 
 use crate::prelude::{Error, KeteResult};
@@ -292,56 +292,32 @@ pub fn inertial_to_noninertial(
     (new_pos, new_vel)
 }
 
-/// Create a unit vector from polar spherical theta and phi angles in radians.
+/// Create a unit vector from latitude and longitude.
 ///
-/// <https://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates>
-pub fn from_polar_spherical(theta: f64, phi: f64) -> UnitVector3<f64> {
-    let (theta_sin, theta_cos) = theta.sin_cos();
-    let (phi_sin, phi_cos) = phi.sin_cos();
-    UnitVector3::<f64>::new_unchecked([theta_sin * phi_cos, theta_sin * phi_sin, theta_cos].into())
-}
-
-/// Create a unit vector from latitude and longitude in units of radians.
-pub fn from_lat_lon(lat: f64, lon: f64) -> UnitVector3<f64> {
-    from_polar_spherical(PI / 2.0 - lat, lon)
-}
-
-/// Create a unit vector from ra and dec in units of radians.
-pub fn from_ra_dec(ra: f64, dec: f64) -> UnitVector3<f64> {
-    from_polar_spherical(PI / 2.0 - dec, ra)
-}
-
-/// Convert a unit vector to polar spherical coordinates.
+/// In the Equatorial frame, lat = dec, lon = ra
 ///
 /// <https://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates>
-pub fn to_polar_spherical(vec: UnitVector3<f64>) -> (f64, f64) {
-    let theta = vec.z.acos();
-    let phi = vec.y.atan2(vec.x) % TAU;
-    (theta, phi)
+#[inline(always)]
+pub fn from_lat_lon(lat: f64, lon: f64) -> [f64; 3] {
+    let (lat_sin, lat_cos) = lat.sin_cos();
+    let (lon_sin, lon_cos) = lon.sin_cos();
+    [lat_cos * lon_cos, lat_cos * lon_sin, lat_sin]
 }
 
-/// Convert a unit vector to latitude and longitude.
+/// Convert a vector to latitude and longitude in the current coordinate frame.
 ///
-/// Input vector needs to be in the [`Frame::Ecliptic`] frame.
-pub fn to_lat_lon(vec: UnitVector3<f64>) -> (f64, f64) {
-    let (mut lat, mut lon) = to_polar_spherical(vec);
-    if lat > PI {
-        lat = TAU - lat;
-        lon += PI
-    }
-    (PI / 2.0 - lat, lon)
-}
-
-/// Convert a unit vector to ra and dec.
+/// In the Equatorial frame, lat = dec, lon = ra
 ///
-/// Input vector needs to be in the [`Frame::Equatorial`] frame.
-pub fn to_ra_dec(vec: UnitVector3<f64>) -> (f64, f64) {
-    let (mut dec, mut ra) = to_polar_spherical(vec);
-    if dec > PI {
-        dec = TAU - dec;
-        ra += PI
+/// <https://en.wikipedia.org/wiki/Spherical_coordinate_system#Cartesian_coordinates>
+#[inline(always)]
+pub fn to_lat_lon(x: f64, y: f64, z: f64) -> (f64, f64) {
+    let r = Vector3::from([x, y, z]).norm();
+    if r < 1e-10 {
+        return (0.0, 0.0);
     }
-    (ra, PI / 2.0 - dec)
+    let lon = (3.0 * FRAC_PI_2 - (z / r).clamp(-1.0, 1.0).acos()).rem_euclid(TAU) - PI;
+    let lat = y.atan2(x).rem_euclid(TAU);
+    (lon, lat)
 }
 
 /// Rotate a collection of vectors around the specified rotation vector.
@@ -363,6 +339,8 @@ pub fn rotate_around(
 
 #[cfg(test)]
 mod tests {
+    use std::f64::consts::{FRAC_PI_2, TAU};
+
     use crate::frames::*;
 
     #[test]
@@ -405,4 +383,51 @@ mod tests {
         assert!((0.2 - vel_return[1]).abs() <= 10.0 * f64::EPSILON);
         assert!((0.3 - vel_return[2]).abs() <= 10.0 * f64::EPSILON);
     }
+
+    #[test]
+    fn test_lat_lon() {
+        // Several cardinal axis checks.
+        let [x, y, z] = from_lat_lon(0.0, 0.0);
+        assert!((x - 1.0).abs() < 10.0 * f64::EPSILON);
+        assert!(y.abs() < 10.0 * f64::EPSILON);
+        assert!(z.abs() < 10.0 * f64::EPSILON);
+
+        let [x, y, z] = from_lat_lon(0.0, FRAC_PI_2);
+        assert!(x.abs() < 10.0 * f64::EPSILON);
+        assert!((y - 1.0).abs() < 10.0 * f64::EPSILON);
+        assert!(z.abs() < 10.0 * f64::EPSILON);
+
+        let [x, y, z] = from_lat_lon(0.0, -FRAC_PI_2);
+        assert!(x.abs() < 10.0 * f64::EPSILON);
+        assert!((y + 1.0).abs() < 10.0 * f64::EPSILON);
+        assert!(z.abs() < 10.0 * f64::EPSILON);
+
+        let [x, y, z] = from_lat_lon(FRAC_PI_2, 0.0);
+        assert!(x.abs() < 10.0 * f64::EPSILON);
+        assert!(y.abs() < 10.0 * f64::EPSILON);
+        assert!((z - 1.0).abs() < 10.0 * f64::EPSILON);
+
+        let [x, y, z] = from_lat_lon(-FRAC_PI_2, 0.0);
+        assert!(x.abs() < 10.0 * f64::EPSILON);
+        assert!(y.abs() < 10.0 * f64::EPSILON);
+        assert!((z + 1.0).abs() < 10.0 * f64::EPSILON);
+
+        // sample conversion from around the sphere
+        for idx in -10..10 {
+            let lat = idx as f64 / 11.0 * FRAC_PI_2;
+            for idy in 0..11 {
+                let lon = idy as f64 / 11. * TAU;
+                let [x, y, z] = from_lat_lon(lat, lon);
+                let (new_lat, new_lon) = to_lat_lon(x, y, z);
+
+                let [x_new, y_new, z_new] = from_lat_lon(new_lat, new_lon);
+                assert!((new_lat - lat).abs() < 10.0 * f64::EPSILON);
+                assert!((new_lon - lon).abs() < 10.0 * f64::EPSILON);
+
+                assert!((x - x_new).abs() < 10.0 * f64::EPSILON);
+                assert!((y - y_new).abs() < 10.0 * f64::EPSILON);
+                assert!((z - z_new).abs() < 10.0 * f64::EPSILON);
+            }
+        }
+    }
 }

src/kete_core/src/simult_states.rs

@@ -2,8 +2,11 @@
 //! These primarily exist to allow for easy read/write to binary formats.
 
 use crate::fov::FOV;
+use crate::frames::to_lat_lon;
 use crate::io::FileIO;
 use crate::prelude::{Error, Frame, KeteResult, State};
+use nalgebra::Vector3;
+use rayon::iter::{IntoParallelRefIterator, ParallelIterator};
 use serde::{Deserialize, Serialize};
 
 /// Collection of [`State`] at the same time.
@@ -73,4 +76,46 @@ impl SimultaneousStates {
             fov,
         })
     }
+
+    /// Compute RA/Dec along with linearized rates.
+    ///
+    /// Returns a vector containing [ra, dec, ra' * cos(dec), dec'], all vectors
+    /// are automatically cast into the equatorial frame.
+    /// The returned RA rate is scaled by cos(dec) so that it is equivalent to a
+    /// linear projection onto the observing plane.
+    pub fn ra_dec_with_rates(&self) -> KeteResult<Vec<[f64; 4]>> {
+        if self.fov.is_none() {
+            return Err(Error::ValueError(
+                "Field of view must be specified for the ra/dec to be computed.".into(),
+            ));
+        }
+        let fov = self.fov.as_ref().unwrap();
+
+        let mut obs = fov.observer().clone();
+        obs.try_change_frame_mut(Frame::Equatorial)?;
+
+        let obs_pos = Vector3::from(obs.pos);
+        let obs_vel = Vector3::from(obs.vel);
+
+        Ok(self
+            .states
+            .par_iter()
+            .map(|state| {
+                let mut state = state.clone();
+                state.try_change_frame_mut(Frame::Equatorial).unwrap();
+                let pos_diff = Vector3::from(state.pos) - obs_pos;
+                let vel_diff = Vector3::from(state.vel) - obs_vel;
+
+                let d_ra = (pos_diff.x * vel_diff.y - pos_diff.y * vel_diff.x)
+                    / (pos_diff.x.powi(2) + pos_diff.y.powi(2));
+                let r2 = pos_diff.norm_squared();
+
+                let d_dec = (vel_diff.z - pos_diff.z * pos_diff.dot(&vel_diff) / r2)
+                    / (r2 - pos_diff.z.powi(2)).sqrt();
+                let (dec, ra) = to_lat_lon(pos_diff[0], pos_diff[1], pos_diff[2]);
+
+                [ra, dec, d_ra * dec.cos(), d_dec]
+            })
+            .collect())
+    }
 }