Struct KfEstimate

pub struct KfEstimate<T: State>
where
    DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>,
    <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy,
    <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,
{
    pub nominal_state: T,
    pub state_deviation: OVector<f64, <T as State>::Size>,
    pub covar: OMatrix<f64, <T as State>::Size, <T as State>::Size>,
    pub covar_bar: OMatrix<f64, <T as State>::Size, <T as State>::Size>,
    pub predicted: bool,
    pub stm: OMatrix<f64, <T as State>::Size, <T as State>::Size>,
}

Kalman filter Estimate

Fields

nominal_state: T

The estimated state

state_deviation: OVector<f64, <T as State>::Size>

The state deviation

covar: OMatrix<f64, <T as State>::Size, <T as State>::Size>

The Covariance of this estimate

covar_bar: OMatrix<f64, <T as State>::Size, <T as State>::Size>

The predicted covariance of this estimate

predicted: bool

Whether or not this is a predicted estimate from a time update, or an estimate from a measurement

stm: OMatrix<f64, <T as State>::Size, <T as State>::Size>

The STM used to compute this Estimate

Implementations

impl<T: State> KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

pub fn from_covar( nominal_state: T, covar: OMatrix<f64, <T as State>::Size, <T as State>::Size>, ) -> Self

Initializes a new filter estimate from the nominal state (not dispersed) and the full covariance

pub fn from_diag( nominal_state: T, diag: OVector<f64, <T as State>::Size>, ) -> Self

Initializes a new filter estimate from the nominal state (not dispersed) and the diagonal of the covariance

impl KfEstimate<Spacecraft>

pub fn disperse_from_diag( nominal_state: Spacecraft, dispersions: Vec<StateDispersion>, seed: Option<u128>, ) -> Result<Self, Box<dyn Error>>

Generates an initial Kalman filter state estimate dispersed from the nominal state using the provided standard deviation parameters.

The resulting estimate will have a diagonal covariance matrix constructed from the variances of each parameter. Limitation: This method may not work correctly for all Keplerian orbital elements, refer to https://github.com/nyx-space/nyx/issues/339 for details.

pub fn to_random_variable(&self) -> Result<MultivariateNormal, Box<dyn Error>>

Builds a multivariate random variable from this estimate’s nominal state and covariance, zero mean.

Examples found in repository

examples/02_jwst_covar_monte_carlo/main.rs (line 132)

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, 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("./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(())
}

pub fn sigma_for(&self, param: StateParameter) -> Result<f64, AstroError>

Returns the 1-sigma uncertainty for a given parameter, in that parameter’s unit

This method uses the OrbitDual structure to compute the estimate in the hyperdual space and rotate the nominal covariance into that space.

pub fn keplerian_covar(&self) -> SMatrix<f64, 6, 6>

Returns the 6x6 covariance (i.e. square of the sigma/uncertainty) of the SMA, ECC, INC, RAAN, AOP, and True Anomaly.

Trait Implementations

impl<T: Clone + State> Clone for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

fn clone(&self) -> KfEstimate<T>

Returns a copy of the value.

fn clone_from(&mut self, source: &Self)

Performs copy-assignment from source.

impl<T: Debug + State> Debug for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<T: State> Display for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<T: State> Estimate<T> for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

fn zeros(nominal_state: T) -> Self

An empty estimate. This is useful if wanting to store an estimate outside the scope of a filtering loop.

fn nominal_state(&self) -> T

The nominal state as reported by the filter dynamics

fn state_deviation(&self) -> OVector<f64, <T as State>::Size>

The state deviation as computed by the filter.

fn covar(&self) -> OMatrix<f64, <T as State>::Size, <T as State>::Size>

The Covariance of this estimate. Will return the predicted covariance if this is a time update/prediction.

fn predicted_covar( &self, ) -> OMatrix<f64, <T as State>::Size, <T as State>::Size>

The predicted covariance of this estimate from the time update

fn predicted(&self) -> bool

Whether or not this is a predicted estimate from a time update, or an estimate from a measurement

fn stm(&self) -> &OMatrix<f64, <T as State>::Size, <T as State>::Size>

The STM used to compute this Estimate

fn set_state_deviation(&mut self, new_state: OVector<f64, <T as State>::Size>)

Sets the state deviation.

fn set_covar( &mut self, new_covar: OMatrix<f64, <T as State>::Size, <T as State>::Size>, )

Sets the Covariance of this estimate

fn epoch(&self) -> Epoch

Epoch of this Estimate

fn set_epoch(&mut self, dt: Epoch)

fn state(&self) -> T

The estimated state

fn within_sigma(&self, sigma: f64) -> bool

Returns whether this estimate is within some bound The 68-95-99.7 rule is a good way to assess whether the filter is operating normally

fn within_3sigma(&self) -> bool

Returns whether this estimate is within 3 sigma, which represent 99.7% for a Normal distribution

impl<T: State> LowerExp for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

fn fmt(&self, f: &mut Formatter<'_>) -> Result

Formats the value using the given formatter.

impl<T: State> Mul<f64> for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

type Output = KfEstimate<T>

The resulting type after applying the * operator.

fn mul(self, rhs: f64) -> Self::Output

Performs the * operation.

impl NavSolution<Spacecraft> for KfEstimate<Spacecraft>

fn orbital_state(&self) -> Orbit

fn expected_state(&self) -> Orbit

Returns the nominal state as computed by the dynamics

impl<T: PartialEq + State> PartialEq for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

fn eq(&self, other: &KfEstimate<T>) -> bool

Tests for self and other values to be equal, and is used by ==.

fn ne(&self, other: &Rhs) -> bool

Tests for !=. The default implementation is almost always sufficient, and should not be overridden without very good reason.

impl<T: Copy + State> Copy for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,

impl<T: State> StructuralPartialEq for KfEstimate<T>

where DefaultAllocator: Allocator<<T as State>::Size> + Allocator<<T as State>::Size, <T as State>::Size> + Allocator<<T as State>::VecLength>, <DefaultAllocator as Allocator<<T as State>::Size>>::Buffer<f64>: Copy, <DefaultAllocator as Allocator<<T as State>::Size, <T as State>::Size>>::Buffer<f64>: Copy,