Initial work on LRO OD example

ChristopherRabotin · ChristopherRabotin · commit 2b875a6c7575 · 2024-08-14T23:25:26.000-06:00
- Add ability to build Trajectory from BSP
- Include line of sight (eclipsing) computation in OD measurements
diff --git a/Cargo.toml b/Cargo.toml
@@ -126,3 +126,8 @@ doc-scrape-examples = true
 name = "03_geo_sk"
 path = "examples/03_geo_analysis/stationkeeping.rs"
 doc-scrape-examples = true
+
+[[example]]
+name = "04_lro_od"
+path = "examples/04_lro_od/main.rs"
+doc-scrape-examples = true
diff --git a/examples/04_lro_od/README.md b/examples/04_lro_od/README.md
@@ -0,0 +1,12 @@
+# Orbit Determination of the Lunar Reconnaissance Orbiter
+
+Spacecraft operations require high fidelity modeling of the orbital dynamics and high fidelity orbit determination.
+
+In this example, you'll learn how to use an "as-flown" (_definitive_) SPICE BSP ephemeris file to simulate orbit determination measurements from ground stations. Then, you'll learn how to set up an orbit determination process in Nyx with high fidelity Moon dynamics and estimate the state of LRO. Finally, you'll learn how to compare two ephemerides in the Radial, In-track, Cross-track (RIC) frame.
+
+To run this example, just execute:
+```sh
+RUST_LOG=info cargo run --example 04_lro_od --release
+```
+
+Building in `release` mode will make the computation significantly faster. Specifying `RUST_LOG=info` will allow you to see all of the information messages happening in ANISE and Nyx throughout the execution of the program.
diff --git a/examples/04_lro_od/dsn-network.yaml b/examples/04_lro_od/dsn-network.yaml
@@ -0,0 +1,59 @@
+- name: DSS-65 Madrid
+  frame:
+    ephemeris_id: 399
+    orientation_id: 399
+    mu_km3_s2: 398600.435436096
+    shape: null
+  elevation_mask_deg: 5.0
+  range_noise_km:
+    bias:
+      tau: 24 h
+      process_noise: 5.0e-3 # 5 m
+  doppler_noise_km_s:
+    bias:
+      tau: 24 h
+      process_noise: 50.0e-6 # 5 cm/s
+  light_time_correction: false
+  latitude_deg: 40.427222
+  longitude_deg: 4.250556
+  height_km: 0.834939
+
+- name: DSS-34 Canberra
+  frame:
+    ephemeris_id: 399
+    orientation_id: 399
+    mu_km3_s2: 398600.435436096
+    shape: null
+  latitude_deg: -35.398333
+  longitude_deg: 148.981944
+  height_km: 0.691750
+  elevation_mask_deg: 5.0
+  range_noise_km:
+    bias:
+      tau: 24 h
+      process_noise: 5.0e-3 # 5 m
+  doppler_noise_km_s:
+    bias:
+      tau: 24 h
+      process_noise: 50.0e-6 # 5 cm/s
+  light_time_correction: false
+
+- name: DSS-13 Goldstone
+  frame:
+    ephemeris_id: 399
+    orientation_id: 399
+    mu_km3_s2: 398600.435436096
+    shape: null
+  latitude_deg: 35.247164
+  longitude_deg: 243.205
+  height_km: 1.071149
+  elevation_mask_deg: 5.0
+  range_noise_km:
+    bias:
+      tau: 24 h
+      process_noise: 5.0e-3 # 5 m
+  doppler_noise_km_s:
+    bias:
+      tau: 24 h
+      process_noise: 50.0e-6 # 5 cm/s
+  light_time_correction: false
diff --git a/examples/04_lro_od/gfx/LRO-ANISE-GUI.png b/examples/04_lro_od/gfx/LRO-ANISE-GUI.png
diff --git a/examples/04_lro_od/main.rs b/examples/04_lro_od/main.rs
@@ -0,0 +1,131 @@
+#![doc = include_str!("./README.md")]
+extern crate log;
+extern crate nyx_space as nyx;
+extern crate pretty_env_logger as pel;
+
+use anise::{
+    almanac::metaload::MetaFile,
+    constants::{
+        celestial_objects::{JUPITER_BARYCENTER, MOON, SUN},
+        frames::{EARTH_J2000, MOON_J2000},
+    },
+};
+use hifitime::{Epoch, TimeUnits, Unit};
+use nyx::{
+    cosmic::{eclipse::EclipseLocator, Aberration, Frame, MetaAlmanac, SrpConfig},
+    dynamics::{guidance::LocalFrame, OrbitalDynamics, SolarPressure, SpacecraftDynamics},
+    io::{ConfigRepr, ExportCfg},
+    mc::MonteCarlo,
+    md::prelude::Traj,
+    od::{
+        msr::RangeDoppler,
+        prelude::{TrackingArcSim, TrkConfig, KF},
+        process::SpacecraftUncertainty,
+        GroundStation, SpacecraftODProcess,
+    },
+    propagators::Propagator,
+    Orbit, Spacecraft, State,
+};
+
+use std::{collections::BTreeMap, error::Error, path::PathBuf, str::FromStr, sync::Arc};
+
+fn main() -> Result<(), Box<dyn Error>> {
+    pel::init();
+    // Dynamics models require planetary constants and ephemerides to be defined.
+    // Let's start by grabbing those by using ANISE's latest MetaAlmanac.
+    // For details, refer to https://github.com/nyx-space/anise/blob/master/data/latest.dhall.
+
+    // Download the latest (of time of writing) LRO definitive ephemeris from the public Nyx Space cloud.
+    // Note that the original file is in _big endian_ format, and my machine is little endian, so I've used the
+    // `bingo` tool from https://naif.jpl.nasa.gov/naif/utilities_PC_Linux_64bit.html to convert the original file
+    // to little endian and upload it to the cloud.
+    // Refer to https://naif.jpl.nasa.gov/pub/naif/pds/data/lro-l-spice-6-v1.0/lrosp_1000/data/spk/?C=M;O=D for original file.
+    let mut lro_def_ephem = MetaFile {
+        uri: "http://public-data.nyxspace.com/nyx/examples/lrorg_2023349_2024075_v01_LE.bsp"
+            .to_string(),
+        crc32: Some(0xe76ce3b5),
+    };
+    lro_def_ephem.process()?;
+
+    // Load this ephem in the general Almanac we're using for this analysis.
+    let almanac = Arc::new(
+        MetaAlmanac::latest()
+            .map_err(Box::new)?
+            .load_from_metafile(lro_def_ephem)?,
+    );
+
+    // Orbit determination requires a Trajectory structure, which can be saved as parquet file.
+    // In our case, the trajectory comes from the BSP file, so we need to build a Trajectory from the almanac directly.
+    // To query the Almanac, we need to build the LRO frame in the J2000 orientation in our case.
+    // Inspecting the LRO BSP in the ANISE GUI shows us that NASA has assigned ID -85 to LRO.
+    let lro_frame = Frame::from_ephem_j2000(-85);
+    // To build the trajectory we need to provide a spacecraft template.
+    let sc_template = Spacecraft::builder()
+        .dry_mass_kg(1018.0) // Launch masses
+        .fuel_mass_kg(900.0)
+        .srp(SrpConfig {
+            // SRP configuration is arbitrary, but we will be estimating it anyway.
+            area_m2: 3.9 * 2.7,
+            cr: 0.9,
+        })
+        .orbit(Orbit::zero(MOON_J2000)) // Setting a zero orbit here because it's just a template
+        .build();
+    // Now we can build the trajectory from the BSP file.
+    // We'll arbitrarily set the tracking arc to 48 hours with a one minute time step.
+    let trajectory = Traj::from_bsp(
+        lro_frame,
+        MOON_J2000,
+        almanac.clone(),
+        sc_template,
+        5.seconds(),
+        Some(Epoch::from_str("2024-01-01 00:00:00 UTC").unwrap()),
+        Some(Epoch::from_str("2024-01-03 00:00:00 UTC").unwrap()),
+        Aberration::LT,
+        Some("LRO".to_string()),
+    )?;
+
+    println!("{trajectory}");
+
+    // Load the Deep Space Network ground stations.
+    // Nyx allows you to build these at runtime but it's pretty static so we can just load them from YAML.
+    let ground_station_file: PathBuf = [
+        env!("CARGO_MANIFEST_DIR"),
+        "examples",
+        "04_lro_od",
+        "dsn-network.yaml",
+    ]
+    .iter()
+    .collect();
+
+    let devices = GroundStation::load_many(ground_station_file).unwrap();
+
+    // Typical OD software requires that you specify your own tracking schedule or you'll have overlapping measurements.
+    // Nyx can build a tracking schedule for you based on the first station with access.
+    let trkconfg_yaml: PathBuf = [
+        env!("CARGO_MANIFEST_DIR"),
+        "examples",
+        "04_lro_od",
+        "tracking-cfg.yaml",
+    ]
+    .iter()
+    .collect();
+
+    let configs: BTreeMap<String, TrkConfig> = TrkConfig::load_named(trkconfg_yaml).unwrap();
+
+    dbg!(&configs);
+
+    // Build the tracking arc simulation to generate a "standard measurement".
+    let mut trk = TrackingArcSim::<Spacecraft, RangeDoppler, _>::with_seed(
+        devices, trajectory, configs, 12345,
+    )
+    .unwrap();
+
+    trk.build_schedule(almanac.clone()).unwrap();
+    let arc = trk.generate_measurements(almanac).unwrap();
+    // Save the simulated tracking data
+    arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
+
+    println!("{arc}");
+
+    Ok(())
+}
diff --git a/examples/04_lro_od/tracking-cfg.yaml b/examples/04_lro_od/tracking-cfg.yaml
@@ -0,0 +1,23 @@
+DSS-65 Madrid:
+  scheduler:
+    handoff: Eager
+    cadence: Continuous
+    min_samples: 10
+    sample_alignment: 10 s
+  sampling: 1 min
+
+DSS-34 Canberra:
+  scheduler:
+    handoff: Eager
+    cadence: Continuous
+    min_samples: 10
+    sample_alignment: 10 s
+  sampling: 1 min
+
+DSS-13 Goldstone:
+  scheduler:
+    handoff: Eager
+    cadence: Continuous
+    min_samples: 10
+    sample_alignment: 10 s
+  sampling: 1 min
diff --git a/src/md/trajectory/sc_traj.rs b/src/md/trajectory/sc_traj.rs
@@ -16,8 +16,11 @@
     along with this program.  If not, see <https://www.gnu.org/licenses/>.
 */
 
