WIP

crates/lox-earth/src/cio.rs

@@ -8,4 +8,4 @@
 
 //! Module cio exposes functions for calculating the Celestial Intermediate Origin (CIO) locator, s.
 
-mod s06;
+pub mod s06;

crates/lox-earth/src/lib.rs

@@ -11,5 +11,6 @@ pub mod cip;
 pub mod coordinate_transformations;
 pub mod nutation;
 pub mod rotation_angle;
+#[allow(dead_code)]
 pub mod tides;
 pub mod tio;

crates/lox-orbits/src/analysis.rs

@@ -7,7 +7,7 @@ use std::f64::consts::PI;
  * License, v. 2.0. If a copy of the MPL was not distributed with this
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
-use lox_bodies::{Origin, RotationalElements, Spheroid};
+use lox_bodies::{MaybeSpheroid, Origin, RotationalElements, Spheroid};
 use lox_math::roots::Brent;
 use lox_math::series::{Series, SeriesError};
 use lox_math::types::units::Radians;
@@ -19,9 +19,12 @@ use lox_time::TimeLike;
 use thiserror::Error;
 
 use crate::events::{find_windows, Window};
-use crate::frames::{BodyFixed, FrameTransformationProvider, Icrf, TryToFrame};
-use crate::ground::GroundLocation;
-use crate::trajectories::Trajectory;
+use crate::frames::{
+    BodyFixed, DynFrame, FrameTransformationProvider, Icrf, ReferenceFrame, TryRotateTo, TryToFrame,
+};
+use crate::ground::{DynGroundLocation, GroundLocation};
+use crate::states::State;
+use crate::trajectories::{DynTrajectory, Trajectory};
 
 #[derive(Debug, Clone, Error, PartialEq)]
 pub enum ElevationMaskError {
@@ -61,9 +64,59 @@ impl ElevationMask {
     }
 }
 
+pub fn elevation_dyn<T: TimeLike + TryToScale<Tdb, P> + Clone, P: FrameTransformationProvider>(
+    time: T,
+    gs: &DynGroundLocation,
+    mask: &ElevationMask,
+    sc: &DynTrajectory<T>,
+    provider: &P,
+) -> Radians {
+    let body_fixed = DynFrame::BodyFixed(gs.origin());
+    let sc = sc.interpolate_at(time.clone());
+    let rot = DynFrame::Icrf.try_rotation(&body_fixed, time, provider);
+    let (r1, v1) = rot.unwrap().rotate_state(sc.position(), sc.velocity());
+    let sc = State::new(sc.time(), r1, v1, sc.origin(), body_fixed);
+    let obs = gs.observables_dyn(sc);
+    obs.elevation() - mask.min_elevation(obs.azimuth())
+}
+
+pub fn visibility_dyn<T: TimeLike + TryToScale<Tdb, P> + Clone, P: FrameTransformationProvider>(
+    times: &[T],
+    gs: &DynGroundLocation,
+    mask: &ElevationMask,
+    sc: &DynTrajectory<T>,
+    provider: &P,
+) -> Vec<Window<T>> {
+    if times.len() < 2 {
+        return vec![];
+    }
+    let start = times.first().unwrap().clone();
+    let end = times.last().unwrap().clone();
+    let times: Vec<f64> = times
+        .iter()
+        .map(|t| (t.clone() - start.clone()).to_decimal_seconds())
+        .collect();
+    let root_finder = Brent::default();
+    find_windows(
+        |t| {
+            elevation_dyn(
+                start.clone() + TimeDelta::from_decimal_seconds(t).unwrap(),
+                gs,
+                mask,
+                sc,
+                provider,
+            )
+        },
+        start.clone(),
+        end.clone(),
+        &times,
+        root_finder,
+    )
+}
+
 pub fn elevation<
     T: TimeLike + TryToScale<Tdb, P> + Clone,
