pub struct ODProcess<'a, D: Dynamics, MsrSize: DimName, Accel: DimName, K: Filter<D::StateType, Accel, MsrSize>, Trk: TrackerSensitivity<D::StateType, D::StateType>>where
D::StateType: Interpolatable + Add<OVector<f64, <D::StateType as State>::Size>, Output = D::StateType>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <D::StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<D::StateType as State>::Size, Accel> + Allocator<Accel, <D::StateType as State>::Size>,{
pub prop: PropInstance<'a, D>,
pub kf: K,
pub devices: BTreeMap<String, Trk>,
pub estimates: Vec<K::Estimate>,
pub residuals: Vec<Option<Residual<MsrSize>>>,
pub ekf_trigger: Option<EkfTrigger>,
pub resid_crit: Option<ResidRejectCrit>,
pub almanac: Arc<Almanac>,
/* private fields */
}
Expand description
An orbit determination process. Note that everything passed to this structure is moved.
Fields§
§prop: PropInstance<'a, D>
PropInstance used for the estimation
kf: K
Kalman filter itself
devices: BTreeMap<String, Trk>
Tracking devices
estimates: Vec<K::Estimate>
Vector of estimates available after a pass
residuals: Vec<Option<Residual<MsrSize>>>
Vector of residuals available after a pass
ekf_trigger: Option<EkfTrigger>
§resid_crit: Option<ResidRejectCrit>
Residual rejection criteria allows preventing bad measurements from affecting the estimation.
almanac: Arc<Almanac>
Implementations§
Source§impl<MsrSize: DimName, Accel: DimName, Trk: TrackerSensitivity<Spacecraft, Spacecraft>> ODProcess<'_, SpacecraftDynamics, MsrSize, Accel, KF<Spacecraft, Accel, MsrSize>, Trk>where
DefaultAllocator: Allocator<MsrSize> + Allocator<MsrSize, <Spacecraft as State>::Size> + Allocator<Const<1>, MsrSize> + Allocator<<Spacecraft as State>::Size> + Allocator<<Spacecraft as State>::Size, <Spacecraft as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<Spacecraft as State>::Size, MsrSize> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<Spacecraft as State>::VecLength> + Allocator<<Spacecraft as State>::Size, Accel> + Allocator<Accel, <Spacecraft as State>::Size>,
impl<MsrSize: DimName, Accel: DimName, Trk: TrackerSensitivity<Spacecraft, Spacecraft>> ODProcess<'_, SpacecraftDynamics, MsrSize, Accel, KF<Spacecraft, Accel, MsrSize>, Trk>where
DefaultAllocator: Allocator<MsrSize> + Allocator<MsrSize, <Spacecraft as State>::Size> + Allocator<Const<1>, MsrSize> + Allocator<<Spacecraft as State>::Size> + Allocator<<Spacecraft as State>::Size, <Spacecraft as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<Spacecraft as State>::Size, MsrSize> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<Spacecraft as State>::VecLength> + Allocator<<Spacecraft as State>::Size, Accel> + Allocator<Accel, <Spacecraft as State>::Size>,
Sourcepub fn to_parquet<P: AsRef<Path>>(
&self,
arc: &TrackingDataArc,
path: P,
cfg: ExportCfg,
) -> Result<PathBuf, ODError>
pub fn to_parquet<P: AsRef<Path>>( &self, arc: &TrackingDataArc, path: P, cfg: ExportCfg, ) -> Result<PathBuf, ODError>
Store the estimates and residuals in a parquet file
Examples found in repository?
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
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 regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
// Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
let mut latest_jwst_ephem = MetaFile {
uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
crc32: None,
};
latest_jwst_ephem.process(true)?;
// 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(latest_jwst_ephem, true)?,
);
// By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
// in the BSP. We need this ID in order to query the ephemeris.
const JWST_NAIF_ID: i32 = -170;
// Let's build a frame in the J2000 orientation centered on the JWST.
const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
// Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
println!("JWST defined from {earliest_epoch} to {latest_epoch}");
// Fetch the state, printing it in the Earth J2000 frame.
let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
println!("{jwst_orbit:x}");
// Build the spacecraft
// SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
// SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
let jwst = Spacecraft::builder()
.orbit(jwst_orbit)
.srp(SrpConfig {
area_m2: 21.197 * 14.162,
cr: 1.56,
})
.dry_mass_kg(6200.0)
.build();
// Build up the spacecraft uncertainty builder.
// We can use the spacecraft uncertainty structure to build this up.
// We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
// in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
// Nyx can also estimate the deviation of the spacecraft parameters.
let jwst_uncertainty = SpacecraftUncertainty::builder()
.nominal(jwst)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.3)
.z_km(1.5)
.vx_km_s(1e-4)
.vy_km_s(0.6e-3)
.vz_km_s(3e-3)
.build();
println!("{jwst_uncertainty}");
// Build the Kalman filter estimate.
// Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
// but this approach requires quite a bit more boilerplate code.
let jwst_estimate = jwst_uncertainty.to_estimate()?;
// Set up the spacecraft dynamics.
// We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
// We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
// Build the propagator set up to use for the whole analysis.
let setup = Propagator::default(dynamics);
// All of the analysis will use this duration.
let prediction_duration = 6.5 * Unit::Day;
// === Covariance mapping ===
// For the covariance mapping / prediction, we'll use the common orbit determination approach.
// This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
let ckf = KF::no_snc(jwst_estimate);
// Build the propagation instance for the OD process.
let prop = setup.with(jwst.with_stm(), almanac.clone());
let mut odp = SpacecraftODProcess::ckf(prop, ckf, BTreeMap::new(), None, almanac.clone());
// Define the prediction step, i.e. how often we want to know the covariance.
let step = 1_i64.minutes();
// Finally, predict, and export the trajectory with covariance to a parquet file.
odp.predict_for(step, prediction_duration)?;
odp.to_parquet(
&TrackingDataArc::default(),
"./02_jwst_covar_map.parquet",
ExportCfg::default(),
)?;
// === Monte Carlo framework ===
// Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
let my_mc = MonteCarlo::new(
jwst, // Nominal state
jwst_estimate.to_random_variable()?,
"02_jwst".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
let num_runs = 5_000;
let rslts = my_mc.run_until_epoch(
setup,
almanac.clone(),
jwst.epoch() + prediction_duration,
num_runs,
);
assert_eq!(rslts.runs.len(), num_runs);
// Finally, export these results, computing the eclipse percentage for all of these results.
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
let umbra_event = eclipse_loc.to_umbra_event();
let penumbra_event = eclipse_loc.to_penumbra_event();
rslts.to_parquet(
"02_jwst_monte_carlo.parquet",
Some(vec![&umbra_event, &penumbra_event]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
More examples
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// ====================== //
// === ALMANAC SET UP === //
// ====================== //
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's MetaAlmanac.
let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
.iter()
.collect();
let meta = data_folder.join("lro-dynamics.dhall");
// Load this ephem in the general Almanac we're using for this analysis.
let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
.map_err(Box::new)?
.process(true)
.map_err(Box::new)?;
let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
moon_pc.mu_km3_s2 = 4902.74987;
almanac.planetary_data.set_by_id(MOON, moon_pc)?;
let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
earth_pc.mu_km3_s2 = 398600.436;
almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
// Save this new kernel for reuse.
// In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
almanac
.planetary_data
.save_as(&data_folder.join("lro-specific.pca"), true)?;
// Lock the almanac (an Arc is a read only structure).
let almanac = Arc::new(almanac);
// 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.96,
})
.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 traj_as_flown = Traj::from_bsp(
lro_frame,
MOON_J2000,
almanac.clone(),
sc_template,
5.seconds(),
Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
Aberration::LT,
Some("LRO".to_string()),
)?;
println!("{traj_as_flown}");
// ====================== //
// === MODEL MATCHING === //
// ====================== //
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
// The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the GRAIL JGGRX model.
let mut jggrx_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jggrx_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Moon principal axes frame.
let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
// let moon_pa_frame = IAU_MOON_FRAME;
let sph_harmonics = Harmonics::from_stor(
almanac.frame_from_uid(moon_pa_frame)?,
HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(sph_harmonics);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth and Moon.
// Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
println!("{dynamics}");
// Now we can build the propagator.
let setup = Propagator::default_dp78(dynamics.clone());
// For reference, let's build the trajectory with Nyx's models from that LRO state.
let (sim_final, traj_as_sim) = setup
.with(*traj_as_flown.first(), almanac.clone())
.until_epoch_with_traj(traj_as_flown.last().epoch())?;
println!("SIM INIT: {:x}", traj_as_flown.first());
println!("SIM FINAL: {sim_final:x}");
// Compute RIC difference between SIM and LRO ephem
let sim_lro_delta = sim_final
.orbit
.ric_difference(&traj_as_flown.last().orbit)?;
println!("{traj_as_sim}");
println!(
"SIM v LRO - RIC Position (m): {:.3}",
sim_lro_delta.radius_km * 1e3
);
println!(
"SIM v LRO - RIC Velocity (m/s): {:.3}",
sim_lro_delta.velocity_km_s * 1e3
);
traj_as_sim.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_sim_truth_error.parquet",
ExportCfg::default(),
)?;
// ==================== //
// === OD SIMULATOR === //
// ==================== //
// After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
// and the truth LRO state.
// Therefore, we will actually run an estimation from a dispersed LRO state.
// The sc_seed is the true LRO state from the BSP.
let sc_seed = *traj_as_flown.first();
// 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_named(ground_station_file)?;
// 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)?;
// Build the tracking arc simulation to generate a "standard measurement".
let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
devices.clone(),
traj_as_flown.clone(),
configs,
)?;
trk.build_schedule(almanac.clone())?;
let arc = trk.generate_measurements(almanac.clone())?;
// Save the simulated tracking data
arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
// We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
println!("{arc}");
// Now that we have simulated measurements, we'll run the orbit determination.
// ===================== //
// === OD ESTIMATION === //
// ===================== //
let sc = SpacecraftUncertainty::builder()
.nominal(sc_seed)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.5)
.z_km(0.5)
.vx_km_s(5e-3)
.vy_km_s(5e-3)
.vz_km_s(5e-3)
.build();
// Build the filter initial estimate, which we will reuse in the filter.
let initial_estimate = sc.to_estimate()?;
println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
let kf = KF::new(
// Increase the initial covariance to account for larger deviation.
initial_estimate,
// Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame.
SNC3::from_diagonal(10 * Unit::Minute, &[1e-12, 1e-12, 1e-12]),
);
// We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
let mut odp = SpacecraftODProcess::ckf(
setup.with(initial_estimate.state().with_stm(), almanac.clone()),
kf,
devices,
Some(ResidRejectCrit::default()),
almanac.clone(),
);
odp.process_arc(&arc)?;
let ric_err = traj_as_flown
.at(odp.estimates.last().unwrap().epoch())?
.orbit
.ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
println!("== RIC at end ==");
println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;
// In our case, we have the truth trajectory from NASA.
// So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
// Export the OD trajectory first.
let od_trajectory = odp.to_traj()?;
// Build the RIC difference.
od_trajectory.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_od_truth_error.parquet",
ExportCfg::default(),
)?;
Ok(())
}
Source§impl<'a, D: Dynamics, MsrSize: DimName, Accel: DimName, K: Filter<D::StateType, Accel, MsrSize>, Trk: TrackerSensitivity<D::StateType, D::StateType>> ODProcess<'a, D, MsrSize, Accel, K, Trk>where
D::StateType: Interpolatable + Add<OVector<f64, <D::StateType as State>::Size>, Output = D::StateType>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <D::StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<D::StateType as State>::Size, Accel> + Allocator<Accel, <D::StateType as State>::Size>,
impl<'a, D: Dynamics, MsrSize: DimName, Accel: DimName, K: Filter<D::StateType, Accel, MsrSize>, Trk: TrackerSensitivity<D::StateType, D::StateType>> ODProcess<'a, D, MsrSize, Accel, K, Trk>where
D::StateType: Interpolatable + Add<OVector<f64, <D::StateType as State>::Size>, Output = D::StateType>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <D::StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<D::StateType as State>::Size, Accel> + Allocator<Accel, <D::StateType as State>::Size>,
Sourcepub fn new(
prop: PropInstance<'a, D>,
kf: K,
devices: BTreeMap<String, Trk>,
ekf_trigger: Option<EkfTrigger>,
resid_crit: Option<ResidRejectCrit>,
almanac: Arc<Almanac>,
) -> Self
pub fn new( prop: PropInstance<'a, D>, kf: K, devices: BTreeMap<String, Trk>, ekf_trigger: Option<EkfTrigger>, resid_crit: Option<ResidRejectCrit>, almanac: Arc<Almanac>, ) -> Self
Initialize a new orbit determination process with an optional trigger to switch from a CKF to an EKF.
Sourcepub fn ekf(
prop: PropInstance<'a, D>,
kf: K,
devices: BTreeMap<String, Trk>,
trigger: EkfTrigger,
resid_crit: Option<ResidRejectCrit>,
almanac: Arc<Almanac>,
) -> Self
pub fn ekf( prop: PropInstance<'a, D>, kf: K, devices: BTreeMap<String, Trk>, trigger: EkfTrigger, resid_crit: Option<ResidRejectCrit>, almanac: Arc<Almanac>, ) -> Self
Initialize a new orbit determination process with an Extended Kalman filter. The switch from a classical KF to an EKF is based on the provided trigger.
Sourcepub fn smooth(
&self,
condition: SmoothingArc,
) -> Result<Vec<K::Estimate>, ODError>
pub fn smooth( &self, condition: SmoothingArc, ) -> Result<Vec<K::Estimate>, ODError>
Allows to smooth the provided estimates. Returns the smoothed estimates or an error.
Estimates must be ordered in chronological order. This function will smooth the estimates from the last in the list to the first one.
Sourcepub fn rms_residual_ratios(&self) -> f64
pub fn rms_residual_ratios(&self) -> f64
Returns the root mean square of the prefit residual ratios
Sourcepub fn iterate_arc(
&mut self,
arc: &TrackingDataArc,
config: IterationConf,
) -> Result<(), ODError>
pub fn iterate_arc( &mut self, arc: &TrackingDataArc, config: IterationConf, ) -> Result<(), ODError>
Allows iterating on the filter solution. Requires specifying a smoothing condition to know where to stop the smoothing.
Sourcepub fn process_arc(&mut self, arc: &TrackingDataArc) -> Result<(), ODError>
pub fn process_arc(&mut self, arc: &TrackingDataArc) -> Result<(), ODError>
Process the provided measurements for this orbit determination process given the associated devices.
§Argument details
- The measurements must be a list mapping the name of the measurement device to the measurement itself.
- The name of all measurement devices must be present in the provided devices, i.e. the key set of
devices
must be a superset of the measurement device names present in the list. - The maximum step size to ensure we don’t skip any measurements.
Examples found in repository?
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// ====================== //
// === ALMANAC SET UP === //
// ====================== //
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's MetaAlmanac.
let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
.iter()
.collect();
let meta = data_folder.join("lro-dynamics.dhall");
// Load this ephem in the general Almanac we're using for this analysis.
let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
.map_err(Box::new)?
.process(true)
.map_err(Box::new)?;
let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
moon_pc.mu_km3_s2 = 4902.74987;
almanac.planetary_data.set_by_id(MOON, moon_pc)?;
let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
earth_pc.mu_km3_s2 = 398600.436;
almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
// Save this new kernel for reuse.
// In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
almanac
.planetary_data
.save_as(&data_folder.join("lro-specific.pca"), true)?;
// Lock the almanac (an Arc is a read only structure).
let almanac = Arc::new(almanac);
// 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.96,
})
.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 traj_as_flown = Traj::from_bsp(
lro_frame,
MOON_J2000,
almanac.clone(),
sc_template,
5.seconds(),
Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
Aberration::LT,
Some("LRO".to_string()),
)?;
println!("{traj_as_flown}");
// ====================== //
// === MODEL MATCHING === //
// ====================== //
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
// The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the GRAIL JGGRX model.
let mut jggrx_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jggrx_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Moon principal axes frame.
let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
// let moon_pa_frame = IAU_MOON_FRAME;
let sph_harmonics = Harmonics::from_stor(
almanac.frame_from_uid(moon_pa_frame)?,
HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(sph_harmonics);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth and Moon.
// Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
println!("{dynamics}");
// Now we can build the propagator.
let setup = Propagator::default_dp78(dynamics.clone());
// For reference, let's build the trajectory with Nyx's models from that LRO state.
let (sim_final, traj_as_sim) = setup
.with(*traj_as_flown.first(), almanac.clone())
.until_epoch_with_traj(traj_as_flown.last().epoch())?;
println!("SIM INIT: {:x}", traj_as_flown.first());
println!("SIM FINAL: {sim_final:x}");
// Compute RIC difference between SIM and LRO ephem
let sim_lro_delta = sim_final
.orbit
.ric_difference(&traj_as_flown.last().orbit)?;
println!("{traj_as_sim}");
println!(
"SIM v LRO - RIC Position (m): {:.3}",
sim_lro_delta.radius_km * 1e3
);
println!(
"SIM v LRO - RIC Velocity (m/s): {:.3}",
sim_lro_delta.velocity_km_s * 1e3
);
traj_as_sim.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_sim_truth_error.parquet",
ExportCfg::default(),
)?;
// ==================== //
// === OD SIMULATOR === //
// ==================== //
// After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
// and the truth LRO state.
// Therefore, we will actually run an estimation from a dispersed LRO state.
// The sc_seed is the true LRO state from the BSP.
let sc_seed = *traj_as_flown.first();
// 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_named(ground_station_file)?;
// 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)?;
// Build the tracking arc simulation to generate a "standard measurement".
let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
devices.clone(),
traj_as_flown.clone(),
configs,
)?;
trk.build_schedule(almanac.clone())?;
let arc = trk.generate_measurements(almanac.clone())?;
// Save the simulated tracking data
arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
// We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
println!("{arc}");
// Now that we have simulated measurements, we'll run the orbit determination.
// ===================== //
// === OD ESTIMATION === //
// ===================== //
let sc = SpacecraftUncertainty::builder()
.nominal(sc_seed)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.5)
.z_km(0.5)
.vx_km_s(5e-3)
.vy_km_s(5e-3)
.vz_km_s(5e-3)
.build();
// Build the filter initial estimate, which we will reuse in the filter.
let initial_estimate = sc.to_estimate()?;
println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
let kf = KF::new(
// Increase the initial covariance to account for larger deviation.
initial_estimate,
// Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame.
SNC3::from_diagonal(10 * Unit::Minute, &[1e-12, 1e-12, 1e-12]),
);
// We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
let mut odp = SpacecraftODProcess::ckf(
setup.with(initial_estimate.state().with_stm(), almanac.clone()),
kf,
devices,
Some(ResidRejectCrit::default()),
almanac.clone(),
);
odp.process_arc(&arc)?;
let ric_err = traj_as_flown
.at(odp.estimates.last().unwrap().epoch())?
.orbit
.ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
println!("== RIC at end ==");
println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;
// In our case, we have the truth trajectory from NASA.
// So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
// Export the OD trajectory first.
let od_trajectory = odp.to_traj()?;
// Build the RIC difference.
od_trajectory.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_od_truth_error.parquet",
ExportCfg::default(),
)?;
Ok(())
}
Sourcepub fn predict_until(
&mut self,
step: Duration,
end_epoch: Epoch,
) -> Result<(), ODError>
pub fn predict_until( &mut self, step: Duration, end_epoch: Epoch, ) -> Result<(), ODError>
Continuously predicts the trajectory until the provided end epoch, with covariance mapping at each step. In other words, this performs a time update.
Sourcepub fn predict_for(
&mut self,
step: Duration,
duration: Duration,
) -> Result<(), ODError>
pub fn predict_for( &mut self, step: Duration, duration: Duration, ) -> Result<(), ODError>
Continuously predicts the trajectory for the provided duration, with covariance mapping at each step. In other words, this performs a time update.
Examples found in repository?
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
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 regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
// Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
let mut latest_jwst_ephem = MetaFile {
uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
crc32: None,
};
latest_jwst_ephem.process(true)?;
// 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(latest_jwst_ephem, true)?,
);
// By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
// in the BSP. We need this ID in order to query the ephemeris.
const JWST_NAIF_ID: i32 = -170;
// Let's build a frame in the J2000 orientation centered on the JWST.
const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
// Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
println!("JWST defined from {earliest_epoch} to {latest_epoch}");
// Fetch the state, printing it in the Earth J2000 frame.
let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
println!("{jwst_orbit:x}");
// Build the spacecraft
// SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
// SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
let jwst = Spacecraft::builder()
.orbit(jwst_orbit)
.srp(SrpConfig {
area_m2: 21.197 * 14.162,
cr: 1.56,
})
.dry_mass_kg(6200.0)
.build();
// Build up the spacecraft uncertainty builder.
// We can use the spacecraft uncertainty structure to build this up.
// We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
// in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
// Nyx can also estimate the deviation of the spacecraft parameters.
let jwst_uncertainty = SpacecraftUncertainty::builder()
.nominal(jwst)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.3)
.z_km(1.5)
.vx_km_s(1e-4)
.vy_km_s(0.6e-3)
.vz_km_s(3e-3)
.build();
println!("{jwst_uncertainty}");
// Build the Kalman filter estimate.
// Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
// but this approach requires quite a bit more boilerplate code.
let jwst_estimate = jwst_uncertainty.to_estimate()?;
// Set up the spacecraft dynamics.
// We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
// We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
// Build the propagator set up to use for the whole analysis.
let setup = Propagator::default(dynamics);
// All of the analysis will use this duration.
let prediction_duration = 6.5 * Unit::Day;
// === Covariance mapping ===
// For the covariance mapping / prediction, we'll use the common orbit determination approach.
// This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
let ckf = KF::no_snc(jwst_estimate);
// Build the propagation instance for the OD process.
let prop = setup.with(jwst.with_stm(), almanac.clone());
let mut odp = SpacecraftODProcess::ckf(prop, ckf, BTreeMap::new(), None, almanac.clone());
// Define the prediction step, i.e. how often we want to know the covariance.
let step = 1_i64.minutes();
// Finally, predict, and export the trajectory with covariance to a parquet file.
odp.predict_for(step, prediction_duration)?;
odp.to_parquet(
&TrackingDataArc::default(),
"./02_jwst_covar_map.parquet",
ExportCfg::default(),
)?;
// === Monte Carlo framework ===
// Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
let my_mc = MonteCarlo::new(
jwst, // Nominal state
jwst_estimate.to_random_variable()?,
"02_jwst".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
let num_runs = 5_000;
let rslts = my_mc.run_until_epoch(
setup,
almanac.clone(),
jwst.epoch() + prediction_duration,
num_runs,
);
assert_eq!(rslts.runs.len(), num_runs);
// Finally, export these results, computing the eclipse percentage for all of these results.
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
let umbra_event = eclipse_loc.to_umbra_event();
let penumbra_event = eclipse_loc.to_penumbra_event();
rslts.to_parquet(
"02_jwst_monte_carlo.parquet",
Some(vec![&umbra_event, &penumbra_event]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
Sourcepub fn to_traj(&self) -> Result<Traj<D::StateType>, NyxError>
pub fn to_traj(&self) -> Result<Traj<D::StateType>, NyxError>
Builds the navigation trajectory for the estimated state only
Examples found in repository?
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// ====================== //
// === ALMANAC SET UP === //
// ====================== //
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's MetaAlmanac.
let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
.iter()
.collect();
let meta = data_folder.join("lro-dynamics.dhall");
// Load this ephem in the general Almanac we're using for this analysis.
let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
.map_err(Box::new)?
.process(true)
.map_err(Box::new)?;
let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
moon_pc.mu_km3_s2 = 4902.74987;
almanac.planetary_data.set_by_id(MOON, moon_pc)?;
let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
earth_pc.mu_km3_s2 = 398600.436;
almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
// Save this new kernel for reuse.
// In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
almanac
.planetary_data
.save_as(&data_folder.join("lro-specific.pca"), true)?;
// Lock the almanac (an Arc is a read only structure).
let almanac = Arc::new(almanac);
// 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.96,
})
.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 traj_as_flown = Traj::from_bsp(
lro_frame,
MOON_J2000,
almanac.clone(),
sc_template,
5.seconds(),
Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
Aberration::LT,
Some("LRO".to_string()),
)?;
println!("{traj_as_flown}");
// ====================== //
// === MODEL MATCHING === //
// ====================== //
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
// The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the GRAIL JGGRX model.
let mut jggrx_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jggrx_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Moon principal axes frame.
let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
// let moon_pa_frame = IAU_MOON_FRAME;
let sph_harmonics = Harmonics::from_stor(
almanac.frame_from_uid(moon_pa_frame)?,
HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(sph_harmonics);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth and Moon.
// Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
println!("{dynamics}");
// Now we can build the propagator.
let setup = Propagator::default_dp78(dynamics.clone());
// For reference, let's build the trajectory with Nyx's models from that LRO state.
let (sim_final, traj_as_sim) = setup
.with(*traj_as_flown.first(), almanac.clone())
.until_epoch_with_traj(traj_as_flown.last().epoch())?;
println!("SIM INIT: {:x}", traj_as_flown.first());
println!("SIM FINAL: {sim_final:x}");
// Compute RIC difference between SIM and LRO ephem
let sim_lro_delta = sim_final
.orbit
.ric_difference(&traj_as_flown.last().orbit)?;
println!("{traj_as_sim}");
println!(
"SIM v LRO - RIC Position (m): {:.3}",
sim_lro_delta.radius_km * 1e3
);
println!(
"SIM v LRO - RIC Velocity (m/s): {:.3}",
sim_lro_delta.velocity_km_s * 1e3
);
traj_as_sim.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_sim_truth_error.parquet",
ExportCfg::default(),
)?;
// ==================== //
// === OD SIMULATOR === //
// ==================== //
// After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
// and the truth LRO state.
// Therefore, we will actually run an estimation from a dispersed LRO state.
// The sc_seed is the true LRO state from the BSP.
let sc_seed = *traj_as_flown.first();
// 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_named(ground_station_file)?;
// 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)?;
// Build the tracking arc simulation to generate a "standard measurement".
let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
devices.clone(),
traj_as_flown.clone(),
configs,
)?;
trk.build_schedule(almanac.clone())?;
let arc = trk.generate_measurements(almanac.clone())?;
// Save the simulated tracking data
arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
// We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
println!("{arc}");
// Now that we have simulated measurements, we'll run the orbit determination.
// ===================== //
// === OD ESTIMATION === //
// ===================== //
let sc = SpacecraftUncertainty::builder()
.nominal(sc_seed)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.5)
.z_km(0.5)
.vx_km_s(5e-3)
.vy_km_s(5e-3)
.vz_km_s(5e-3)
.build();
// Build the filter initial estimate, which we will reuse in the filter.
let initial_estimate = sc.to_estimate()?;
println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
let kf = KF::new(
// Increase the initial covariance to account for larger deviation.
initial_estimate,
// Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame.
SNC3::from_diagonal(10 * Unit::Minute, &[1e-12, 1e-12, 1e-12]),
);
// We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
let mut odp = SpacecraftODProcess::ckf(
setup.with(initial_estimate.state().with_stm(), almanac.clone()),
kf,
devices,
Some(ResidRejectCrit::default()),
almanac.clone(),
);
odp.process_arc(&arc)?;
let ric_err = traj_as_flown
.at(odp.estimates.last().unwrap().epoch())?
.orbit
.ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
println!("== RIC at end ==");
println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;
// In our case, we have the truth trajectory from NASA.
// So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
// Export the OD trajectory first.
let od_trajectory = odp.to_traj()?;
// Build the RIC difference.
od_trajectory.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_od_truth_error.parquet",
ExportCfg::default(),
)?;
Ok(())
}
Source§impl<'a, D: Dynamics, MsrSize: DimName, Accel: DimName, K: Filter<D::StateType, Accel, MsrSize>, Trk: TrackerSensitivity<D::StateType, D::StateType>> ODProcess<'a, D, MsrSize, Accel, K, Trk>where
D::StateType: Interpolatable + Add<OVector<f64, <D::StateType as State>::Size>, Output = D::StateType>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <D::StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<D::StateType as State>::Size, Accel> + Allocator<Accel, <D::StateType as State>::Size>,
impl<'a, D: Dynamics, MsrSize: DimName, Accel: DimName, K: Filter<D::StateType, Accel, MsrSize>, Trk: TrackerSensitivity<D::StateType, D::StateType>> ODProcess<'a, D, MsrSize, Accel, K, Trk>where
D::StateType: Interpolatable + Add<OVector<f64, <D::StateType as State>::Size>, Output = D::StateType>,
<DefaultAllocator as Allocator<<D::StateType as State>::VecLength>>::Buffer<f64>: Send,
DefaultAllocator: Allocator<<D::StateType as State>::Size> + Allocator<<D::StateType as State>::VecLength> + Allocator<MsrSize> + Allocator<MsrSize, <D::StateType as State>::Size> + Allocator<MsrSize, MsrSize> + Allocator<<D::StateType as State>::Size, <D::StateType as State>::Size> + Allocator<Accel> + Allocator<Accel, Accel> + Allocator<<D::StateType as State>::Size, Accel> + Allocator<Accel, <D::StateType as State>::Size>,
Sourcepub fn ckf(
prop: PropInstance<'a, D>,
kf: K,
devices: BTreeMap<String, Trk>,
resid_crit: Option<ResidRejectCrit>,
almanac: Arc<Almanac>,
) -> Self
pub fn ckf( prop: PropInstance<'a, D>, kf: K, devices: BTreeMap<String, Trk>, resid_crit: Option<ResidRejectCrit>, almanac: Arc<Almanac>, ) -> Self
Examples found in repository?
26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165
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 regularly update of the James Webb Space Telescope reconstucted (or definitive) ephemeris.
// Refer to https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/aareadme.txt for details.
let mut latest_jwst_ephem = MetaFile {
uri: "https://naif.jpl.nasa.gov/pub/naif/JWST/kernels/spk/jwst_rec.bsp".to_string(),
crc32: None,
};
latest_jwst_ephem.process(true)?;
// 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(latest_jwst_ephem, true)?,
);
// By loading this ephemeris file in the ANISE GUI or ANISE CLI, we can find the NAIF ID of the JWST
// in the BSP. We need this ID in order to query the ephemeris.
const JWST_NAIF_ID: i32 = -170;
// Let's build a frame in the J2000 orientation centered on the JWST.
const JWST_J2000: Frame = Frame::from_ephem_j2000(JWST_NAIF_ID);
// Since the ephemeris file is updated regularly, we'll just grab the latest state in the ephem.
let (earliest_epoch, latest_epoch) = almanac.spk_domain(JWST_NAIF_ID)?;
println!("JWST defined from {earliest_epoch} to {latest_epoch}");
// Fetch the state, printing it in the Earth J2000 frame.
let jwst_orbit = almanac.transform(JWST_J2000, EARTH_J2000, latest_epoch, None)?;
println!("{jwst_orbit:x}");
// Build the spacecraft
// SRP area assumed to be the full sunshield and mass if 6200.0 kg, c.f. https://webb.nasa.gov/content/about/faqs/facts.html
// SRP Coefficient of reflectivity assumed to be that of Kapton, i.e. 2 - 0.44 = 1.56, table 1 from https://amostech.com/TechnicalPapers/2018/Poster/Bengtson.pdf
let jwst = Spacecraft::builder()
.orbit(jwst_orbit)
.srp(SrpConfig {
area_m2: 21.197 * 14.162,
cr: 1.56,
})
.dry_mass_kg(6200.0)
.build();
// Build up the spacecraft uncertainty builder.
// We can use the spacecraft uncertainty structure to build this up.
// We start by specifying the nominal state (as defined above), then the uncertainty in position and velocity
// in the RIC frame. We could also specify the Cr, Cd, and mass uncertainties, but these aren't accounted for until
// Nyx can also estimate the deviation of the spacecraft parameters.
let jwst_uncertainty = SpacecraftUncertainty::builder()
.nominal(jwst)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.3)
.z_km(1.5)
.vx_km_s(1e-4)
.vy_km_s(0.6e-3)
.vz_km_s(3e-3)
.build();
println!("{jwst_uncertainty}");
// Build the Kalman filter estimate.
// Note that we could have used the KfEstimate structure directly (as seen throughout the OD integration tests)
// but this approach requires quite a bit more boilerplate code.
let jwst_estimate = jwst_uncertainty.to_estimate()?;
// Set up the spacecraft dynamics.
// We'll use the point masses of the Earth, Sun, Jupiter (barycenter, because it's in the DE440), and the Moon.
// We'll also enable solar radiation pressure since the James Webb has a huge and highly reflective sun shield.
let orbital_dyn = OrbitalDynamics::point_masses(vec![MOON, SUN, JUPITER_BARYCENTER]);
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
// Build the propagator set up to use for the whole analysis.
let setup = Propagator::default(dynamics);
// All of the analysis will use this duration.
let prediction_duration = 6.5 * Unit::Day;
// === Covariance mapping ===
// For the covariance mapping / prediction, we'll use the common orbit determination approach.
// This is done by setting up a spacecraft OD process, and predicting for the analysis duration.
let ckf = KF::no_snc(jwst_estimate);
// Build the propagation instance for the OD process.
let prop = setup.with(jwst.with_stm(), almanac.clone());
let mut odp = SpacecraftODProcess::ckf(prop, ckf, BTreeMap::new(), None, almanac.clone());
// Define the prediction step, i.e. how often we want to know the covariance.
let step = 1_i64.minutes();
// Finally, predict, and export the trajectory with covariance to a parquet file.
odp.predict_for(step, prediction_duration)?;
odp.to_parquet(
&TrackingDataArc::default(),
"./02_jwst_covar_map.parquet",
ExportCfg::default(),
)?;
// === Monte Carlo framework ===
// Nyx comes with a complete multi-threaded Monte Carlo frame. It's blazing fast.
let my_mc = MonteCarlo::new(
jwst, // Nominal state
jwst_estimate.to_random_variable()?,
"02_jwst".to_string(), // Scenario name
None, // No specific seed specified, so one will be drawn from the computer's entropy.
);
let num_runs = 5_000;
let rslts = my_mc.run_until_epoch(
setup,
almanac.clone(),
jwst.epoch() + prediction_duration,
num_runs,
);
assert_eq!(rslts.runs.len(), num_runs);
// Finally, export these results, computing the eclipse percentage for all of these results.
// For all of the resulting trajectories, we'll want to compute the percentage of penumbra and umbra.
let eclipse_loc = EclipseLocator::cislunar(almanac.clone());
let umbra_event = eclipse_loc.to_umbra_event();
let penumbra_event = eclipse_loc.to_penumbra_event();
rslts.to_parquet(
"02_jwst_monte_carlo.parquet",
Some(vec![&umbra_event, &penumbra_event]),
ExportCfg::default(),
almanac,
)?;
Ok(())
}
More examples
33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 56 57 58 59 60 61 62 63 64 65 66 67 68 69 70 71 72 73 74 75 76 77 78 79 80 81 82 83 84 85 86 87 88 89 90 91 92 93 94 95 96 97 98 99 100 101 102 103 104 105 106 107 108 109 110 111 112 113 114 115 116 117 118 119 120 121 122 123 124 125 126 127 128 129 130 131 132 133 134 135 136 137 138 139 140 141 142 143 144 145 146 147 148 149 150 151 152 153 154 155 156 157 158 159 160 161 162 163 164 165 166 167 168 169 170 171 172 173 174 175 176 177 178 179 180 181 182 183 184 185 186 187 188 189 190 191 192 193 194 195 196 197 198 199 200 201 202 203 204 205 206 207 208 209 210 211 212 213 214 215 216 217 218 219 220 221 222 223 224 225 226 227 228 229 230 231 232 233 234 235 236 237 238 239 240 241 242 243 244 245 246 247 248 249 250 251 252 253 254 255 256 257 258 259 260 261 262 263 264 265 266 267 268 269 270 271 272 273 274 275 276 277 278 279 280 281 282 283 284 285 286 287 288 289 290 291 292
fn main() -> Result<(), Box<dyn Error>> {
pel::init();
// ====================== //
// === ALMANAC SET UP === //
// ====================== //
// Dynamics models require planetary constants and ephemerides to be defined.
// Let's start by grabbing those by using ANISE's MetaAlmanac.
let data_folder: PathBuf = [env!("CARGO_MANIFEST_DIR"), "examples", "04_lro_od"]
.iter()
.collect();
let meta = data_folder.join("lro-dynamics.dhall");
// Load this ephem in the general Almanac we're using for this analysis.
let mut almanac = MetaAlmanac::new(meta.to_string_lossy().to_string())
.map_err(Box::new)?
.process(true)
.map_err(Box::new)?;
let mut moon_pc = almanac.planetary_data.get_by_id(MOON)?;
moon_pc.mu_km3_s2 = 4902.74987;
almanac.planetary_data.set_by_id(MOON, moon_pc)?;
let mut earth_pc = almanac.planetary_data.get_by_id(EARTH)?;
earth_pc.mu_km3_s2 = 398600.436;
almanac.planetary_data.set_by_id(EARTH, earth_pc)?;
// Save this new kernel for reuse.
// In an operational context, this would be part of the "Lock" process, and should not change throughout the mission.
almanac
.planetary_data
.save_as(&data_folder.join("lro-specific.pca"), true)?;
// Lock the almanac (an Arc is a read only structure).
let almanac = Arc::new(almanac);
// 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.96,
})
.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 traj_as_flown = Traj::from_bsp(
lro_frame,
MOON_J2000,
almanac.clone(),
sc_template,
5.seconds(),
Some(Epoch::from_str("2024-01-01 00:00:00 UTC")?),
Some(Epoch::from_str("2024-01-02 00:00:00 UTC")?),
Aberration::LT,
Some("LRO".to_string()),
)?;
println!("{traj_as_flown}");
// ====================== //
// === MODEL MATCHING === //
// ====================== //
// Set up the spacecraft dynamics.
// Specify that the orbital dynamics must account for the graviational pull of the Earth and the Sun.
// The gravity of the Moon will also be accounted for since the spaceraft in a lunar orbit.
let mut orbital_dyn = OrbitalDynamics::point_masses(vec![EARTH, SUN, JUPITER_BARYCENTER]);
// We want to include the spherical harmonics, so let's download the gravitational data from the Nyx Cloud.
// We're using the GRAIL JGGRX model.
let mut jggrx_meta = MetaFile {
uri: "http://public-data.nyxspace.com/nyx/models/Luna_jggrx_1500e_sha.tab.gz".to_string(),
crc32: Some(0x6bcacda8), // Specifying the CRC32 avoids redownloading it if it's cached.
};
// And let's download it if we don't have it yet.
jggrx_meta.process(true)?;
// Build the spherical harmonics.
// The harmonics must be computed in the body fixed frame.
// We're using the long term prediction of the Moon principal axes frame.
let moon_pa_frame = MOON_PA_FRAME.with_orient(31008);
// let moon_pa_frame = IAU_MOON_FRAME;
let sph_harmonics = Harmonics::from_stor(
almanac.frame_from_uid(moon_pa_frame)?,
HarmonicsMem::from_shadr(&jggrx_meta.uri, 80, 80, true)?,
);
// Include the spherical harmonics into the orbital dynamics.
orbital_dyn.accel_models.push(sph_harmonics);
// We define the solar radiation pressure, using the default solar flux and accounting only
// for the eclipsing caused by the Earth and Moon.
// Note that by default, enabling the SolarPressure model will also enable the estimation of the coefficient of reflectivity.
let srp_dyn = SolarPressure::new(vec![EARTH_J2000, MOON_J2000], almanac.clone())?;
// Finalize setting up the dynamics, specifying the force models (orbital_dyn) separately from the
// acceleration models (SRP in this case). Use `from_models` to specify multiple accel models.
let dynamics = SpacecraftDynamics::from_model(orbital_dyn, srp_dyn);
println!("{dynamics}");
// Now we can build the propagator.
let setup = Propagator::default_dp78(dynamics.clone());
// For reference, let's build the trajectory with Nyx's models from that LRO state.
let (sim_final, traj_as_sim) = setup
.with(*traj_as_flown.first(), almanac.clone())
.until_epoch_with_traj(traj_as_flown.last().epoch())?;
println!("SIM INIT: {:x}", traj_as_flown.first());
println!("SIM FINAL: {sim_final:x}");
// Compute RIC difference between SIM and LRO ephem
let sim_lro_delta = sim_final
.orbit
.ric_difference(&traj_as_flown.last().orbit)?;
println!("{traj_as_sim}");
println!(
"SIM v LRO - RIC Position (m): {:.3}",
sim_lro_delta.radius_km * 1e3
);
println!(
"SIM v LRO - RIC Velocity (m/s): {:.3}",
sim_lro_delta.velocity_km_s * 1e3
);
traj_as_sim.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_sim_truth_error.parquet",
ExportCfg::default(),
)?;
// ==================== //
// === OD SIMULATOR === //
// ==================== //
// After quite some time trying to exactly match the model, we still end up with an oscillatory difference on the order of 150 meters between the propagated state
// and the truth LRO state.
// Therefore, we will actually run an estimation from a dispersed LRO state.
// The sc_seed is the true LRO state from the BSP.
let sc_seed = *traj_as_flown.first();
// 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_named(ground_station_file)?;
// 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)?;
// Build the tracking arc simulation to generate a "standard measurement".
let mut trk = TrackingArcSim::<Spacecraft, GroundStation>::new(
devices.clone(),
traj_as_flown.clone(),
configs,
)?;
trk.build_schedule(almanac.clone())?;
let arc = trk.generate_measurements(almanac.clone())?;
// Save the simulated tracking data
arc.to_parquet_simple("./04_lro_simulated_tracking.parquet")?;
// We'll note that in our case, we have continuous coverage of LRO when the vehicle is not behind the Moon.
println!("{arc}");
// Now that we have simulated measurements, we'll run the orbit determination.
// ===================== //
// === OD ESTIMATION === //
// ===================== //
let sc = SpacecraftUncertainty::builder()
.nominal(sc_seed)
.frame(LocalFrame::RIC)
.x_km(0.5)
.y_km(0.5)
.z_km(0.5)
.vx_km_s(5e-3)
.vy_km_s(5e-3)
.vz_km_s(5e-3)
.build();
// Build the filter initial estimate, which we will reuse in the filter.
let initial_estimate = sc.to_estimate()?;
println!("== FILTER STATE ==\n{sc_seed:x}\n{initial_estimate}");
let kf = KF::new(
// Increase the initial covariance to account for larger deviation.
initial_estimate,
// Until https://github.com/nyx-space/nyx/issues/351, we need to specify the SNC in the acceleration of the Moon J2000 frame.
SNC3::from_diagonal(10 * Unit::Minute, &[1e-12, 1e-12, 1e-12]),
);
// We'll set up the OD process to reject measurements whose residuals are move than 3 sigmas away from what we expect.
let mut odp = SpacecraftODProcess::ckf(
setup.with(initial_estimate.state().with_stm(), almanac.clone()),
kf,
devices,
Some(ResidRejectCrit::default()),
almanac.clone(),
);
odp.process_arc(&arc)?;
let ric_err = traj_as_flown
.at(odp.estimates.last().unwrap().epoch())?
.orbit
.ric_difference(&odp.estimates.last().unwrap().orbital_state())?;
println!("== RIC at end ==");
println!("RIC Position (m): {}", ric_err.radius_km * 1e3);
println!("RIC Velocity (m/s): {}", ric_err.velocity_km_s * 1e3);
odp.to_parquet(&arc, "./04_lro_od_results.parquet", ExportCfg::default())?;
// In our case, we have the truth trajectory from NASA.
// So we can compute the RIC state difference between the real LRO ephem and what we've just estimated.
// Export the OD trajectory first.
let od_trajectory = odp.to_traj()?;
// Build the RIC difference.
od_trajectory.ric_diff_to_parquet(
&traj_as_flown,
"./04_lro_od_truth_error.parquet",
ExportCfg::default(),
)?;
Ok(())
}
Auto Trait Implementations§
impl<'a, D, MsrSize, Accel, K, Trk> Freeze for ODProcess<'a, D, MsrSize, Accel, K, Trk>
impl<'a, D, MsrSize, Accel, K, Trk> !RefUnwindSafe for ODProcess<'a, D, MsrSize, Accel, K, Trk>
impl<'a, D, MsrSize, Accel, K, Trk> !Send for ODProcess<'a, D, MsrSize, Accel, K, Trk>
impl<'a, D, MsrSize, Accel, K, Trk> !Sync for ODProcess<'a, D, MsrSize, Accel, K, Trk>
impl<'a, D, MsrSize, Accel, K, Trk> !Unpin for ODProcess<'a, D, MsrSize, Accel, K, Trk>
impl<'a, D, MsrSize, Accel, K, Trk> !UnwindSafe for ODProcess<'a, D, MsrSize, Accel, K, Trk>
Blanket Implementations§
Source§impl<T> BorrowMut<T> for Twhere
T: ?Sized,
impl<T> BorrowMut<T> for Twhere
T: ?Sized,
Source§fn borrow_mut(&mut self) -> &mut T
fn borrow_mut(&mut self) -> &mut T
§impl<T> Instrument for T
impl<T> Instrument for T
§fn instrument(self, span: Span) -> Instrumented<Self>
fn instrument(self, span: Span) -> Instrumented<Self>
§fn in_current_span(self) -> Instrumented<Self>
fn in_current_span(self) -> Instrumented<Self>
Source§impl<T> IntoEither for T
impl<T> IntoEither for T
Source§fn into_either(self, into_left: bool) -> Either<Self, Self>
fn into_either(self, into_left: bool) -> Either<Self, Self>
self
into a Left
variant of Either<Self, Self>
if into_left
is true
.
Converts self
into a Right
variant of Either<Self, Self>
otherwise. Read moreSource§fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
fn into_either_with<F>(self, into_left: F) -> Either<Self, Self>
self
into a Left
variant of Either<Self, Self>
if into_left(&self)
returns true
.
Converts self
into a Right
variant of Either<Self, Self>
otherwise. Read more§impl<T> Pointable for T
impl<T> Pointable for T
§impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
impl<SS, SP> SupersetOf<SS> for SPwhere
SS: SubsetOf<SP>,
§fn to_subset(&self) -> Option<SS>
fn to_subset(&self) -> Option<SS>
self
from the equivalent element of its
superset. Read more§fn is_in_subset(&self) -> bool
fn is_in_subset(&self) -> bool
self
is actually part of its subset T
(and can be converted to it).§fn to_subset_unchecked(&self) -> SS
fn to_subset_unchecked(&self) -> SS
self.to_subset
but without any property checks. Always succeeds.§fn from_subset(element: &SS) -> SP
fn from_subset(element: &SS) -> SP
self
to the equivalent element of its superset.