+use anise::astro::Aberration;
 use anise::constants::orientations::J2000;
+use anise::errors::AlmanacError;
 use anise::prelude::{Almanac, Frame, Orbit};
+use hifitime::TimeSeries;
 use snafu::ResultExt;
 
 use super::TrajError;
@@ -40,6 +43,41 @@ use std::sync::Arc;
 use std::time::Instant;
 
 impl Traj<Spacecraft> {
+    /// Builds a new trajectory built from the SPICE BSP (SPK) file loaded in the provided Almanac, provided the start and stop epochs.
+    ///
+    /// If the start and stop epochs are not provided, then the full domain of the trajectory will be used.
+    pub fn from_bsp(
+        target_frame: Frame,
+        observer_frame: Frame,
+        almanac: Arc<Almanac>,
+        sc_template: Spacecraft,
+        step: Duration,
+        start_epoch: Option<Epoch>,
+        end_epoch: Option<Epoch>,
+        ab_corr: Option<Aberration>,
+        name: Option<String>,
+    ) -> Result<Self, AlmanacError> {
+        let (domain_start, domain_end) =
+            almanac
+                .spk_domain(target_frame.ephemeris_id)
+                .map_err(|e| AlmanacError::Ephemeris {
+                    action: "could not fetch domain",
+                    source: Box::new(e),
+                })?;
+
+        let start_epoch = start_epoch.unwrap_or(domain_start);
+        let end_epoch = end_epoch.unwrap_or(domain_end);
+
+        let time_series = TimeSeries::inclusive(start_epoch, end_epoch, step);
+        let mut states = Vec::with_capacity(time_series.len());
+        for epoch in time_series {
+            let orbit = almanac.transform(target_frame, observer_frame, epoch, ab_corr)?;
+
+            states.push(sc_template.with_orbit(orbit));
+        }
+
+        Ok(Self { name, states })
+    }
     /// Allows converting the source trajectory into the (almost) equivalent trajectory in another frame
     #[allow(clippy::map_clone)]
     pub fn to_frame(&self, new_frame: Frame, almanac: Arc<Almanac>) -> Result<Self, NyxError> {
@@ -129,11 +167,6 @@ impl Traj<Spacecraft> {
         traj.to_parquet(path, events, cfg, almanac)
     }
 
-    /// Convert this spacecraft trajectory into an Orbit trajectory, loosing all references to the spacecraft
-    pub fn downcast(&self) -> Self {
-        unimplemented!()
-    }
-
     /// Initialize a new spacecraft trajectory from the path to a CCSDS OEM file.
     ///
     /// CCSDS OEM only contains the orbit information but Nyx builds spacecraft trajectories.
diff --git a/src/md/trajectory/traj.rs b/src/md/trajectory/traj.rs
@@ -599,7 +599,7 @@ where
                 f,
                 "Trajectory {}in {} from {} to {} ({}, or {:.3} s) [{} states]",
                 match &self.name {
-                    Some(name) => format!("of {name}"),
+                    Some(name) => format!("of {name} "),
                     None => String::new(),
                 },
                 self.first().frame(),
diff --git a/src/od/ground_station.rs b/src/od/ground_station.rs
@@ -22,7 +22,8 @@ use anise::prelude::{Almanac, Frame, Orbit};
 
 use super::msr::RangeDoppler;
 use super::noise::StochasticNoise;
-use super::{ODAlmanacSnafu, ODError, ODTrajSnafu, TrackingDeviceSim};
+use super::{ODAlmanacSnafu, ODError, ODPlanetaryDataSnafu, ODTrajSnafu, TrackingDeviceSim};
+use crate::cosmic::eclipse::{line_of_sight, EclipseState};
 use crate::errors::EventError;
 use crate::io::ConfigRepr;
 use crate::md::prelude::{Interpolatable, Traj};
@@ -157,11 +158,9 @@ impl GroundStation {
     }
 
     /// Computes the azimuth and elevation of the provided object seen from this ground station, both in degrees.
-    /// Also returns the ground station's orbit in the frame of the receiver
+    /// This is a shortcut to almanac.azimuth_elevation_range_sez.
     pub fn azimuth_elevation_of(&self, rx: Orbit, almanac: &Almanac) -> AlmanacResult<AzElRange> {
-        almanac
-            .clone()
-            .azimuth_elevation_range_sez(rx, self.to_orbit(rx.epoch, almanac).unwrap())
+        almanac.azimuth_elevation_range_sez(rx, self.to_orbit(rx.epoch, almanac).unwrap())
     }
 
     /// Return this ground station as an orbit in its current frame
@@ -256,6 +255,27 @@ impl TrackingDeviceSim<Spacecraft, RangeDoppler> for GroundStation {
                     return Ok(None);
                 }
 
+                // If the frame of the trajectory is different from that of the ground station, then check that there is no eclipse.
+                if !self.frame.ephem_origin_match(rx_0.frame()) {
+                    let observer = self.to_orbit(rx_0.orbit.epoch, &almanac).unwrap();
+                    if line_of_sight(
+                        observer,
+                        rx_0.orbit,
+                        almanac
+                            .frame_from_uid(rx_0.frame())
+                            .context(ODPlanetaryDataSnafu {
+                                action: "computing line of sight",
+                            })?,
+                        &almanac,
+                    )
+                    .context(ODAlmanacSnafu {
+                        action: "computing line of sight",
+                    })? == EclipseState::Umbra
+                    {
+                        return Ok(None);
+                    }
+                }
+
                 // Noises are computed at the midpoint of the integration time.
                 let (timestamp_noise_s, range_noise_km, doppler_noise_km_s) =
                     self.noises(epoch - integration_time * 0.5, rng)?;
diff --git a/src/od/mod.rs b/src/od/mod.rs
@@ -26,6 +26,7 @@ use crate::propagators::PropagationError;
 use crate::time::Epoch;
 use crate::Orbit;
 pub use crate::{State, TimeTagged};
+use anise::almanac::planetary::PlanetaryDataError;
 use anise::errors::AlmanacError;
 use hifitime::Duration;
 use snafu::prelude::Snafu;
@@ -204,6 +205,12 @@ pub enum ODError {
         source: Box<AlmanacError>,
         action: &'static str,
     },
+    #[snafu(display("OD failed due to planetary data in Almanac: {action} {source}"))]
+    ODPlanetaryData {
+        #[snafu(source(from(PlanetaryDataError, Box::new)))]
+        source: Box<PlanetaryDataError>,
+        action: &'static str,
+    },
     #[snafu(display("not enough residuals to {action}"))]
     ODNoResiduals { action: &'static str },
 }