-    O: Origin + Spheroid + RotationalElements + Clone,
+    O: Origin + MaybeSpheroid + RotationalElements + Clone,
     P: FrameTransformationProvider,
 >(
     time: T,
@@ -136,14 +189,23 @@ mod tests {
     fn test_elevation() {
         let gs = ground_station_trajectory();
         let sc = spacecraft_trajectory();
+        let mask = ElevationMask::with_fixed_elevation(0.0);
         let expected: Vec<Radians> = include_str!("../../../data/elevation.csv")
             .lines()
             .map(|line| line.parse::<f64>().unwrap().to_radians())
             .collect();
         let actual: Vec<Radians> = gs
             .times()
             .iter()
-            .map(|t| elevation(*t, &gs, &sc, &NoOpFrameTransformationProvider))
+            .map(|t| {
+                elevation(
+                    *t,
+                    &location(),
+                    &mask,
+                    &sc,
+                    &NoOpFrameTransformationProvider,
+                )
+            })
             .collect();
         for (actual, expected) in actual.iter().zip(expected.iter()) {
             assert_close!(actual, expected, 1e-1);

crates/lox-orbits/src/frames.rs

@@ -8,14 +8,17 @@
 
 use std::{convert::Infallible, str::FromStr};
 
+use crate::frames::iau::icrf_to_bodyfixed;
+use crate::frames::iers::{cirf_to_tirf, icrf_to_cirf, tirf_to_itrf};
+use crate::rotations::Rotation;
 use glam::{DMat3, DVec3};
 use lox_bodies::{DynOrigin, MaybeRotationalElements, Origin, RotationalElements, Spheroid};
 use lox_math::types::units::Seconds;
-use lox_time::transformations::OffsetProvider;
+use lox_time::prelude::Tdb;
+use lox_time::transformations::{OffsetProvider, TryToScale};
+use lox_time::TimeLike;
 use thiserror::Error;
 
-use crate::rotations::Rotation;
-
 mod iau;
 mod iers;
 
@@ -205,10 +208,12 @@ impl ReferenceFrame for DynFrame {
 
 fn parse_iau_frame(s: &str) -> Option<DynFrame> {
     let (prefix, origin) = s.split_once("_")?;
+    dbg!(prefix);
+    dbg!(origin);
     if prefix.to_lowercase() != "iau" {
         return None;
     }
-    let origin: DynOrigin = origin.parse().ok()?;
+    let origin: DynOrigin = origin.to_lowercase().parse().ok()?;
     let _ = origin.maybe_rotational_elements(0.0)?;
     Some(DynFrame::BodyFixed(origin))
 }
@@ -237,33 +242,106 @@ impl FromStr for DynFrame {
     }
 }
 
-pub trait RotateTo<T: ReferenceFrame> {
-    type Time;
-
-    fn rotation(&self, frame: &T, time: Self::Time) -> Rotation;
-}
-
-pub trait TryRotateTo<T: ReferenceFrame> {
-    type Time;
+pub trait TryRotateTo<R: ReferenceFrame, P: FrameTransformationProvider> {
     type Error;
 
-    fn try_rotation(&self, frame: &T) -> Result<Rotation, Self::Error>;
+    fn try_rotation<T: TimeLike + Clone>(
+        &self,
+        frame: &R,
+        time: T,
+        provider: &P,
+    ) -> Result<Rotation, Self::Error>;
 }
 
-impl<T: ReferenceFrame, U: RotateTo<T>> TryRotateTo<T> for U {
-    type Error = Infallible;
-
-    fn try_rotation(&self, frame: &T) -> Result<Rotation, Self::Error> {
-        Ok(self.rotation(frame))
-    }
-}
-
-impl TryRotateTo<DynFrame> for DynFrame {
-    type Error = ();
-
-    fn try_rotation(&self, frame: &DynFrame) -> Result<Rotation, Self::Error> {
+impl<P: FrameTransformationProvider> TryRotateTo<DynFrame, P> for DynFrame {
+    // FIXME
+    type Error = String;
+
+    fn try_rotation<T: TimeLike + Clone>(
+        &self,
+        frame: &DynFrame,
+        time: T,
+        provider: &P,
+    ) -> Result<Rotation, Self::Error> {
+        // FIXME
+        let seconds_j2000 = time.seconds_since_j2000();
+        let centuries_j2000 = time.centuries_since_j2000();
         match self {
-            _ => todo!(),
+            DynFrame::Icrf => match frame {
+                DynFrame::Icrf => Ok(Rotation::IDENTITY),
+                DynFrame::Cirf => Ok(icrf_to_cirf(centuries_j2000)),
+                DynFrame::Tirf => {
+                    Ok(icrf_to_cirf(centuries_j2000).compose(&cirf_to_tirf(seconds_j2000)))
+                }
+                DynFrame::Itrf => Ok(icrf_to_cirf(centuries_j2000)
+                    .compose(&cirf_to_tirf(seconds_j2000))
+                    .compose(&tirf_to_itrf(centuries_j2000))),
+                DynFrame::BodyFixed(target) => icrf_to_bodyfixed(target, seconds_j2000).ok_or(
+                    format!("no rotational parameters available for {}", target.name()),
+                ),
+            },
+            DynFrame::Cirf => match frame {
+                DynFrame::Icrf => Ok(icrf_to_cirf(centuries_j2000).transpose()),
+                DynFrame::Cirf => Ok(Rotation::IDENTITY),
+                DynFrame::Tirf => Ok(cirf_to_tirf(seconds_j2000)),
+                DynFrame::Itrf => {
+                    Ok(cirf_to_tirf(seconds_j2000).compose(&tirf_to_itrf(centuries_j2000)))
+                }
+                DynFrame::BodyFixed(_) => Ok(self
+                    .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                    .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?)),
+            },
+            DynFrame::Tirf => match frame {
+                DynFrame::Icrf => Ok(cirf_to_tirf(seconds_j2000)
+                    .transpose()
+                    .compose(&icrf_to_cirf(centuries_j2000).transpose())),
+                DynFrame::Cirf => Ok(cirf_to_tirf(seconds_j2000).transpose()),
+                DynFrame::Tirf => Ok(Rotation::IDENTITY),
+                DynFrame::Itrf => Ok(tirf_to_itrf(centuries_j2000)),
+                DynFrame::BodyFixed(_) => Ok(self
+                    .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                    .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?)),
+            },
+            DynFrame::Itrf => match frame {
+                DynFrame::Icrf => Ok(tirf_to_itrf(centuries_j2000)
+                    .transpose()
+                    .compose(&cirf_to_tirf(seconds_j2000).transpose())
+                    .compose(&icrf_to_cirf(centuries_j2000).transpose())),
+                DynFrame::Cirf => Ok(tirf_to_itrf(centuries_j2000)
+                    .transpose()
+                    .compose(&cirf_to_tirf(seconds_j2000).transpose())),
+                DynFrame::Tirf => Ok(tirf_to_itrf(centuries_j2000).transpose()),
+                DynFrame::Itrf => Ok(Rotation::IDENTITY),
+                DynFrame::BodyFixed(_) => Ok(self
+                    .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                    .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?)),
+            },
+            DynFrame::BodyFixed(origin) => match frame {
+                DynFrame::Icrf => Ok(icrf_to_bodyfixed(origin, seconds_j2000)
+                    .ok_or(format!(
+                        "no rotational parameters available for {}",
+                        origin.name()
+                    ))?
+                    .transpose()),
+                DynFrame::Cirf => Ok(self
+                    .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                    .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?)),
+                DynFrame::Tirf => Ok(self
+                    .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                    .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?)),
+                DynFrame::Itrf => Ok(self
+                    .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                    .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?)),
+                DynFrame::BodyFixed(target) => {
+                    if origin == target {
+                        Ok(Rotation::IDENTITY)
+                    } else {
+                        Ok(self
+                            .try_rotation(&DynFrame::Icrf, time.clone(), provider)?
+                            .compose(&DynFrame::Icrf.try_rotation(frame, time, provider)?))
+                    }
+                }
+            },
         }
     }
 }

crates/lox-orbits/src/frames/iau.rs

@@ -0,0 +1,47 @@
+/*
+ * Copyright (c) 2024. Helge Eichhorn and the LOX contributors
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, you can obtain one at https://mozilla.org/MPL/2.0/.
+ */
+use crate::frames::{BodyFixed, FrameTransformationProvider, Icrf, TryRotateTo};
+use crate::rotations::Rotation;
+use glam::{DMat3, DVec3};
+use lox_bodies::{MaybeRotationalElements, RotationalElements};
+use lox_math::types::units::Seconds;
+use lox_time::prelude::Tdb;
+use lox_time::TimeLike;
+
+pub(crate) fn icrf_to_bodyfixed<T: MaybeRotationalElements>(
+    body: &T,
+    seconds: Seconds,
+) -> Option<Rotation> {
+    let (Some((right_ascension, declination, prime_meridian)), Some(rates)) = (
+        body.maybe_rotational_elements(seconds),
+        body.maybe_rotational_element_rates(seconds),
+    ) else {
+        return None;
+    };
+    let m1 = DMat3::from_rotation_z(-right_ascension);
+    let m2 = DMat3::from_rotation_x(-declination);
+    let m3 = DMat3::from_rotation_z(-prime_meridian);
+    let m = m3 * m2 * m1;
+    let v = DVec3::from(rates);
+    Some(Rotation::new(m).with_angular_velocity(v))
+}
+
+impl<O: RotationalElements, P: FrameTransformationProvider> TryRotateTo<BodyFixed<O>, P> for Icrf {
+    type Error = ();
+
+    fn try_rotation<T: TimeLike>(
+        &self,
+        frame: &BodyFixed<O>,
+        time: T,
+        provider: &P,
+    ) -> Result<Rotation, Self::Error> {
+        // FIXME
+        let seconds = time.seconds_since_j2000();
+        Ok(icrf_to_bodyfixed(&frame.0, seconds).unwrap())
+    }
+}

crates/lox-orbits/src/frames/iers.rs

@@ -5,11 +5,38 @@
  * License, v. 2.0. If a copy of the MPL was not distributed with this
  * file, you can obtain one at https://mozilla.org/MPL/2.0/.
  */
-use crate::frames::{Cirf, Icrf, RotateTo};
 use crate::rotations::Rotation;
+use glam::{DMat3, DVec2, DVec3};
+use lox_bodies::{Earth, RotationalElements};
+use lox_earth::cio::s06::s;
+use lox_earth::cip::xy06::xy;
+use lox_earth::coordinate_transformations::{
+    celestial_to_intermediate_frame_of_date_matrix, polar_motion_matrix,
+};
+use lox_earth::rotation_angle::RotationAngle;
+use lox_earth::tio::sp_00;
+use lox_math::constants::f64::time::SECONDS_PER_DAY;
 
-impl RotateTo<Cirf> for Icrf {
-    fn rotation(&self, _: &Cirf) -> Rotation {
-        todo!()
-    }
+pub fn icrf_to_cirf(centuries: f64) -> Rotation {
+    // TODO: Add IERS corrections
+    let cip_coords = xy(centuries);
+    let cio_locator = s(centuries, cip_coords);
+    let m = celestial_to_intermediate_frame_of_date_matrix(cip_coords, cio_locator);
+    Rotation::new(m)
+}
+
+pub fn cirf_to_tirf(seconds: f64) -> Rotation {
+    let era = Earth::rotation_angle_00(seconds / SECONDS_PER_DAY);
+    let rate = Earth.prime_meridian_dot(seconds);
+    let m = DMat3::from_rotation_z(-era);
+    let v = DVec3::new(0.0, 0.0, rate);
+    Rotation::new(m).with_angular_velocity(v)
+}
+
+pub fn tirf_to_itrf(centuries: f64) -> Rotation {
+    // TODO: Add IERS corrections
+    let pole_coords: DVec2 = (0.0, 0.0).into();
+    let tio_locator = sp_00(centuries);
+    let m = polar_motion_matrix(pole_coords, tio_locator);
+    Rotation::new(m)
 }