Skip to main content

fusion_nav/eskf/
predict.rs

1//! Propagation, equations (9)–(22), and the coast of (22′) across a gap.
2//!
3//! Entry point: [`Eskf::predict`].
4
5use crate::health::Propagation;
6use crate::propagate::{self, ImuSample, Propagated, propagate};
7
8use super::Eskf;
9
10impl Eskf {
11    /// Propagate the nominal state and covariance across one IMU sample. Equations (9)–(22).
12    ///
13    /// The hot path, called at IMU rate. The step is `Δt`, the time from the filter's
14    /// [`time`](Self::time) to the sample's; the filter never reads a clock. The state
15    /// integrates the sample's increments over their own intervals, and `Δt` is what
16    /// everything that is about time passing reads: the health timers, and the test for a gap.
17    ///
18    /// A `Δt` longer than [`Config::max_predict_dt`](crate::Config::max_predict_dt) is not
19    /// integrated: one IMU sample cannot describe a long interval, and propagating it anyway
20    /// would put a number in the state that looks like an estimate and is not. The filter
21    /// coasts across it instead, equation (22′), as [`Propagation::Coasted`]: position moves
22    /// on the estimated velocity and the covariance grows by what
23    /// [`Config::coast`](crate::Config::coast) allows, so the first fix after the gap is
24    /// judged against an uncertainty that grew with it. With coasting off the step is refused,
25    /// as [`Propagation::StepTooLong`]. Either way the timers advance, so [`Status`](crate::Status)
26    /// degrades on schedule. It is `Δt` that is tested rather than an interval, because a gap is
27    /// time no sample describes: a driver that integrated across a stall hands over an increment
28    /// that does describe it, and a logger that dropped samples hands over one that does not.
29    ///
30    /// A sample no later than the filter's time is refused before the timers move at all,
31    /// and leaves the clock where it was.
32    ///
33    /// A sample carrying a NaN or an infinity is refused too, as
34    /// [`Propagation::NotFinite`]: propagating it would put the NaN in the quaternion and
35    /// then in the covariance, where nothing reports it and it never leaves. So is one whose
36    /// integration interval is under a microsecond or longer than `max_predict_dt`, as
37    /// [`Propagation::InvalidInterval`]. The timers advance in both cases, since `Δt` was fine
38    /// and only the sample was not.
39    ///
40    /// A step that comes out of (11)–(14) or (22) non-finite is discarded rather than
41    /// stored, as [`Propagation::StateNotFinite`], and the state and the covariance are
42    /// discarded together: a finite sample is not enough to guarantee a finite result, since
43    /// f32 has a finite range and both `a_n Δt` and `F P Fᵀ` can leave it. Committing one
44    /// half would leave the filter reporting an estimate whose uncertainty describes a
45    /// different step.
46    ///
47    /// The covariance only grows here. (22) adds `Q` and (20) spreads what is already there;
48    /// nothing in propagation takes uncertainty back out, which is the measurement update's
49    /// job. So an unaided filter's [`Validity`](crate::Validity) flags go false in the order their
50    /// variances cross [`Config::accuracy`](crate::Config::accuracy), and [`Status`](crate::Status)
51    /// does not follow: it reads the aiding timers and an alignment that has already latched.
52    /// The one exception is the position hold of [`Config::hold`](crate::Config::hold), which a
53    /// committed step runs while nothing aids the filter.
54    pub fn predict(&mut self, imu: ImuSample) -> Propagation {
55        let outcome = self.propagate_or_coast(imu);
56        // After the step returns rather than inside it, so that the hold's update sits beside
57        // the propagation's frame instead of above it ([measured]).
58        //
59        // [measured]: https://github.com/wboayue/fusion-nav/blob/main/DESIGN.md#measured-cost-by-function
60        match outcome {
61            Propagation::Propagated => self.step_yaw_estimator(imu),
62            Propagation::Coasted { .. } => self.restart_yaw_estimator(),
63            // The clock did not move: nothing was missed.
64            Propagation::NotInitialized | Propagation::InvalidStep { .. } => return outcome,
65            // It did, and the yaw estimator integrated none of it.
66            Propagation::StepTooLong { .. }
67            | Propagation::NotFinite
68            | Propagation::InvalidInterval { .. }
69            | Propagation::StateNotFinite => {
70                self.restart_yaw_estimator();
71                return outcome;
72            }
73        }
74        self.hold_if_unaided();
75        outcome
76    }
77
78    /// [`predict`](Self::predict) without the hold: refuse, coast or integrate the step.
79    #[inline(never)]
80    fn propagate_or_coast(&mut self, imu: ImuSample) -> Propagation {
81        if !self.initialized {
82            return Propagation::NotInitialized;
83        }
84        let dt = imu.time.since(self.time);
85        // Rejected ahead of the bookkeeping: a negative `dt` would wind the timers back.
86        if !dt.is_usable_step() {
87            return self.refuse_step(Propagation::InvalidStep { dt });
88        }
89
90        // Past here the time genuinely passed, so the clock and the health bookkeeping are
91        // real even when the propagation itself is refused.
92        self.time = imu.time;
93        self.diagnostics.advance(dt);
94
95        let limit = self.config.max_predict_dt;
96        if dt > limit {
97            // Noted here, whatever becomes of the step: nothing else records how far the
98            // interval ran, and a coast discarded as non-finite is a gap all the same.
99            self.diagnostics.propagation.note_gap(dt);
100            let Some(coast) = self.config.coast else {
101                return self.refuse_step(Propagation::StepTooLong { dt, limit });
102            };
103            // The sample is not read, so a non-finite one does not stop a coast.
104            let coasted = propagate::coast(
105                *self.estimate.state(),
106                self.covariance,
107                self.offset,
108                dt,
109                &self.config,
110                &coast,
111            );
112            return self.commit_step(coasted, Propagation::Coasted { dt });
113        }
114
115        // Tested after the gap rather than before it, so that a gap is still measured
116        // whatever the sample holds. A sensor producing NaN produces it again on the next
117        // step, where the count picks it up.
118        if !imu.is_finite() {
119            return self.refuse_step(Propagation::NotFinite);
120        }
121        if let Some(interval) = imu.unusable_interval() {
122            return self.refuse_step(Propagation::InvalidInterval { interval });
123        }
124        // An interval past the limit is the gap test's case arriving as one increment, a
125        // driver reporting milliseconds as seconds as readily as a stall: one sample does not
126        // describe it, whichever of the two it was.
127        if imu.longest_interval() > limit {
128            let interval = imu.longest_interval();
129            return self.refuse_step(Propagation::InvalidInterval { interval });
130        }
131        let propagated = propagate(
132            *self.estimate.state(),
133            self.covariance,
134            self.offset,
135            imu,
136            &self.config,
137        );
138        self.commit_step(propagated, Propagation::Propagated)
139    }
140
141    /// Commit a propagated or coasted step, or discard it whole if it came out non-finite, and
142    /// record the outcome.
143    ///
144    /// Propagated into a local first, a coast as much as a step: (11)–(14) and (22) can both
145    /// overflow f32 on finite input, and a state written before it is checked is one the
146    /// filter has already published. One tail for both, so neither can commit a state without
147    /// its covariance or skip [`note_alignment`](Self::note_alignment). Taken by each branch
148    /// rather than after them: carrying the step out of a branch as a value measured a
149    /// 960-byte copy of `P` in `predict`'s frame.
150    fn commit_step(&mut self, propagated: Propagated, outcome: Propagation) -> Propagation {
151        if !propagated.is_finite() {
152            return self.refuse_step(Propagation::StateNotFinite);
153        }
154        self.estimate.step(self.time, propagated.state);
155        self.angular_rate = propagated.omega;
156        self.commit_covariance(propagated.covariance, propagated.offset);
157        self.note_alignment();
158        self.diagnostics.propagation.record(outcome);
159        outcome
160    }
161
162    /// Record a refused step and hand the outcome back, as [`refuse`](super::fuse::refuse) does for
163    /// a measurement.
164    ///
165    /// Same reason: one place maps an outcome to what
166    /// [`PropagationHealth`](crate::PropagationHealth) counts, so a guard added to
167    /// [`predict`](Self::predict) cannot forget to count it. Whether the timers moved is
168    /// the guard's business, not this one's.
169    fn refuse_step(&mut self, outcome: Propagation) -> Propagation {
170        self.diagnostics.propagation.record(outcome);
171        outcome
172    }
173}
174
175#[cfg(test)]
176mod tests {
177
178    use crate::config::{Coast, Config, GRAVITY};
179    use crate::eskf::Eskf;
180    use crate::eskf::fixtures::*;
181
182    use crate::health::{Fusion, Propagation};
183    use crate::init::Alignment;
184    use crate::init::tests::still;
185
186    use crate::observation::mag::tests::attitude_of;
187
188    use crate::propagate::ImuSample;
189    use crate::state::{Covariance, ErrorState, STATES, State};
190    use crate::units::{
191        Acceleration, AngularRate, Attitude, Position, PositionNoise, Seconds, Velocity,
192    };
193
194    use nalgebra::Vector3;
195
196    fn elapsed(filter: &Eskf) -> f32 {
197        filter
198            .diagnostics()
199            .baro_altitude
200            .time_since_accepted
201            .expect("baro has been accepted")
202            .as_secs()
203    }
204
205    #[test]
206    fn predict_before_initialize_is_refused() {
207        let mut filter = Eskf::default();
208        assert_eq!(
209            filter.step(ImuSample::default(), DT),
210            Propagation::NotInitialized
211        );
212    }
213
214    #[test]
215    fn a_normal_step_propagates_and_advances_the_timers() {
216        let mut filter = aided();
217        assert_eq!(
218            filter.step(ImuSample::default(), DT),
219            Propagation::Propagated
220        );
221        assert_eq!(elapsed(&filter), DT.as_secs());
222    }
223
224    /// The `Eskf`-level check that (16)–(22) are wired at all: a step grows the uncertainty it
225    /// was initialized with. Without them a static start would hold `Initialization`'s sigmas
226    /// for the whole flight.
227    #[test]
228    fn a_step_grows_the_covariance() {
229        let mut filter = unheld();
230        let before = *filter.covariance();
231
232        assert_eq!(
233            filter.step(ImuSample::default(), DT),
234            Propagation::Propagated
235        );
236
237        let after = filter.covariance();
238        for state in [
239            ErrorState::PositionNorth,
240            ErrorState::VelocityNorth,
241            ErrorState::AttitudeX,
242            ErrorState::GyroBiasX,
243        ] {
244            assert!(
245                after.variance(state) > before.variance(state),
246                "{state:?}: {} did not grow from {}",
247                after.variance(state),
248                before.variance(state),
249            );
250        }
251    }
252
253    /// A refused step commits neither half. The state and the covariance advance together or
254    /// not at all: a state stored beside the covariance of a different step reports an estimate
255    /// whose uncertainty describes something else, which is worse than the refusal it replaces.
256    ///
257    /// The seed is finite and its diagonal positive, so `initialize_from` accepts it, and one
258    /// step of `F P Fᵀ` then leaves f32's range — the covariance's version of the overflow
259    /// (11)–(14) already had.
260    #[test]
261    fn an_overflowing_covariance_commits_neither_half() {
262        let mut filter = Eskf::default();
263        let seed = State {
264            velocity: Velocity::ned(1.0, 2.0, 3.0),
265            ..State::default()
266        };
267        let enormous = Covariance::from_matrix(
268            crate::state::CovarianceMatrix::from_diagonal_element(f32::MAX),
269        );
270        assert_eq!(filter.seed(seed, enormous), Ok(Alignment::Seeded));
271
272        assert_eq!(
273            filter.step(ImuSample::default(), DT),
274            Propagation::StateNotFinite
275        );
276        assert_eq!(filter.state().velocity, seed.velocity);
277        assert_eq!(*filter.covariance(), enormous);
278        assert_eq!(filter.diagnostics().propagation.refused_state_not_finite, 1);
279    }
280
281    /// The `Eskf`-level check that (11) is wired at all: a sample reporting *no* specific
282    /// force is a vehicle in free fall, whatever its attitude, because `a_n = R(q̂) 0 + g`
283    /// is gravity in any frame. One step and the estimate is falling at `γ Δt`.
284    ///
285    /// Written on the default sample rather than a plausible one for that reason — the
286    /// answer does not depend on what `aided()` happened to level to.
287    #[test]
288    fn a_step_with_no_specific_force_leaves_the_estimate_falling() {
289        let mut filter = aided();
290        assert_eq!(
291            filter.step(ImuSample::default(), DT),
292            Propagation::Propagated
293        );
294        let velocity = filter.state().velocity.vector();
295        let free_fall = GRAVITY * DT.as_secs();
296        assert!((velocity.z - free_fall).abs() < 1e-6, "{velocity:?}");
297        assert!(
298            velocity.x.abs() < 1e-6 && velocity.y.abs() < 1e-6,
299            "{velocity:?}"
300        );
301    }
302
303    /// A finite sample whose propagation overflows f32 is refused, and the estimate the
304    /// filter keeps is the last one that was a number.
305    ///
306    /// It takes accumulation rather than one step — `f32::MAX` of specific force is
307    /// `1.7e36` of velocity over 5 ms — so the loop runs until the refusal rather than
308    /// asserting on a step count, and the comparison is against the state immediately
309    /// before it.
310    #[test]
311    fn a_propagation_that_overflows_is_refused_and_the_estimate_is_left_alone() {
312        let mut filter = aided();
313        let imu = ImuSample::reading(
314            AngularRate::body(0.0, 0.0, 0.0),
315            Acceleration::body(f32::MAX, 0.0, -GRAVITY),
316        );
317        assert!(imu.is_finite());
318
319        let mut before = filter.state();
320        let mut steps = 0;
321        loop {
322            match filter.step(imu, DT) {
323                Propagation::Propagated => {
324                    before = filter.state();
325                    steps += 1;
326                    assert!(steps < 10_000, "never overflowed");
327                }
328                Propagation::StateNotFinite => break,
329                other => panic!("unexpected outcome {other:?}"),
330            }
331        }
332
333        assert_eq!(
334            filter.state(),
335            before,
336            "a poisoned state reached the filter"
337        );
338        assert!(filter.state().is_finite());
339        assert_eq!(
340            filter.diagnostics().propagation.refused_state_not_finite,
341            1,
342            "the refusal went uncounted"
343        );
344        // The refused step counts too: the time passed. Compared with a tolerance because
345        // the timer accumulates 0.005 a couple of hundred times in f32 while the
346        // right-hand side multiplies once.
347        let expected = DT.as_secs() * (steps + 1) as f32;
348        assert!(
349            (elapsed(&filter) - expected).abs() < 1e-3,
350            "{} vs {expected}: a refused step still happened in real time",
351            elapsed(&filter)
352        );
353    }
354
355    #[test]
356    fn a_sample_not_after_the_clock_is_refused_without_moving_the_timers_or_the_clock() {
357        let mut filter = aided();
358        let clock = filter.time();
359        for bad in [0.0, -0.01] {
360            let dt = Seconds::from_secs(bad);
361            assert_eq!(
362                filter.step(ImuSample::default(), dt),
363                Propagation::InvalidStep { dt },
364                "dt of {bad} should be refused"
365            );
366            assert_eq!(elapsed(&filter), 0.0, "dt of {bad} moved the timers");
367            assert_eq!(filter.time(), clock, "dt of {bad} moved the clock");
368        }
369    }
370
371    /// An interval is the sample's, not time passing: refused, but the clock and the timers
372    /// move as they do for a sample carrying a NaN.
373    #[test]
374    fn a_sample_with_an_unusable_interval_is_refused_after_the_time_passes() {
375        let mut filter = aided();
376        let before = filter.state();
377        let time = filter.time().expect("initialized").after(DT);
378        let backwards = Seconds::from_secs(-0.005);
379        let imu = ImuSample {
380            angle_interval: backwards,
381            ..still().imu.timed(time, DT)
382        };
383        assert_eq!(
384            filter.predict(imu),
385            Propagation::InvalidInterval {
386                interval: backwards
387            }
388        );
389        assert_eq!(filter.state(), before);
390        assert_eq!(filter.time(), Some(time));
391        assert!((elapsed(&filter) - DT.as_secs()).abs() < 1e-6);
392        assert_eq!(filter.diagnostics().propagation.refused_invalid, 1);
393    }
394
395    /// A gap is time no sample describes, and only the timestamps see it: a logger that
396    /// dropped samples hands over an increment integrated over one IMU interval, a second
397    /// after the last one.
398    #[test]
399    fn a_gap_is_read_off_the_timestamps_not_the_interval() {
400        let mut filter = aided();
401        let gap = Seconds::from_secs(1.0);
402        let time = filter.time().expect("initialized").after(gap);
403        assert_eq!(
404            filter.predict(still().imu.timed(time, DT)),
405            Propagation::Coasted { dt: gap }
406        );
407    }
408
409    /// Coasting off, so a long step is refused.
410    fn not_coasting() -> Config {
411        Config {
412            coast: None,
413            ..Config::default()
414        }
415    }
416
417    #[test]
418    fn a_step_over_the_limit_is_refused_with_coasting_off_but_the_time_still_passes() {
419        let mut filter = aided_with(not_coasting());
420        let before = filter.state();
421        // The worst SD-card dropout in the bundled corpus.
422        let dt = Seconds::from_secs(1.304);
423        assert!(matches!(
424            filter.step(ImuSample::default(), dt),
425            Propagation::StepTooLong { .. }
426        ));
427        assert_eq!(filter.state(), before, "a refused step moves nothing");
428        assert_eq!(
429            elapsed(&filter),
430            1.304,
431            "a refused step still happened in real time"
432        );
433    }
434
435    /// Level and still, seeded with a tilt the propagation can leak into velocity.
436    fn seeded_level(config: Config) -> Eskf {
437        let mut filter = Eskf::new(config).unwrap();
438        let mut sigmas = [0.5f32; STATES];
439        sigmas[ErrorState::AttitudeX.index()] = 0.05;
440        sigmas[ErrorState::AttitudeY.index()] = 0.05;
441        let _ = filter
442            .seed(State::default(), Covariance::from_sigmas(sigmas))
443            .expect("a sane seed");
444        filter
445    }
446
447    #[test]
448    fn a_vehicle_at_rest_under_its_own_gravity_does_not_fall() {
449        // (11): the specific force a still vehicle reads is the site's `γ`.
450        let reading = ImuSample::reading(
451            AngularRate::body(0.0, 0.0, 0.0),
452            Acceleration::body(0.0, 0.0, -SITE_GRAVITY),
453        );
454        let fallen = |config| {
455            let mut filter = seeded_level(config);
456            for _ in 0..100 {
457                assert!(
458                    filter
459                        .step(reading, Seconds::from_secs(0.01))
460                        .is_propagated()
461                );
462            }
463            filter.state().velocity.to_array()[2]
464        };
465        assert!(fallen(at_site()).abs() < 1e-5, "{}", fallen(at_site()));
466        let off = fallen(Config::default());
467        assert!((off - (GRAVITY - SITE_GRAVITY)).abs() < 1e-4, "{off}");
468    }
469
470    #[test]
471    fn a_coast_leaks_tilt_into_velocity_through_the_gravity_configured() {
472        // (22′) with no noise of its own: the horizontal velocity variance a coast adds is the
473        // tilt's, through (17)'s `γ`, so it scales as `γ²`.
474        let quiet = |config: Config| Config {
475            imu: crate::config::ImuNoise {
476                gyro_white: 0.0,
477                accel_white: 0.0,
478                gyro_bias_walk: 0.0,
479                accel_bias_walk: 0.0,
480            },
481            coast: Some(Coast {
482                acceleration: 0.0,
483                rotation: 0.0,
484            }),
485            ..config
486        };
487        let grown = |config| {
488            let mut filter = seeded_level(quiet(config));
489            let before = filter.covariance().variance(ErrorState::VelocityNorth);
490            assert!(matches!(
491                filter.step(still().imu, Seconds::from_secs(1.2)),
492                Propagation::Coasted { .. }
493            ));
494            filter.covariance().variance(ErrorState::VelocityNorth) - before
495        };
496        let ratio = grown(at_site()) / grown(Config::default());
497        let expected = (SITE_GRAVITY / GRAVITY) * (SITE_GRAVITY / GRAVITY);
498        // 1.1e-4 off, the rest of `F` over the gap; read against one `γ`, the ratio would be 1.
499        assert!(
500            (ratio - expected).abs() < 3e-4,
501            "{ratio} against {expected}"
502        );
503    }
504
505    #[test]
506    fn a_gap_is_coasted_on_the_estimated_velocity_and_the_time_still_passes() {
507        let mut filter = flying(Config::default());
508        let before = filter.state();
509        let gap = Seconds::from_secs(1.2);
510        assert_eq!(
511            filter.step(still().imu, gap),
512            Propagation::Coasted { dt: gap }
513        );
514
515        let after = filter.state();
516        let moved = after.position.vector() - before.position.vector();
517        assert!((moved.x - 24.0).abs() < 1e-3, "north {} m", moved.x);
518        assert!(moved.y.abs() < 1e-6 && moved.z.abs() < 1e-6, "{moved:?}");
519        assert_eq!(
520            after.velocity, before.velocity,
521            "nothing measured it change"
522        );
523        assert_eq!(after.attitude, before.attitude, "nothing measured it turn");
524        assert_eq!(elapsed(&filter), 1.2, "the gap happened in real time");
525
526        let propagation = filter.diagnostics().propagation;
527        assert_eq!(propagation.coasted, 1);
528        assert_eq!(propagation.refused_too_long, 0, "a coast is not a refusal");
529        assert_eq!(propagation.longest_gap, Some(gap));
530    }
531
532    #[test]
533    fn a_coast_discarded_as_non_finite_still_measures_the_gap() {
534        // Finite, so `Config::validate` passes it, and squared into `P` it overflows.
535        let mut filter = flying(Config {
536            coast: Some(Coast {
537                acceleration: f32::MAX,
538                rotation: 0.1,
539            }),
540            ..Config::default()
541        });
542        let before = filter.state();
543        let gap = Seconds::from_secs(1.2);
544        assert_eq!(filter.step(still().imu, gap), Propagation::StateNotFinite);
545        assert_eq!(filter.state(), before, "a discarded coast commits nothing");
546        let propagation = filter.diagnostics().propagation;
547        assert_eq!(propagation.longest_gap, Some(gap));
548        assert_eq!(propagation.coasted, 0);
549        assert_eq!(propagation.refused_state_not_finite, 1);
550    }
551
552    #[test]
553    fn a_coast_grows_velocity_by_the_unmeasured_acceleration() {
554        let acceleration = 3.0;
555        let gap = Seconds::from_secs(1.2);
556        let variances = |acceleration| {
557            let mut filter = flying(Config {
558                coast: Some(Coast {
559                    acceleration,
560                    rotation: 0.0,
561                }),
562                hold: None,
563                ..Config::default()
564            });
565            let _ = filter.step(still().imu, gap);
566            let p = *filter.covariance().as_matrix();
567            (
568                p[(
569                    ErrorState::VelocityNorth.index(),
570                    ErrorState::VelocityNorth.index(),
571                )],
572                p[(
573                    ErrorState::PositionNorth.index(),
574                    ErrorState::PositionNorth.index(),
575                )],
576            )
577        };
578        let (velocity, position) = variances(acceleration);
579        let (velocity_q, position_q) = variances(0.0);
580        // Exactly the white-noise integral over the gap: `a² Δt` on velocity, `a² Δt³ / 3`
581        // on position, whatever the step count.
582        let expected = acceleration * acceleration * gap.as_secs();
583        assert!(
584            ((velocity - velocity_q) - expected).abs() < 1e-3 * expected,
585            "velocity grew {} over (22) alone, expected {expected}",
586            velocity - velocity_q
587        );
588        let integrated = expected * gap.as_secs() * gap.as_secs() / 3.0;
589        let grown = position - position_q;
590        assert!(
591            (grown - integrated).abs() < 1e-3 * integrated,
592            "position grew {grown}, expected {integrated}"
593        );
594    }
595
596    #[test]
597    fn a_coast_grows_attitude_by_the_unmeasured_rotation() {
598        let rotation = 0.1;
599        let gap = Seconds::from_secs(1.2);
600        let variance = |rotation| {
601            let mut filter = flying(Config {
602                coast: Some(Coast {
603                    acceleration: 0.0,
604                    rotation,
605                }),
606                ..Config::default()
607            });
608            let _ = filter.step(still().imu, gap);
609            let p = *filter.covariance().as_matrix();
610            p[(ErrorState::AttitudeZ.index(), ErrorState::AttitudeZ.index())]
611        };
612        // The attitude block of `F` is the identity at `ω = 0`, so the density lands whole.
613        let grown = variance(rotation) - variance(0.0);
614        let expected = rotation * rotation * gap.as_secs();
615        assert!(
616            (grown - expected).abs() < 1e-3 * expected,
617            "attitude grew {grown} over (22) alone, expected {expected}"
618        );
619    }
620
621    #[test]
622    fn a_coast_walks_the_barometric_offset_across_the_gap() {
623        // (30′)'s offset drifts whether or not the IMU was logged.
624        let mut filter = flying(Config::default());
625        let before = filter.offset.variance;
626        let gap = Seconds::from_secs(1.2);
627        let _ = filter.step(still().imu, gap);
628        let walk = filter.config.baro_offset_walk;
629        let expected = walk * walk * gap.as_secs();
630        let grown = filter.offset.variance - before;
631        assert!(
632            (grown - expected).abs() < 1e-3 * expected,
633            "offset variance grew {grown}, expected {expected}"
634        );
635    }
636
637    #[test]
638    fn the_first_fix_after_a_coasted_gap_is_accepted_where_a_refused_one_is_rejected() {
639        // Where the vehicle is after 1.2 s at 20 m/s north.
640        let fix = Position::ned(24.0, 0.0, -100.0);
641        let noise = PositionNoise::horizontal_vertical(1.0, 1.0);
642        let gap = Seconds::from_secs(1.2);
643
644        let mut coasting = flying(Config::default());
645        let _ = coasting.step(still().imu, gap);
646        assert!(
647            coasting
648                .fuse_gnss_position(coasting.now(), fix, noise, Position::zero())
649                .horizontal
650                .is_accepted()
651        );
652
653        let mut refusing = flying(not_coasting());
654        let _ = refusing.step(still().imu, gap);
655        assert!(matches!(
656            refusing
657                .fuse_gnss_position(refusing.now(), fix, noise, Position::zero())
658                .horizontal,
659            Fusion::Rejected { .. }
660        ));
661    }
662
663    #[test]
664    fn a_non_finite_sample_is_refused_but_the_time_still_passes() {
665        let level = Acceleration::body(0.0, 0.0, -GRAVITY);
666        let still = AngularRate::body(0.0, 0.0, 0.0);
667        for (name, imu) in [
668            (
669                "NaN gyro",
670                ImuSample::reading(AngularRate::body(f32::NAN, 0.0, 0.0), level),
671            ),
672            (
673                "NaN accel",
674                ImuSample::reading(still, Acceleration::body(0.0, f32::NAN, -GRAVITY)),
675            ),
676            (
677                "infinite accel",
678                ImuSample::reading(still, Acceleration::body(0.0, 0.0, f32::INFINITY)),
679            ),
680        ] {
681            let mut filter = aided();
682            let before = filter.state();
683            assert_eq!(filter.step(imu, DT), Propagation::NotFinite, "{name}");
684            assert_eq!(filter.state(), before, "{name} reached the state");
685            assert_eq!(
686                elapsed(&filter),
687                DT.as_secs(),
688                "{name}: a refused step still happened in real time"
689            );
690            assert_eq!(
691                filter.diagnostics().propagation.refused_not_finite,
692                1,
693                "{name} went uncounted"
694            );
695        }
696    }
697
698    #[test]
699    fn a_long_step_carrying_a_non_finite_sample_still_measures_the_gap() {
700        // Both refusals apply; the gap is reported because nothing else records how far
701        // the interval ran, while the sensor fault recurs on the next step.
702        // A coast reads no sample, so the NaN does not stop it either.
703        let dt = Seconds::from_secs(1.304);
704        let imu = ImuSample::reading(
705            AngularRate::body(f32::NAN, 0.0, 0.0),
706            Acceleration::body(0.0, 0.0, -GRAVITY),
707        );
708        for (config, coasted) in [(Config::default(), true), (not_coasting(), false)] {
709            let mut filter = aided_with(config);
710            let outcome = filter.step(imu, dt);
711            assert_eq!(matches!(outcome, Propagation::Coasted { .. }), coasted);
712            assert_eq!(matches!(outcome, Propagation::StepTooLong { .. }), !coasted);
713            assert!(
714                filter
715                    .state()
716                    .position
717                    .vector()
718                    .iter()
719                    .all(|x| x.is_finite())
720            );
721            let propagation = filter.diagnostics().propagation;
722            assert_eq!(propagation.longest_gap, Some(dt));
723            assert_eq!(propagation.refused_not_finite, 0);
724        }
725    }
726
727    #[test]
728    fn propagation_refusals_are_counted_and_the_worst_gap_kept() {
729        let mut filter = aided_with(not_coasting());
730        for bad in [0.0, -0.01] {
731            assert!(
732                !filter
733                    .step(ImuSample::default(), Seconds::from_secs(bad))
734                    .is_propagated()
735            );
736        }
737        for gap in [0.34, 1.304, 0.5] {
738            assert!(
739                !filter
740                    .step(ImuSample::default(), Seconds::from_secs(gap))
741                    .is_propagated()
742            );
743        }
744
745        let propagation = filter.diagnostics().propagation;
746        assert_eq!(propagation.refused_invalid, 2);
747        assert_eq!(propagation.refused_too_long, 3);
748        assert_eq!(
749            propagation.longest_gap.map(Seconds::as_secs),
750            Some(1.304),
751            "the worst gap, not the last"
752        );
753    }
754
755    /// An interval longer than any step the filter integrates is refused, whatever the
756    /// timestamps say: one sample does not describe it.
757    #[test]
758    fn an_interval_longer_than_the_limit_is_refused() {
759        let mut filter = aided();
760        let before = filter.state();
761        let limit = filter.config().max_predict_dt;
762        let time = filter.now().after(DT);
763        let long = Seconds::from_secs(limit.as_secs() * 25.0);
764        let imu = ImuSample {
765            velocity_interval: long,
766            ..still().imu.timed(time, DT)
767        };
768        assert_eq!(
769            filter.predict(imu),
770            Propagation::InvalidInterval { interval: long }
771        );
772        assert_eq!(filter.state(), before);
773    }
774
775    #[test]
776    fn the_angular_rate_is_the_last_integrated_sample_less_the_bias() {
777        let bias = AngularRate::body(0.01, -0.02, 0.005);
778        let mut filter = seeded(Attitude::level(), bias);
779        assert_eq!(filter.angular_rate(), None, "no step integrated yet");
780
781        let gyro = AngularRate::body(0.3, 0.1, -0.2);
782        let imu = ImuSample::reading(gyro, Acceleration::body(0.0, 0.0, -GRAVITY));
783        assert!(filter.step(imu, DT).is_propagated());
784        let omega = filter.angular_rate().expect("a step was integrated");
785        // Corrected as an increment and divided back out, so equal to rounding.
786        let expected = gyro.vector() - bias.vector();
787        assert!((omega.vector() - expected).norm() < 1e-6, "{omega:?}");
788
789        // A refused sample leaves the last rate, as it leaves the state.
790        let broken = imu.with_gyro(AngularRate::body(f32::NAN, 0.0, 0.0));
791        assert_eq!(filter.step(broken, DT), Propagation::NotFinite);
792        assert_eq!(filter.angular_rate(), Some(omega));
793
794        // A gap is coasted on no sample, so there is no rate until the next one.
795        let gap = Seconds::from_secs(1.0);
796        assert_eq!(filter.step(imu, gap), Propagation::Coasted { dt: gap });
797        assert_eq!(filter.angular_rate(), None);
798        assert!(filter.step(imu, DT).is_propagated());
799        assert!(filter.angular_rate().is_some());
800
801        // A fresh start forgets a rate measured against the last one's bias.
802        let _ = filter
803            .initialize_over(&[still(); 8], Seconds::from_secs(0.25))
804            .expect("a 2 s window of stillness");
805        assert_eq!(filter.angular_rate(), None);
806    }
807
808    #[test]
809    fn the_angular_rate_corrects_a_lever_arm_on_a_rolled_vehicle() {
810        // Nose east and rolled right a quarter turn, so the right wing points down. Worked by
811        // hand rather than through the rotation, so a transposed R cannot agree with itself:
812        // an antenna 1 m forward sits 1 m east of the IMU; a rate about body right is a rate
813        // about down, which turns east toward south, so that antenna moves south.
814        use core::f32::consts::FRAC_PI_2;
815        let mut filter = seeded(attitude_of(FRAC_PI_2, 0.0, FRAC_PI_2), AngularRate::zero());
816        let imu = ImuSample::reading(
817            AngularRate::body(0.0, 0.5, 0.0),
818            Acceleration::body(0.0, -GRAVITY, 0.0),
819        );
820        assert!(filter.step(imu, DT).is_propagated());
821
822        let r = Vector3::new(1.0, 0.0, 0.0);
823        let rotation = filter.state().attitude.quaternion();
824        let omega = filter.angular_rate().expect("a step was integrated");
825        // To within the 5 mrad the one step turned it.
826        let arm = rotation * r;
827        assert!((arm - Vector3::new(0.0, 1.0, 0.0)).norm() < 1e-2, "{arm}");
828        let lever = rotation * omega.vector().cross(&r);
829        assert!(
830            (lever - Vector3::new(-0.5, 0.0, 0.0)).norm() < 1e-2,
831            "{lever}"
832        );
833    }
834}