From fc5fa3d83327b5f35a83d53689f59f99478a853e Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sat, 27 Jan 2024 11:14:55 -0700 Subject: [PATCH 1/4] Initial implementation sun_angle for sun probe earth angle calcs --- anise/src/almanac/mod.rs | 1 + anise/src/almanac/solar.rs | 91 ++++++++++++++++++++++++++++++++++++++ 2 files changed, 92 insertions(+) create mode 100644 anise/src/almanac/solar.rs diff --git a/anise/src/almanac/mod.rs b/anise/src/almanac/mod.rs index 447df7ba..e3a6fe8c 100644 --- a/anise/src/almanac/mod.rs +++ b/anise/src/almanac/mod.rs @@ -38,6 +38,7 @@ pub const MAX_PLANETARY_DATA: usize = 64; pub mod aer; pub mod bpc; pub mod planetary; +pub mod solar; pub mod spk; pub mod transform; diff --git a/anise/src/almanac/solar.rs b/anise/src/almanac/solar.rs new file mode 100644 index 00000000..aabb6716 --- /dev/null +++ b/anise/src/almanac/solar.rs @@ -0,0 +1,91 @@ +/* + * ANISE Toolkit + * Copyright (C) 2021-2023 Christopher Rabotin et al. (cf. AUTHORS.md) + * 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/. + * + * Documentation: https://nyxspace.com/ + */ + +use crate::{constants::frames::SUN_J2000, ephemerides::EphemerisError, prelude::Frame, NaifId}; + +use super::Almanac; + +use hifitime::Epoch; + +#[cfg(feature = "python")] +use pyo3::prelude::*; + +#[cfg_attr(feature = "python", pymethods)] +impl Almanac { + /// Returns the angle (between 0 and 180 degrees) between the observer and the Sun, and the observer and the target body ID. + /// This computes the Sun Probe Earth angle (SPE) if the probe is in a loaded, its ID is the "observer_id", and the target is set to its central body. + /// + /// # Geometry + /// If the SPE is greater than 90 degrees, then the celestial object below the probe is in sunlight. + /// + /// ## Sunrise at nadir + /// ```text + /// Sun + /// | \ + /// | \ + /// | \ + /// Obs. -- Target + /// ``` + /// ## Sun high at nadir + /// ```text + /// Sun + /// \ + /// \ __ θ > 90 + /// \ \ + /// Obs. ---------- Target + /// ``` + /// + /// ## Sunset at nadir + /// ```text + /// Sun + /// / + /// / __ θ < 90 + /// / / + /// Obs. -- Target + /// ``` + /// + /// # Algorithm + /// 1. Compute the position of the Sun as seen from the observer + /// 2. Compute the position of the target as seen from the observer + /// 3. Return the arccosine of the dot product of the norms of these vectors. + pub fn sun_angle_deg( + &self, + target_id: NaifId, + observer_id: NaifId, + epoch: Epoch, + ) -> Result { + let obs_to_sun = + self.translate_geometric(SUN_J2000, Frame::from_ephem_j2000(observer_id), epoch)?; + let obs_to_target = self.translate_geometric( + Frame::from_ephem_j2000(observer_id), + Frame::from_ephem_j2000(target_id), + epoch, + )?; + + Ok(obs_to_sun + .r_hat() + .dot(&obs_to_target.r_hat()) + .acos() + .to_degrees()) + } + + /// Convenience function that calls `sun_angle_deg` with the provided frames instead of the ephemeris ID. + pub fn sun_angle_deg_from_frame( + &self, + target: Frame, + observer: Frame, + epoch: Epoch, + ) -> Result { + self.sun_angle_deg(target.ephemeris_id, observer.ephemeris_id, epoch) + } +} + +#[cfg(test)] +mod ut_solar {} From 9a425479eddc5c84e3ed9145b50c37ff018dd117 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sat, 27 Jan 2024 15:16:31 -0700 Subject: [PATCH 2/4] Add test for sun_angle (SPE) calculation Next up: a Python tutorial --- anise/src/almanac/solar.rs | 86 +++++++++++++++++++++++++++++++++++++- 1 file changed, 84 insertions(+), 2 deletions(-) diff --git a/anise/src/almanac/solar.rs b/anise/src/almanac/solar.rs index aabb6716..6898d7db 100644 --- a/anise/src/almanac/solar.rs +++ b/anise/src/almanac/solar.rs @@ -64,8 +64,8 @@ impl Almanac { let obs_to_sun = self.translate_geometric(SUN_J2000, Frame::from_ephem_j2000(observer_id), epoch)?; let obs_to_target = self.translate_geometric( - Frame::from_ephem_j2000(observer_id), Frame::from_ephem_j2000(target_id), + Frame::from_ephem_j2000(observer_id), epoch, )?; @@ -88,4 +88,86 @@ impl Almanac { } #[cfg(test)] -mod ut_solar {} +mod ut_solar { + use crate::{ + constants::{ + celestial_objects::EARTH, + frames::{EARTH_J2000, IAU_EARTH_FRAME, SUN_J2000}, + usual_planetary_constants::MEAN_EARTH_ANGULAR_VELOCITY_DEG_S, + }, + prelude::*, + }; + + /// Load a BSP of a spacecraft trajectory, compute the sun elevation at different points on the surface below that spacecraft + /// and ensure that it matches with the SPE calculation. + #[test] + fn verify_geometry() { + let ctx = Almanac::default() + .load("../data/de440s.bsp") + .and_then(|ctx| ctx.load("../data/gmat-hermite.bsp")) + .and_then(|ctx| ctx.load("../data/pck11.pca")) + .unwrap(); + + let epoch = Epoch::from_gregorian_hms(2000, 1, 1, 12, 0, 0, TimeScale::UTC); + + let sc_id = -10000001; + + let my_sc_j2k = Frame::from_ephem_j2000(sc_id); + + // Grab the state in the J2000 frame + let state = ctx.transform(my_sc_j2k, EARTH_J2000, epoch, None).unwrap(); + + // We'll check at four different points in the orbit + for epoch in TimeSeries::inclusive( + epoch, + epoch + state.period().unwrap(), + 0.05 * state.period().unwrap(), + ) { + let spe_deg = ctx.sun_angle_deg(EARTH, sc_id, epoch).unwrap(); + assert_eq!( + spe_deg, + ctx.sun_angle_deg_from_frame(EARTH_J2000, my_sc_j2k, epoch) + .unwrap() + ); + + let iau_earth = ctx.frame_from_uid(IAU_EARTH_FRAME).unwrap(); + + // Compute this state in the body fixed frame. + let state_bf = ctx + .transform(my_sc_j2k, IAU_EARTH_FRAME, epoch, None) + .unwrap(); + // Build a local point on the surface below this spacecraft + let nadir_surface_point = Orbit::try_latlongalt( + state_bf.latitude_deg().unwrap(), + state_bf.longitude_deg(), + 0.0, + MEAN_EARTH_ANGULAR_VELOCITY_DEG_S, + epoch, + iau_earth, + ) + .unwrap(); + // Fetch the state of the sun at this time. + let sun_state = ctx + .transform(SUN_J2000, IAU_EARTH_FRAME, epoch, None) + .unwrap(); + + // Compute the Sun elevation from that point. + let sun_elevation_deg = ctx + .azimuth_elevation_range_sez(sun_state, nadir_surface_point) + .unwrap() + .elevation_deg; + + println!( + "{epoch}\tsun el = {sun_elevation_deg:.3} deg\tlat = {:.3} deg\tlong = {:.3} deg\t-->\tSPE = {spe_deg:.3} deg", + nadir_surface_point.latitude_deg().unwrap(), + nadir_surface_point.longitude_deg() + ); + + // Test: the sun elevation from the ground should be about 90 degrees _less_ than the SPE + // The "about" is because the SPE effectively assumes a spherical celestial object, but the frame + // from the Almanac accounts for the flattening when computing the exact location below the vehicle. + // Hence, there is a small difference (we accept up to 0.05 degrees of difference). + assert!((sun_elevation_deg + 90.0 - spe_deg).abs() < 5e-2) + } + } +} From 2944c0e7b0a79c909177d1834615f603859ed66b Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 28 Jan 2024 18:34:01 -0700 Subject: [PATCH 3/4] Add validity domain via `spk_domain` and `bpc_domain`, also available in Python Fix #179 --- anise/src/almanac/bpc.rs | 52 ++++++++++++++++- anise/src/almanac/spk.rs | 58 ++++++++++++++++++- anise/src/astro/orbit.rs | 2 +- anise/src/constants.rs | 10 ++-- anise/src/ephemerides/translations.rs | 4 +- anise/src/math/rotation/quaternion.rs | 2 +- anise/src/naif/daf/data_types.rs | 5 ++ anise/src/naif/pck/mod.rs | 20 +++++++ anise/src/naif/spk/summary.rs | 35 +++++++---- anise/src/orientations/rotations.rs | 2 +- anise/src/structure/planetocentric/mod.rs | 2 +- .../ephemerides/parent_translation_verif.rs | 16 ++++- 12 files changed, 181 insertions(+), 27 deletions(-) diff --git a/anise/src/almanac/bpc.rs b/anise/src/almanac/bpc.rs index 22b9b34e..bee4ab7c 100644 --- a/anise/src/almanac/bpc.rs +++ b/anise/src/almanac/bpc.rs @@ -10,10 +10,14 @@ use hifitime::Epoch; -use crate::naif::daf::DAFError; +#[cfg(feature = "python")] +use pyo3::prelude::*; + +use crate::naif::daf::NAIFSummaryRecord; use crate::naif::pck::BPCSummaryRecord; use crate::naif::BPC; use crate::orientations::OrientationError; +use crate::{naif::daf::DAFError, NaifId}; use super::{Almanac, MAX_LOADED_BPCS}; @@ -172,6 +176,52 @@ impl Almanac { } } +#[cfg_attr(feature = "python", pymethods)] +impl Almanac { + /// Returns a vector of the summaries whose ID matches the desired `id`, in the order in which they will be used, i.e. in reverse loading order. + pub fn bpc_summaries(&self, id: NaifId) -> Result, OrientationError> { + let mut summaries = vec![]; + + for maybe_bpc in self.bpc_data.iter().take(self.num_loaded_bpc()).rev() { + let bpc = maybe_bpc.as_ref().unwrap(); + if let Ok((summary, _)) = bpc.summary_from_id(id) { + // NOTE: We're iterating backward, so the correct BPC number is "total loaded" minus "current iteration". + summaries.push(*summary); + } + } + + if summaries.is_empty() { + // If we're reached this point, there is no relevant summary + Err(OrientationError::BPC { + action: "searching for BPC summary", + source: DAFError::SummaryIdError { kind: "BPC", id }, + }) + } else { + Ok(summaries) + } + } + + /// Returns the applicable domain of the request id, i.e. start and end epoch that the provided id has loaded data. + pub fn bpc_domain(&self, id: NaifId) -> Result<(Epoch, Epoch), OrientationError> { + let summaries = self.bpc_summaries(id)?; + + // We know that the summaries is non-empty because if it is, the previous function call returns an error. + let start = summaries + .iter() + .min_by_key(|summary| summary.start_epoch()) + .unwrap() + .start_epoch(); + + let end = summaries + .iter() + .max_by_key(|summary| summary.end_epoch()) + .unwrap() + .end_epoch(); + + Ok((start, end)) + } +} + #[cfg(test)] mod ut_almanac_bpc { use crate::prelude::{Almanac, Epoch}; diff --git a/anise/src/almanac/spk.rs b/anise/src/almanac/spk.rs index 824cb822..eee13c6f 100644 --- a/anise/src/almanac/spk.rs +++ b/anise/src/almanac/spk.rs @@ -10,10 +10,14 @@ use hifitime::Epoch; -use crate::ephemerides::EphemerisError; +#[cfg(feature = "python")] +use pyo3::prelude::*; + use crate::naif::daf::DAFError; +use crate::naif::daf::NAIFSummaryRecord; use crate::naif::spk::summary::SPKSummaryRecord; use crate::naif::SPK; +use crate::{ephemerides::EphemerisError, NaifId}; use log::error; use super::{Almanac, MAX_LOADED_SPKS}; @@ -45,7 +49,9 @@ impl Almanac { me.spk_data[data_idx] = Some(spk); Ok(me) } +} +impl Almanac { pub fn num_loaded_spk(&self) -> usize { let mut count = 0; for maybe in &self.spk_data { @@ -122,7 +128,7 @@ impl Almanac { }) } - /// Returns the summary given the name of the summary record. + /// Returns the most recently loaded summary by its name, if any with that ID are available pub fn spk_summary_from_name( &self, name: &str, @@ -152,7 +158,7 @@ impl Almanac { }) } - /// Returns the summary given the name of the summary record if that summary has data defined at the requested epoch + /// Returns the most recently loaded summary by its ID, if any with that ID are available pub fn spk_summary( &self, id: i32, @@ -180,6 +186,52 @@ impl Almanac { } } +#[cfg_attr(feature = "python", pymethods)] +impl Almanac { + /// Returns a vector of the summaries whose ID matches the desired `id`, in the order in which they will be used, i.e. in reverse loading order. + pub fn spk_summaries(&self, id: NaifId) -> Result, EphemerisError> { + let mut summaries = vec![]; + for maybe_spk in self.spk_data.iter().take(self.num_loaded_spk()).rev() { + let spk = maybe_spk.as_ref().unwrap(); + if let Ok((summary, _)) = spk.summary_from_id(id) { + // NOTE: We're iterating backward, so the correct SPK number is "total loaded" minus "current iteration". + summaries.push(*summary); + } + } + + if summaries.is_empty() { + error!("Almanac: No summary {id} valid"); + // If we're reached this point, there is no relevant summary + Err(EphemerisError::SPK { + action: "searching for SPK summary", + source: DAFError::SummaryIdError { kind: "SPK", id }, + }) + } else { + Ok(summaries) + } + } + + /// Returns the applicable domain of the request id, i.e. start and end epoch that the provided id has loaded data. + pub fn spk_domain(&self, id: NaifId) -> Result<(Epoch, Epoch), EphemerisError> { + let summaries = self.spk_summaries(id)?; + + // We know that the summaries is non-empty because if it is, the previous function call returns an error. + let start = summaries + .iter() + .min_by_key(|summary| summary.start_epoch()) + .unwrap() + .start_epoch(); + + let end = summaries + .iter() + .max_by_key(|summary| summary.end_epoch()) + .unwrap() + .end_epoch(); + + Ok((start, end)) + } +} + #[cfg(test)] mod ut_almanac_spk { use crate::{ diff --git a/anise/src/astro/orbit.rs b/anise/src/astro/orbit.rs index fbed1838..53bf1307 100644 --- a/anise/src/astro/orbit.rs +++ b/anise/src/astro/orbit.rs @@ -675,7 +675,7 @@ impl Orbit { /// Returns the true anomaly in degrees between 0 and 360.0 /// /// NOTE: This function will emit a warning stating that the TA should be avoided if in a very near circular orbit - /// Code from https://github.com/ChristopherRabotin/GMAT/blob/80bde040e12946a61dae90d9fc3538f16df34190/src/gmatutil/util/StateConversionUtil.cpp#L6835 + /// Code from /// /// LIMITATION: For an orbit whose true anomaly is (very nearly) 0.0 or 180.0, this function may return either 0.0 or 180.0 with a very small time increment. /// This is due to the precision of the cosine calculation: if the arccosine calculation is out of bounds, the sign of the cosine of the true anomaly is used diff --git a/anise/src/constants.rs b/anise/src/constants.rs index bc985717..315e5e98 100644 --- a/anise/src/constants.rs +++ b/anise/src/constants.rs @@ -89,7 +89,7 @@ pub mod orientations { /// The rotation from B1950 to J2000 is /// \[ -z \] \[ theta \] \[ -zeta \] /// 3 2 3 - /// The values for z, theta, and zeta are computed from the formulas given in table 5 of [5]. + /// The values for z, theta, and zeta are computed from the formulas given in table 5 of \[5\]. /// z = 1153.04066200330" /// theta = 1002.26108439117" /// zeta = 1152.84248596724" @@ -105,7 +105,7 @@ pub mod orientations { /// \[ 0.53155" \] /// 3 /// - /// In [2], Standish uses two separate rotations, + /// In \[2\], Standish uses two separate rotations, /// /// \[ 0.00073" \] P \[ 0.5316" \] /// 3 3 @@ -204,7 +204,7 @@ pub mod orientations { pub const IAU_NEPTUNE: NaifId = 799; pub const IAU_URANUS: NaifId = 899; - /// Angle between J2000 to solar system ecliptic J2000 ([ECLIPJ2000]), in radians (about 23.43929 degrees). Apply this rotation about the X axis ([r1]) + /// Angle between J2000 to solar system ecliptic J2000 ([ECLIPJ2000]), in radians (about 23.43929 degrees). Apply this rotation about the X axis (R1) pub const J2000_TO_ECLIPJ2000_ANGLE_RAD: f64 = 0.40909280422232897; /// Given the frame ID, try to return a human name @@ -280,7 +280,7 @@ pub mod frames { /// Typical planetary constants that aren't found in SPICE input files. pub mod usual_planetary_constants { /// Mean angular velocity of the Earth in deg/s - /// Source: G. Xu and Y. Xu, "GPS", DOI 10.1007/978-3-662-50367-6_2, 2016 (confirmed by https://hpiers.obspm.fr/eop-pc/models/constants.html) + /// Source: G. Xu and Y. Xu, "GPS", DOI 10.1007/978-3-662-50367-6_2, 2016 (confirmed by ) pub const MEAN_EARTH_ANGULAR_VELOCITY_DEG_S: f64 = 0.004178079012116429; /// Mean angular velocity of the Moon in deg/s, computed from hifitime: /// ```py @@ -288,7 +288,7 @@ pub mod usual_planetary_constants { /// >>> tau/moon_period.to_seconds() /// 2.661698975163682e-06 /// ``` - /// Source: https://www.britannica.com/science/month#ref225844 via https://en.wikipedia.org/w/index.php?title=Lunar_day&oldid=1180701337 + /// Source: via pub const MEAN_MOON_ANGULAR_VELOCITY_DEG_S: f64 = 2.661_698_975_163_682e-6; } diff --git a/anise/src/ephemerides/translations.rs b/anise/src/ephemerides/translations.rs index ec21fcae..c476fbf8 100644 --- a/anise/src/ephemerides/translations.rs +++ b/anise/src/ephemerides/translations.rs @@ -188,7 +188,7 @@ impl Almanac { /// Translates the provided Cartesian state into the requested observer frame /// - /// **WARNING:** This function only performs the translation and no rotation _whatsoever_. Use the [transform_to] function instead to include rotations. + /// **WARNING:** This function only performs the translation and no rotation _whatsoever_. Use the `transform_to` function instead to include rotations. #[allow(clippy::too_many_arguments)] pub fn translate_to( &self, @@ -205,7 +205,7 @@ impl Almanac { impl Almanac { /// Translates a state with its origin (`to_frame`) and given its units (distance_unit, time_unit), returns that state with respect to the requested frame /// - /// **WARNING:** This function only performs the translation and no rotation _whatsoever_. Use the [transform_state_to] function instead to include rotations. + /// **WARNING:** This function only performs the translation and no rotation _whatsoever_. Use the `transform_state_to` function instead to include rotations. #[allow(clippy::too_many_arguments)] pub fn translate_state_to( &self, diff --git a/anise/src/math/rotation/quaternion.rs b/anise/src/math/rotation/quaternion.rs index 2eb01bb9..8089ae00 100644 --- a/anise/src/math/rotation/quaternion.rs +++ b/anise/src/math/rotation/quaternion.rs @@ -20,7 +20,7 @@ use snafu::ensure; use super::EPSILON_RAD; -/// Quaternion will always be a unit quaternion in ANISE, cf. [EulerParameter]. +/// Quaternion will always be a unit quaternion in ANISE, cf. EulerParameter. /// /// In ANISE, Quaternions use exclusively the Hamiltonian convenstion. pub type Quaternion = EulerParameter; diff --git a/anise/src/naif/daf/data_types.rs b/anise/src/naif/daf/data_types.rs index 1753c1f0..69b451d4 100644 --- a/anise/src/naif/daf/data_types.rs +++ b/anise/src/naif/daf/data_types.rs @@ -16,8 +16,13 @@ use crate::naif::daf::DatatypeSnafu; use super::DAFError; +#[cfg(feature = "python")] +use pyo3::prelude::*; + +#[cfg_attr(feature = "python", pyclass)] #[derive(Copy, Clone, Debug, PartialEq)] #[repr(u8)] + pub enum DataType { Type1ModifiedDifferenceArray = 1, Type2ChebyshevTriplet = 2, diff --git a/anise/src/naif/pck/mod.rs b/anise/src/naif/pck/mod.rs index 6c13b14f..827779a0 100644 --- a/anise/src/naif/pck/mod.rs +++ b/anise/src/naif/pck/mod.rs @@ -15,8 +15,13 @@ use crate::{ use hifitime::Epoch; use zerocopy::{AsBytes, FromBytes, FromZeroes}; +#[cfg(feature = "python")] +use pyo3::prelude::*; + use super::daf::DafDataType; +#[cfg_attr(feature = "python", pyclass)] +#[cfg_attr(feature = "python", pyo3(module = "anise.internals"))] #[derive(Clone, Copy, Debug, Default, AsBytes, FromZeroes, FromBytes)] #[repr(C)] pub struct BPCSummaryRecord { @@ -39,6 +44,21 @@ impl BPCSummaryRecord { } } +#[cfg(feature = "python")] +#[pymethods] +impl BPCSummaryRecord { + /// Returns the start epoch of this BPC Summary + pub fn start_epoch(&self) -> Epoch { + ::start_epoch(self) + } + + /// Returns the start epoch of this BPC Summary + + pub fn end_epoch(&self) -> Epoch { + ::end_epoch(self) + } +} + impl NAIFRecord for BPCSummaryRecord {} impl NAIFSummaryRecord for BPCSummaryRecord { diff --git a/anise/src/naif/spk/summary.rs b/anise/src/naif/spk/summary.rs index f6bb9e67..ebdf944c 100644 --- a/anise/src/naif/spk/summary.rs +++ b/anise/src/naif/spk/summary.rs @@ -12,12 +12,17 @@ use core::fmt; use hifitime::Epoch; use zerocopy::{AsBytes, FromBytes, FromZeroes}; +#[cfg(feature = "python")] +use pyo3::prelude::*; + use crate::{ ephemerides::EphemerisError, naif::daf::{DafDataType, NAIFRecord, NAIFSummaryRecord}, prelude::{Frame, FrameUid}, }; +#[cfg_attr(feature = "python", pyclass)] +#[cfg_attr(feature = "python", pyo3(module = "anise.internals"))] #[derive(Clone, Copy, Debug, Default, AsBytes, FromZeroes, FromBytes, PartialEq)] #[repr(C)] pub struct SPKSummaryRecord { @@ -32,13 +37,6 @@ pub struct SPKSummaryRecord { } impl SPKSummaryRecord { - pub fn data_type(&self) -> Result { - DafDataType::try_from(self.data_type_i).map_err(|source| EphemerisError::SPK { - action: "converting data type from i32", - source, - }) - } - /// Returns the target frame UID of this summary pub fn target_frame_uid(&self) -> FrameUid { FrameUid { @@ -54,6 +52,16 @@ impl SPKSummaryRecord { orientation_id: self.frame_id, } } +} + +#[cfg_attr(feature = "python", pymethods)] +impl SPKSummaryRecord { + pub fn data_type(&self) -> Result { + DafDataType::try_from(self.data_type_i).map_err(|source| EphemerisError::SPK { + action: "converting data type from i32", + source, + }) + } /// Returns the target frame UID of this summary pub fn target_frame(&self) -> Frame { @@ -65,9 +73,16 @@ impl SPKSummaryRecord { Frame::from(self.center_frame_uid()) } - #[cfg(feature = "spkezr_validation")] - pub fn spice_name(&self) -> Result<&'static str, EphemerisError> { - Self::id_to_spice_name(self.target_id) + /// Returns the start epoch of this SPK Summary + #[cfg(feature = "python")] + pub fn start_epoch(&self) -> Epoch { + ::start_epoch(self) + } + + /// Returns the start epoch of this SPK Summary + #[cfg(feature = "python")] + pub fn end_epoch(&self) -> Epoch { + ::end_epoch(self) } /// Converts the provided ID to its human name. diff --git a/anise/src/orientations/rotations.rs b/anise/src/orientations/rotations.rs index 2e2e915a..75a5bebd 100644 --- a/anise/src/orientations/rotations.rs +++ b/anise/src/orientations/rotations.rs @@ -28,7 +28,7 @@ impl Almanac { /// This function only performs the rotation and no translation whatsoever. Use the `transform_from_to` function instead to include rotations. /// /// # Note - /// This function performs a recursion of no more than twice the [MAX_TREE_DEPTH]. + /// This function performs a recursion of no more than twice the MAX_TREE_DEPTH. pub fn rotate_from_to( &self, from_frame: Frame, diff --git a/anise/src/structure/planetocentric/mod.rs b/anise/src/structure/planetocentric/mod.rs index 806fe9d4..58091e4c 100644 --- a/anise/src/structure/planetocentric/mod.rs +++ b/anise/src/structure/planetocentric/mod.rs @@ -234,7 +234,7 @@ impl PlanetaryData { /// Computes the rotation to the parent frame, including its time derivative. /// - /// Source: https://naif.jpl.nasa.gov/pub/naif/toolkit_docs/C/req/rotation.html#Working%20with%20RA,%20Dec%20and%20Twist + /// Source: pub fn rotation_to_parent(&self, epoch: Epoch, system: &Self) -> PhysicsResult { if self.pole_declination.is_none() && self.prime_meridian.is_none() diff --git a/anise/tests/ephemerides/parent_translation_verif.rs b/anise/tests/ephemerides/parent_translation_verif.rs index 7e688766..2e0f4eaa 100644 --- a/anise/tests/ephemerides/parent_translation_verif.rs +++ b/anise/tests/ephemerides/parent_translation_verif.rs @@ -15,13 +15,25 @@ use anise::file2heap; use anise::math::Vector3; use anise::prelude::*; -const ZEROS: &[u8] = &[0; 2048]; -/// Test that we can load data from a static pointer to it. +const ZEROS: &[u8] = &[0; 256]; +/// Test that we can load data from a static pointer to it, even if there is less than one record length #[test] fn invalid_load_from_static() { assert!(SPK::from_static(&ZEROS).is_err()); } +#[test] +fn de400_domain() { + let almanac = Almanac::new("../data/de440s.bsp").unwrap(); + + assert!(almanac.spk_domain(-1012).is_err()); + assert!(almanac.spk_domain(399).is_ok()); + + // No BPC loaded, so it should error. + assert!(almanac.bpc_domain(-1).is_err()); + assert!(almanac.bpc_domain(399).is_err()); +} + #[test] fn de440s_parent_translation_verif() { let _ = pretty_env_logger::try_init(); From 5591b2cd1a5632537d5d418a9591677ab63bdb44 Mon Sep 17 00:00:00 2001 From: Christopher Rabotin Date: Sun, 28 Jan 2024 21:30:47 -0700 Subject: [PATCH 4/4] Update tutorials for SPE --- ...ing Azimuth Elevation and Range data.ipynb | 2 +- ...e kernels and text planetary kernels.ipynb | 60 +- .../Tutorial 06 - Sun probe Earth angle.ipynb | 620 ++++++++++++++++++ anise/src/almanac/bpc.rs | 1 + anise/src/almanac/solar.rs | 2 +- anise/src/almanac/spk.rs | 1 + 6 files changed, 682 insertions(+), 4 deletions(-) create mode 100644 anise-py/tutorials/Tutorial 06 - Sun probe Earth angle.ipynb diff --git a/anise-py/tutorials/Tutorial 04 - Computing Azimuth Elevation and Range data.ipynb b/anise-py/tutorials/Tutorial 04 - Computing Azimuth Elevation and Range data.ipynb index c18b1722..dfc0e95a 100644 --- a/anise-py/tutorials/Tutorial 04 - Computing Azimuth Elevation and Range data.ipynb +++ b/anise-py/tutorials/Tutorial 04 - Computing Azimuth Elevation and Range data.ipynb @@ -22,7 +22,7 @@ "id": "645bb375-58da-42a2-8107-2af6bbb38b5f", "metadata": {}, "source": [ - "## Loading the latest orientation and planteray data\n", + "## Loading the latest orientation and planetary data\n", "\n", "In this tutorial, we're using the latest Earth orientation parameters from JPL, and the latest planetary ephemerides, constants, and high precision Moon rotation parameters. These can be downloaded automatically by ANISE using the MetaAlmanac class (refer to tutorial #02 for details)." ] diff --git a/anise-py/tutorials/Tutorial 05 - Using frame kernels and text planetary kernels.ipynb b/anise-py/tutorials/Tutorial 05 - Using frame kernels and text planetary kernels.ipynb index e0ecb5f0..83dc439a 100644 --- a/anise-py/tutorials/Tutorial 05 - Using frame kernels and text planetary kernels.ipynb +++ b/anise-py/tutorials/Tutorial 05 - Using frame kernels and text planetary kernels.ipynb @@ -246,14 +246,70 @@ "source": [ "## Euler Parameter ANISE kernel\n", "\n", - "Euler parameters are a normalized quaternion. In fact, in the ANISE code in Rust, `EulerParameter` is an alias for `Quaternion`. Euler parameters are useful for defining fixed rotations, e.g. the rotation from the Moon Principal Axes frame to the Moon Mean Earth frame." + "Euler parameters are a normalized quaternion. In fact, in the ANISE code in Rust, `EulerParameter` is an alias for `Quaternion`. Euler parameters are useful for defining fixed rotations, e.g. the rotation from the Moon Principal Axes frame to the Moon Mean Earth frame.\n", + "\n", + "Euler Parameter ANISE kernels (EPA) can be used to store these fixed rotations between frames, and reference them either by their ID or by their name. **They are the equivalent of SPICE's `FK` text files.**\n", + "\n", + "Until [#175](https://github.com/nyx-space/anise/issues/175), rotation data is _not_ exposed to Python. \n", + "\n", + "However, it's possible to convert an FK file into the EPA file using the following function." ] }, { "cell_type": "code", - "execution_count": null, + "execution_count": 1, "id": "263507b1-8f58-44e5-b0b3-4abcf0a00b92", "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\u001b[0;31mSignature:\u001b[0m\n", + "\u001b[0mconvert_fk\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0;34m\u001b[0m\n", + "\u001b[0;34m\u001b[0m \u001b[0mfk_file_path\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", + "\u001b[0;34m\u001b[0m \u001b[0manise_output_path\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", + "\u001b[0;34m\u001b[0m \u001b[0mshow_comments\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", + "\u001b[0;34m\u001b[0m \u001b[0moverwrite\u001b[0m\u001b[0;34m=\u001b[0m\u001b[0;32mNone\u001b[0m\u001b[0;34m,\u001b[0m\u001b[0;34m\u001b[0m\n", + "\u001b[0;34m\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mDocstring:\u001b[0m\n", + "Converts a KPL/FK file, that defines frame constants like fixed rotations, and frame name to ID mappings into the EulerParameterDataSet equivalent ANISE file.\n", + "KPL/FK files must be converted into \"PCA\" (Planetary Constant ANISE) files before being loaded into ANISE.\n", + "\u001b[0;31mType:\u001b[0m builtin_function_or_method" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "from anise.utils import convert_fk\n", + "convert_fk?" + ] + }, + { + "cell_type": "markdown", + "id": "e26f92f9-9fa6-4fb7-b142-8bcaf0451b28", + "metadata": {}, + "source": [ + "On the Nyx Space Cloud, you'll find the `moon.fk` file, which includes the mapping between the Moon PA and Moon ME frame. The latest Almanac also includes the high precision Moon ME frame. Hence, with default data, you can rotate an object from any frame into the high precision Moon ME frame.\n", + "\n", + "## Exercise\n", + "\n", + "Using tutorial 04, compute the azimuth, elevation, and range of the [Shackleton](https://en.wikipedia.org/wiki/Shackleton_(crater)) crater on the Moon to the city of Paris, France.\n", + "\n", + "Here are the general steps:\n", + "\n", + "1. Load the latest Almanac, and check (by printing it) that it includes both EPA and PCA data. Else, load the moon_fk.epa file from the Nyx Space Cloud using a MetaFile with the URL .\n", + "2. Define a time series over a year with a granularity of 12 hours. This crater is on the South Pole of the Moon, and its visibility is often below the horizon of an object as far north as Paris.\n", + "3. For each epoch, define Paris as an `Orbit` instance from its longitude and latitde (recall that the constants include the mean Earth angular rotation rate), in the IAU_EARTH frame. Also build the crater in the IAU_MOON frame.\n", + "4. Finally, call the AER function of the Almanac with each epoch to compute the AER data. Plot it!" + ] + }, + { + "cell_type": "code", + "execution_count": null, + "id": "0a91851b-188a-4dc8-9e8c-40276012dbd1", + "metadata": {}, "outputs": [], "source": [] } diff --git a/anise-py/tutorials/Tutorial 06 - Sun probe Earth angle.ipynb b/anise-py/tutorials/Tutorial 06 - Sun probe Earth angle.ipynb new file mode 100644 index 00000000..0ec4f33f --- /dev/null +++ b/anise-py/tutorials/Tutorial 06 - Sun probe Earth angle.ipynb @@ -0,0 +1,620 @@ +{ + "cells": [ + { + "cell_type": "markdown", + "id": "142de0bb-6053-4122-9777-ea9db5dd61ce", + "metadata": {}, + "source": [ + "# ANISE\n", + "\n", + "ANISE is a modern rewrite of NAIF SPICE, written in Rust and providing interfaces to other languages including Python.\n", + "\n", + "Evidently, this tutorial applies to the Python usage of ANISE.\n", + "\n", + "## Goal\n", + "By the end of this tutorial, you should know how to build a data frame containing the Sun probe Earth angle of a given spacecraft BSP and plot that information. Your exercise will be to confirm that this calculation is correct by computing the Sun elevation at nadir below the spacecraft as detailed in tutorial 04.## Loading the latest orientation and planetary data\n", + "\n", + "Let's start by installing ANISE: `pip install anise`" + ] + }, + { + "cell_type": "markdown", + "id": "20351a7f-e6b0-4d0c-ac6b-0cfabb6d76ca", + "metadata": {}, + "source": [ + "## Load a BSP containing a spacecraft trajectory\n", + "\n", + "In this tutorial, we will use the `gmat-hermite.bsp` file which contains some trajectory data build in GMAT and used in ANISE for validation purposes. Although the Sun probe Earth angle _typically_ applies to spacecraft trajectories, the calculation works for any two objects." + ] + }, + { + "cell_type": "code", + "execution_count": 2, + "id": "59c15415-c54d-465f-962b-585e870d296b", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "=== SPK #0 ===\n", + "┌─────────────┬──────────────────────┬─────────────┬───────────────────────────────────┬───────────────────────────────────┬────────────┬──────────────────────┐\n", + "│ Name │ Target │ Center │ Start epoch │ End epoch │ Duration │ Interpolation kind │\n", + "├─────────────┼──────────────────────┼─────────────┼───────────────────────────────────┼───────────────────────────────────┼────────────┼──────────────────────┤\n", + "│ SPK_SEGMENT │ body -10000001 J2000 │ Earth J2000 │ 2000-01-01T12:00:32.183927328 TDB │ 2000-01-01T15:20:32.183931556 TDB │ 3 h 20 min │ Hermite Unequal Step │\n", + "└─────────────┴──────────────────────┴─────────────┴───────────────────────────────────┴───────────────────────────────────┴────────────┴──────────────────────┘\n", + "=== SPK #1 ===\n", + "┌────────────────┬─────────────────────────────┬─────────��─────────────────────┬───────────────────────────────────┬───────────────────────────────────┬─────────────┬────────────────────┐\n", + "│ Name │ Target │ Center │ Start epoch │ End epoch │ Duration │ Interpolation kind │\n", + "├────────────────┼─────────────────────────────┼───────────────────────────────┼───────────────────────────────────┼───────────────────��───────────────┼─────────────┼────────────────────┤\n", + "│ DE-0440LE-0440 │ Mercury J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Venus J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Earth-Moon Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Mars Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Jupiter Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Saturn Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Uranus Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Neptune Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Pluto Barycenter J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Sun J2000 │ Solar System Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Moon J2000 │ Earth-Moon Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ Earth J2000 │ Earth-Moon Barycenter J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ body 199 J2000 │ Mercury J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "│ DE-0440LE-0440 │ body 299 J2000 │ Venus J2000 │ 1849-12-26T00:00:00.000046818 TDB │ 2150-01-21T23:59:59.999955683 TDB │ 109600 days │ Chebyshev Triplet │\n", + "└────────────────┴─────────────────────────────┴───────────────────────────────┴───────────────────────────────────┴───────────────────────────────────┴─────────────┴────────────────────┘\n" + ] + } + ], + "source": [ + "from anise import MetaAlmanac\n", + "almanac = MetaAlmanac.latest().load(\"../../data/gmat-hermite.bsp\")\n", + "almanac.describe(spk=True)" + ] + }, + { + "cell_type": "markdown", + "id": "9f1d744b-461a-4a38-9d24-b5e989747355", + "metadata": {}, + "source": [ + "We have successfully loaded two BSP files, one with the planetary data and the other with spacecraft data. We now know that this spacecraft has the ID `-10000001`. We can actually query the Almanac for the precise start and stop epochs (returned as hifitime `Epoch` objects) of this ID in our loaded files." + ] + }, + { + "cell_type": "code", + "execution_count": 3, + "id": "c85fc5ec-7b54-4c00-9a1a-9722d7a80cd1", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "2000-01-01T11:59:28.000000021 UTC 2000-01-01T15:19:28.000000226 UTC\n" + ] + } + ], + "source": [ + "start_epoch, stop_epoch = almanac.spk_domain(-10000001)\n", + "print(start_epoch, stop_epoch)" + ] + }, + { + "cell_type": "markdown", + "id": "414d228b-1b69-4352-b415-0ebcd71a669c", + "metadata": {}, + "source": [ + "The Sun Probe Earth angle is the angle between a probe and the Sun and the probe the Earth. It allows one to know whether the point exactly nadir (i.e. below) the spacecraft is illuminated by the Sun or not. This is helpful information if the spacecraft carries a visible light camera and needs to take pictures when it's daytime below the spacecraft.\n", + "\n", + "Let's look at the signature of this function." + ] + }, + { + "cell_type": "code", + "execution_count": 7, + "id": "40176b6b-dfb8-4d9d-a6c2-93cc5c8cd751", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "\u001b[0;31mSignature:\u001b[0m \u001b[0malmanac\u001b[0m\u001b[0;34m.\u001b[0m\u001b[0msun_angle_deg\u001b[0m\u001b[0;34m(\u001b[0m\u001b[0mtarget_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mobserver_id\u001b[0m\u001b[0;34m,\u001b[0m \u001b[0mepoch\u001b[0m\u001b[0;34m)\u001b[0m\u001b[0;34m\u001b[0m\u001b[0;34m\u001b[0m\u001b[0m\n", + "\u001b[0;31mDocstring:\u001b[0m\n", + "Returns the angle (between 0 and 180 degrees) between the observer and the Sun, and the observer and the target body ID.\n", + "This computes the Sun Probe Earth angle (SPE) if the probe is in a loaded, its ID is the \"observer_id\", and the target is set to its central body.\n", + "\n", + "# Geometry\n", + "If the SPE is greater than 90 degrees, then the celestial object below the probe is in sunlight.\n", + "\n", + "## Sunrise at nadir\n", + "```text\n", + "Sun\n", + " | \\ \n", + " | \\\n", + " | \\\n", + " Obs. -- Target\n", + "```\n", + "## Sun high at nadir\n", + "```text\n", + "Sun\n", + " \\ \n", + " \\ __ θ > 90\n", + " \\ \\\n", + " Obs. ---------- Target\n", + "```\n", + "\n", + "## Sunset at nadir\n", + "```text\n", + " Sun\n", + " / \n", + " / __ θ < 90\n", + " / /\n", + " Obs. -- Target\n", + "```\n", + "\n", + "# Algorithm\n", + "1. Compute the position of the Sun as seen from the observer\n", + "2. Compute the position of the target as seen from the observer\n", + "3. Return the arccosine of the dot product of the norms of these vectors.\n", + "\u001b[0;31mType:\u001b[0m builtin_function_or_method" + ] + }, + "metadata": {}, + "output_type": "display_data" + } + ], + "source": [ + "almanac.sun_angle_deg?" + ] + }, + { + "cell_type": "markdown", + "id": "b36bb3ec-8f8a-4298-90ea-2579c3f74ea5", + "metadata": {}, + "source": [ + "One will note that this function is generic to what the \"probe\" SPK ID is (\"observer\") and what its central object should be instead of Earth (\"target\").\n", + "\n", + "The other crucial point here is that this is one of the few functions where the object _ID_ is required instead of a frame. This is because the Almanac will compute everything in the J2000 frame. Don't worry, if you have frame objects instead, you may use the `sun_angle_deg_from_frame` function instead.\n", + "\n", + "Let's see what is the SPE of our spacecraft at the start of the trajectory.\n" + ] + }, + { + "cell_type": "code", + "execution_count": 6, + "id": "0366c5b9-cf4d-473e-9385-3a425fd554bf", + "metadata": {}, + "outputs": [ + { + "data": { + "text/plain": [ + "83.87312777296376" + ] + }, + "execution_count": 6, + "metadata": {}, + "output_type": "execute_result" + } + ], + "source": [ + "from anise.astro.constants import CelestialObjects\n", + "almanac.sun_angle_deg(-10000001, CelestialObjects.EARTH, start_epoch)" + ] + }, + { + "cell_type": "markdown", + "id": "c18f00cf-fdc7-4050-b8b4-4f951e6e11e6", + "metadata": {}, + "source": [ + "A angle of less than 90 degrees means that the nadir point is in the darkness. Let's look at the evolution of the SPE over the duration of the trajectory." + ] + }, + { + "cell_type": "markdown", + "id": "8d8a8d44-b433-4955-a0e5-81988decafce", + "metadata": { + "jp-MarkdownHeadingCollapsed": true + }, + "source": [ + "## Package installation for plotting" + ] + }, + { + "cell_type": "code", + "execution_count": 12, + "id": "b91b712a-3fe7-41a9-9e7d-78cabb756b9d", + "metadata": {}, + "outputs": [ + { + "name": "stdout", + "output_type": "stream", + "text": [ + "Requirement already satisfied: polars[plot] in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (0.20.4)\n", + "Requirement already satisfied: hvplot in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (0.9.1)\n", + "Collecting geoviews\n", + " Downloading geoviews-1.11.0-py2.py3-none-any.whl (511 kB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m511.2/511.2 kB\u001b[0m \u001b[31m690.4 kB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m kB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m:01\u001b[0m\n", + "\u001b[?25hRequirement already satisfied: bokeh>=1.0.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (3.3.3)\n", + "Requirement already satisfied: colorcet>=2 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (3.0.1)\n", + "Requirement already satisfied: holoviews>=1.11.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (1.18.1)\n", + "Requirement already satisfied: pandas in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (2.1.4)\n", + "Requirement already satisfied: numpy>=1.15 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (1.26.3)\n", + "Requirement already satisfied: packaging in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (23.2)\n", + "Requirement already satisfied: panel>=0.11.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (1.3.6)\n", + "Requirement already satisfied: param<3.0,>=1.12.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from hvplot) (2.0.1)\n", + "Collecting cartopy>=0.18.0 (from geoviews)\n", + " Downloading Cartopy-0.22.0-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (11.9 MB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m11.9/11.9 MB\u001b[0m \u001b[31m14.1 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0mm eta \u001b[36m0:00:01\u001b[0m[36m0:00:01\u001b[0m\n", + "\u001b[?25hCollecting shapely (from geoviews)\n", + " Downloading shapely-2.0.2-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (2.5 MB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m2.5/2.5 MB\u001b[0m \u001b[31m54.5 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m31m72.7 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\n", + "\u001b[?25hCollecting pyproj (from geoviews)\n", + " Downloading pyproj-3.6.1-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (8.6 MB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m8.6/8.6 MB\u001b[0m \u001b[31m13.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0mm eta \u001b[36m0:00:01\u001b[0m0:01\u001b[0m:01\u001b[0m\n", + "\u001b[?25hRequirement already satisfied: xyzservices in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from geoviews) (2023.10.1)\n", + "Requirement already satisfied: Jinja2>=2.9 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from bokeh>=1.0.0->hvplot) (3.1.2)\n", + "Requirement already satisfied: contourpy>=1 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from bokeh>=1.0.0->hvplot) (1.2.0)\n", + "Requirement already satisfied: pillow>=7.1.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from bokeh>=1.0.0->hvplot) (10.2.0)\n", + "Requirement already satisfied: PyYAML>=3.10 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from bokeh>=1.0.0->hvplot) (6.0.1)\n", + "Requirement already satisfied: tornado>=5.1 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from bokeh>=1.0.0->hvplot) (6.4)\n", + "Collecting matplotlib>=3.4 (from cartopy>=0.18.0->geoviews)\n", + " Downloading matplotlib-3.8.2-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (11.6 MB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m11.6/11.6 MB\u001b[0m \u001b[31m47.7 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0mm eta \u001b[36m0:00:01\u001b[0m[36m0:00:01\u001b[0m\n", + "\u001b[?25hCollecting pyshp>=2.1 (from cartopy>=0.18.0->geoviews)\n", + " Downloading pyshp-2.3.1-py2.py3-none-any.whl (46 kB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m46.5/46.5 kB\u001b[0m \u001b[31m14.8 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", + "\u001b[?25hRequirement already satisfied: pyct>=0.4.4 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from colorcet>=2->hvplot) (0.5.0)\n", + "Requirement already satisfied: pyviz-comms>=0.7.4 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from holoviews>=1.11.0->hvplot) (3.0.0)\n", + "Requirement already satisfied: python-dateutil>=2.8.2 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from pandas->hvplot) (2.8.2)\n", + "Requirement already satisfied: pytz>=2020.1 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from pandas->hvplot) (2023.3.post1)\n", + "Requirement already satisfied: tzdata>=2022.1 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from pandas->hvplot) (2023.4)\n", + "Requirement already satisfied: markdown in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (3.5.2)\n", + "Requirement already satisfied: markdown-it-py in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (3.0.0)\n", + "Requirement already satisfied: linkify-it-py in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (2.0.2)\n", + "Requirement already satisfied: mdit-py-plugins in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (0.4.0)\n", + "Requirement already satisfied: requests in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (2.31.0)\n", + "Requirement already satisfied: tqdm>=4.48.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (4.66.1)\n", + "Requirement already satisfied: bleach in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (6.1.0)\n", + "Requirement already satisfied: typing-extensions in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from panel>=0.11.0->hvplot) (4.9.0)\n", + "Requirement already satisfied: certifi in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from pyproj->geoviews) (2023.11.17)\n", + "Requirement already satisfied: MarkupSafe>=2.0 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from Jinja2>=2.9->bokeh>=1.0.0->hvplot) (2.1.3)\n", + "Collecting cycler>=0.10 (from matplotlib>=3.4->cartopy>=0.18.0->geoviews)\n", + " Downloading cycler-0.12.1-py3-none-any.whl (8.3 kB)\n", + "Collecting fonttools>=4.22.0 (from matplotlib>=3.4->cartopy>=0.18.0->geoviews)\n", + " Downloading fonttools-4.47.2-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (4.9 MB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m4.9/4.9 MB\u001b[0m \u001b[31m49.5 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m31m56.4 MB/s\u001b[0m eta \u001b[36m0:00:01\u001b[0m\n", + "\u001b[?25hCollecting kiwisolver>=1.3.1 (from matplotlib>=3.4->cartopy>=0.18.0->geoviews)\n", + " Downloading kiwisolver-1.4.5-cp311-cp311-manylinux_2_17_x86_64.manylinux2014_x86_64.whl (1.4 MB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m1.4/1.4 MB\u001b[0m \u001b[31m47.4 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", + "\u001b[?25hCollecting pyparsing>=2.3.1 (from matplotlib>=3.4->cartopy>=0.18.0->geoviews)\n", + " Downloading pyparsing-3.1.1-py3-none-any.whl (103 kB)\n", + "\u001b[2K \u001b[38;2;114;156;31m━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━━\u001b[0m \u001b[32m103.1/103.1 kB\u001b[0m \u001b[31m42.5 MB/s\u001b[0m eta \u001b[36m0:00:00\u001b[0m\n", + "\u001b[?25hRequirement already satisfied: six>=1.5 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from python-dateutil>=2.8.2->pandas->hvplot) (1.16.0)\n", + "Requirement already satisfied: webencodings in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from bleach->panel>=0.11.0->hvplot) (0.5.1)\n", + "Requirement already satisfied: uc-micro-py in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from linkify-it-py->panel>=0.11.0->hvplot) (1.0.2)\n", + "Requirement already satisfied: mdurl~=0.1 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from markdown-it-py->panel>=0.11.0->hvplot) (0.1.2)\n", + "Requirement already satisfied: charset-normalizer<4,>=2 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from requests->panel>=0.11.0->hvplot) (3.3.2)\n", + "Requirement already satisfied: idna<4,>=2.5 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from requests->panel>=0.11.0->hvplot) (3.6)\n", + "Requirement already satisfied: urllib3<3,>=1.21.1 in /home/chris/Workspace/nyx-space/anise/anise-py/.venv/lib64/python3.11/site-packages (from requests->panel>=0.11.0->hvplot) (2.1.0)\n", + "Installing collected packages: shapely, pyshp, pyproj, pyparsing, kiwisolver, fonttools, cycler, matplotlib, cartopy, geoviews\n", + "Successfully installed cartopy-0.22.0 cycler-0.12.1 fonttools-4.47.2 geoviews-1.11.0 kiwisolver-1.4.5 matplotlib-3.8.2 pyparsing-3.1.1 pyproj-3.6.1 pyshp-2.3.1 shapely-2.0.2\n", + "\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m A new release of pip is available: \u001b[0m\u001b[31;49m23.1.2\u001b[0m\u001b[39;49m -> \u001b[0m\u001b[32;49m23.3.2\u001b[0m\n", + "\u001b[1m[\u001b[0m\u001b[34;49mnotice\u001b[0m\u001b[1;39;49m]\u001b[0m\u001b[39;49m To update, run: \u001b[0m\u001b[32;49mpip install --upgrade pip\u001b[0m\n", + "Note: you may need to restart the kernel to use updated packages.\n" + ] + } + ], + "source": [ + "%pip install \"polars[plot]\" hvplot geoviews # geoviews for geographic data" + ] + }, + { + "cell_type": "markdown", + "id": "70a3ddf6-18d7-4f94-94c4-14d73018b2ea", + "metadata": {}, + "source": [ + "## Evolution of SPE over time\n", + "\n", + "We'll plot the change in Sun probe Earth angle over time. We'll also grab the latitude and longitude data so we can plot the position of the spacecraft above the Earth at those times (as a separate plot)." + ] + }, + { + "cell_type": "code", + "execution_count": 38, + "id": "55ec0bf7-2ec8-480f-89be-cd698d0d1e55", + "metadata": {}, + "outputs": [], + "source": [ + "from anise.time import TimeSeries, Unit, Epoch\n", + "from anise.astro import Frame\n", + "from anise.astro.constants import Frames, Orientations\n", + "\n", + "import polars as pl\n", + "from datetime import datetime\n", + "\n", + "def hifitime_to_datetime(e: Epoch) -> datetime:\n", + " return datetime.fromisoformat(str(e).replace(\" UTC\", \"\")[:23])\n", + "\n", + "epochs = []\n", + "spe_deg = []\n", + "lat_deg = []\n", + "long_deg = []\n", + "\n", + "# Let's be sure to load the PCK data, which includes the frame information such as the shape of\n", + "# the ellipsoid and how to compute the rotation of the body fixed frames\n", + "almanac = almanac.load(\"../../data/pck08.pca\")\n", + "\n", + "SC_ID = -10000001\n", + "SC_J2K = Frame(SC_ID, Orientations.J2000)\n", + "\n", + "for epoch in TimeSeries(start_epoch, stop_epoch, Unit.Minute*1, inclusive=True):\n", + " epochs += [hifitime_to_datetime(epoch)]\n", + " spe_deg += [almanac.sun_angle_deg(CelestialObjects.EARTH, SC_ID, epoch)]\n", + " # Grab position of the spacecraft in the IAU Earth frame\n", + " sc_iau_earth = almanac.transform(SC_J2K, Frames.IAU_EARTH_FRAME, epoch)\n", + " lat_deg += [sc_iau_earth.latitude_deg()]\n", + " long_deg += [sc_iau_earth.longitude_deg()]" + ] + }, + { + "cell_type": "code", + "execution_count": 39, + "id": "4f6a34a7-ddae-430d-a138-e3111e7b696c", + "metadata": {}, + "outputs": [], + "source": [ + "# Build the data frame\n", + "df = pl.DataFrame(\n", + " {\n", + " 'Epoch': epochs,\n", + " 'Sun Probe Earth angle (deg)': spe_deg,\n", + " 'Latitude (deg)': lat_deg,\n", + " 'Longitude (deg)': long_deg\n", + " }\n", + ")" + ] + }, + { + "cell_type": "code", + "execution_count": 43, + "id": "614ff2ee-593c-41a6-8932-72671c430ac7", + "metadata": {}, + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Curve [Epoch] (Sun Probe Earth angle (deg),Latitude (deg),Longitude (deg))" + ] + }, + "execution_count": 43, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p2696" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "import hvplot.polars\n", + "from bokeh.models.formatters import DatetimeTickFormatter\n", + "\n", + "formatter = DatetimeTickFormatter(days='%d/%m') # Make the world a better place\n", + "\n", + "df.hvplot(x=\"Epoch\", y=\"Sun Probe Earth angle (deg)\", xformatter=formatter, \n", + " title=\"Sun probe Earth angle over time\", hover_cols=[\"Latitude (deg)\", \"Longitude (deg)\"])" + ] + }, + { + "cell_type": "code", + "execution_count": 41, + "id": "abb43538-5c41-4481-96dc-eefcb36764a2", + "metadata": {}, + "outputs": [ + { + "data": {}, + "metadata": {}, + "output_type": "display_data" + }, + { + "data": { + "application/vnd.holoviews_exec.v0+json": "", + "text/html": [ + "
\n", + "
\n", + "
\n", + "" + ], + "text/plain": [ + ":Overlay\n", + " .WMTS.I :WMTS [Longitude,Latitude]\n", + " .Points.I :Points [Longitude (deg),Latitude (deg)] (Epoch,Sun Probe Earth angle (deg))" + ] + }, + "execution_count": 41, + "metadata": { + "application/vnd.holoviews_exec.v0+json": { + "id": "p2520" + } + }, + "output_type": "execute_result" + } + ], + "source": [ + "df.hvplot.points('Longitude (deg)', 'Latitude (deg)', geo=True, color='red', tiles='ESRI', xlim=(0, 360), ylim=(-60, 60), hover_cols=[\"Epoch\", \"Sun Probe Earth angle (deg)\"])" + ] + }, + { + "cell_type": "markdown", + "id": "4e64369b-ae0e-4977-aec7-a0722e9a06a5", + "metadata": {}, + "source": [ + "## Exercise\n", + "\n", + "The goal of this exercise is to show that the SPE is essentially the Sun elevation angle from the position exactly nadir of the vehicle plus 90 degrees (or so, since the Earth isn't a perfect sphere). It also shows you how to combine the different functionality you've seen with the other tutorials to solve similar problem.\n", + "\n", + "_Note:_ this is the `verify_geometry` Rust test, but rebuilt in Python.\n", + "\n", + "1. For the whole spacecraft trajectory, build the locality exactly nadir of it at each epoch (remember that a `TimeSeries` instance can only be used once, so you'll need to rebuild a new one). For this step, initialize a new `Orbit` instance from the latitude and longitude constructor using the position of the spacecraft in the IAU Earth frame as we did above.\n", + "2. At each epoch, grab the state of the Sun as seen from the Earth (both can be in the J2000 frame since the AER computation will transform the states into the correct frames anyway).\n", + "3. Call the azimuth, elevtion, and range function of your loaded Almanac, and store the elevation of the Sun as seen from the point exactly nadir of the spacecraft.\n", + "4. Plot the elevation data in degrees compared to the SPE angle calculated above." + ] + } + ], + "metadata": { + "kernelspec": { + "display_name": ".venv", + "language": "python", + "name": ".venv" + }, + "language_info": { + "codemirror_mode": { + "name": "ipython", + "version": 3 + }, + "file_extension": ".py", + "mimetype": "text/x-python", + "name": "python", + "nbconvert_exporter": "python", + "pygments_lexer": "ipython3", + "version": "3.11.4" + } + }, + "nbformat": 4, + "nbformat_minor": 5 +} diff --git a/anise/src/almanac/bpc.rs b/anise/src/almanac/bpc.rs index bee4ab7c..63d9c3f0 100644 --- a/anise/src/almanac/bpc.rs +++ b/anise/src/almanac/bpc.rs @@ -13,6 +13,7 @@ use hifitime::Epoch; #[cfg(feature = "python")] use pyo3::prelude::*; +#[cfg(not(feature = "python"))] use crate::naif::daf::NAIFSummaryRecord; use crate::naif::pck::BPCSummaryRecord; use crate::naif::BPC; diff --git a/anise/src/almanac/solar.rs b/anise/src/almanac/solar.rs index 6898d7db..6bf0e03f 100644 --- a/anise/src/almanac/solar.rs +++ b/anise/src/almanac/solar.rs @@ -20,7 +20,7 @@ use pyo3::prelude::*; #[cfg_attr(feature = "python", pymethods)] impl Almanac { /// Returns the angle (between 0 and 180 degrees) between the observer and the Sun, and the observer and the target body ID. - /// This computes the Sun Probe Earth angle (SPE) if the probe is in a loaded, its ID is the "observer_id", and the target is set to its central body. + /// This computes the Sun Probe Earth angle (SPE) if the probe is in a loaded SPK, its ID is the "observer_id", and the target is set to its central body. /// /// # Geometry /// If the SPE is greater than 90 degrees, then the celestial object below the probe is in sunlight. diff --git a/anise/src/almanac/spk.rs b/anise/src/almanac/spk.rs index eee13c6f..8954e6de 100644 --- a/anise/src/almanac/spk.rs +++ b/anise/src/almanac/spk.rs @@ -14,6 +14,7 @@ use hifitime::Epoch; use pyo3::prelude::*; use crate::naif::daf::DAFError; +#[cfg(not(feature = "python"))] use crate::naif::daf::NAIFSummaryRecord; use crate::naif::spk::summary::SPKSummaryRecord; use crate::naif::SPK;