From d8170211000b3040d2a9ba09c18175aeef0559b8 Mon Sep 17 00:00:00 2001 From: Dar Dahlen Date: Fri, 25 Oct 2024 08:44:52 -0700 Subject: [PATCH] Add RA/DEC + rates to `SimultaneousStates` (#144) Add RA/DEC + rates to SimultaneousStates --- CHANGELOG.md | 8 ++ src/kete/rust/simult_states.rs | 28 +++++- src/kete/rust/vector.rs | 3 +- src/kete_core/src/frames/definitions.rs | 111 +++++++++++++++--------- src/kete_core/src/simult_states.rs | 45 ++++++++++ 5 files changed, 149 insertions(+), 46 deletions(-) diff --git a/CHANGELOG.md b/CHANGELOG.md index b13998e..1419bb2 100644 --- a/CHANGELOG.md +++ b/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] diff --git a/src/kete/rust/simult_states.rs b/src/kete/rust/simult_states.rs index 5d30238..59cdd87 100644 --- a/src/kete/rust/simult_states.rs +++ b/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> { + 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() { diff --git a/src/kete/rust/vector.rs b/src/kete/rust/vector.rs index 3cbd41c..38274b4 100644 --- a/src/kete/rust/vector.rs +++ b/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 } diff --git a/src/kete_core/src/frames/definitions.rs b/src/kete_core/src/frames/definitions.rs index ddc264c..8e2ba65 100644 --- a/src/kete_core/src/frames/definitions.rs +++ b/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. /// -/// -pub fn from_polar_spherical(theta: f64, phi: f64) -> UnitVector3 { - let (theta_sin, theta_cos) = theta.sin_cos(); - let (phi_sin, phi_cos) = phi.sin_cos(); - UnitVector3::::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 { - 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 { - from_polar_spherical(PI / 2.0 - dec, ra) -} - -/// Convert a unit vector to polar spherical coordinates. +/// In the Equatorial frame, lat = dec, lon = ra /// /// -pub fn to_polar_spherical(vec: UnitVector3) -> (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) { - 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) { - let (mut dec, mut ra) = to_polar_spherical(vec); - if dec > PI { - dec = TAU - dec; - ra += PI +/// +#[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); + } + } + } } diff --git a/src/kete_core/src/simult_states.rs b/src/kete_core/src/simult_states.rs index bc08dbf..f9736d7 100644 --- a/src/kete_core/src/simult_states.rs +++ b/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> { + 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()) + } }