Fix find_arcs and iteration convergence

src/md/events/search.rs

@@ -355,12 +355,39 @@ where
     /// - Handles edge cases where the trajectory starts or ends with a rising or falling edge.
     /// - Prints debug information for each event and arc.
     ///
-
+    /// ## Note
+    /// If no zero crossing happens in the trajectory, i.e. the there is "event is true" _and_ "event is false",
+    /// then this function checks whether the event is true at the start and end of the trajectory. If so, it means
+    /// that there is a single arc that spans the whole trajectory.
+    ///
     pub fn find_arcs<E>(&self, event: &E) -> Result<Vec<EventArc<S>>, NyxError>
     where
         E: EventEvaluator<S>,
     {
-        let mut events = self.find(event)?;
+        let mut events = match self.find(event) {
+            Ok(events) => events,
+            Err(_) => {
+                // We haven't found the start or end of an arc, i.e. no zero crossing on the event.
+                // However, if the trajectory start and end are above the event value, then we found an arc.
+                let first_eval = event.eval(self.first());
+                let last_eval = event.eval(self.last());
+                if first_eval > 0.0 && last_eval > 0.0 {
+                    // No event crossing found, but from the start until the end of the trajectory, we're in the same arc
+                    // because the evaluation of the event is above the zero crossing.
+                    // Hence, there's a single arc, and it's from start until the end of the trajectory.
+                    vec![
+                        EventDetails::new(*self.first(), first_eval, event, self)?,
+                        EventDetails::new(*self.last(), last_eval, event, self)?,
+                    ]
+                } else {
+                    return Err(NyxError::from(TrajError::EventNotFound {
+                        start: self.first().epoch(),
+                        end: self.last().epoch(),
+                        event: format!("{event}"),
+                    }));
+                }
+            }
+        };
         events.sort_by_key(|event| event.state.epoch());
 
         // Now, let's pair the events.

src/od/ground_station.rs

@@ -405,12 +405,12 @@ impl fmt::Display for GroundStation {
     fn fmt(&self, f: &mut fmt::Formatter) -> fmt::Result {
         write!(
             f,
-            "[{}] {} (lat.: {:.4} deg    long.: {:.4} deg    alt.: {:.3} m)",
-            self.frame,
+            "{} (lat.: {:.4} deg    long.: {:.4} deg    alt.: {:.3} m) [{}]",
             self.name,
             self.latitude_deg,
             self.longitude_deg,
             self.height_km * 1e3,
+            self.frame,
         )
     }
 }

src/od/process/conf.rs

@@ -56,9 +56,9 @@ impl fmt::Display for SmoothingArc {
 pub struct IterationConf {
     /// The number of measurements to account for in the iteration
     pub smoother: SmoothingArc,
-    /// The absolute tolerance of the RMS postfit residual
+    /// The absolute tolerance of the RMS prefit residual ratios
     pub absolute_tol: f64,
-    /// The relative tolerance between the latest RMS postfit residual and the best RMS postfit residual so far
+    /// The relative tolerance between the latest RMS prefit residual ratios and the best RMS prefit residual ratios so far
     pub relative_tol: f64,
     /// The maximum number of iterations to allow (will raise an error if the filter has not converged after this many iterations)
     pub max_iterations: usize,
@@ -149,8 +149,8 @@ impl Default for IterationConf {
     fn default() -> Self {
         Self {
             smoother: SmoothingArc::All,
-            absolute_tol: 1e-2,
-            relative_tol: 1e-3,
+            absolute_tol: 1e-1,
+            relative_tol: 1e-2,
             max_iterations: 15,
             max_divergences: 3,
             force_failure: false,

src/od/process/mod.rs

@@ -314,7 +314,7 @@ where
             }
 
             info!("***************************");
-            info!("*** Iteration number {} ***", iter_cnt);
+            info!("*** Iteration number {iter_cnt:02} ***");
             info!("***************************");
 
             // First, smooth the estimates
@@ -331,27 +331,32 @@ where
 
             // Compute the new RMS
             let new_rms = self.rms_residual_ratios();
-            let cur_rel_rms = (new_rms - best_rms).abs() / best_rms;
-            if cur_rel_rms < config.relative_tol {
+            let cur_rms_num = (new_rms - previous_rms).abs();
+            let cur_rel_rms = cur_rms_num / previous_rms;
+            if cur_rel_rms < config.relative_tol || cur_rms_num < config.absolute_tol * best_rms {
                 info!("*****************");
                 info!("*** CONVERGED ***");
                 info!("*****************");
                 info!(
                     "New residual RMS: {:.5}\tPrevious RMS: {:.5}\tBest RMS: {:.5}",
                     new_rms, previous_rms, best_rms
                 );
-                info!(
-                    "Filter converged to relative tolerance ({:.2e} < {:.2e}) after {} iterations",
-                    cur_rel_rms, config.relative_tol, iter_cnt
-                );
-
+                if cur_rel_rms < config.relative_tol {
+                    info!(
+                        "Filter converged on relative tolerance ({:.2e} < {:.2e}) after {} iterations",
+                        cur_rel_rms, config.relative_tol, iter_cnt
+                    );
+                } else {
+                    info!(
+                        "Filter converged on relative change ({:.2e} < {:.2e} * {:.2e}) after {} iterations",
+                        cur_rms_num, config.absolute_tol, best_rms, iter_cnt
+                    );
+                }
                 break;
-            }
-
-            if new_rms > previous_rms {
+            } else if new_rms > previous_rms {
                 warn!(
-                    "New residual RMS: {:.5}\tPrevious RMS: {:.5}\tBest RMS: {:.5}",
-                    new_rms, previous_rms, best_rms
+                    "New residual RMS: {:.5}\tPrevious RMS: {:.5}\tBest RMS: {:.5} ({cur_rel_rms:.2e} > {:.2e})",
+                    new_rms, previous_rms, best_rms, config.relative_tol
                 );
                 divergence_cnt += 1;
                 previous_rms = new_rms;
@@ -371,8 +376,8 @@ where
                 }
             } else {
                 info!(
-                    "New residual RMS: {:.5}\tPrevious RMS: {:.5}\tBest RMS: {:.5}",
-                    new_rms, previous_rms, best_rms
+                    "New residual RMS: {:.5}\tPrevious RMS: {:.5}\tBest RMS: {:.5} ({cur_rel_rms:.2e} > {:.2e})",
+                    new_rms, previous_rms, best_rms, config.relative_tol
                 );
                 // Reset the counter
                 divergence_cnt = 0;
@@ -411,8 +416,7 @@ where
         let mut devices = arc.rebuild_devices::<S, Dev>(self.cosm.clone()).unwrap();
 
         let measurements = &arc.measurements;
-        let step_size = Duration::ZERO;
-        match arc.min_duration_sep() {
+        let step_size = match arc.min_duration_sep() {
             Some(step_size) => step_size,
             None => {
                 return Err(NyxError::CustomError(

src/od/simulator/arc.rs

@@ -279,51 +279,56 @@ impl TrackingArcSim<Orbit, RangeDoppler, GroundStation> {
                 // Convert the trajectory into the ground station frame.
                 let traj = self.trajectory.to_frame(device.frame, cosm.clone())?;
 
-                let elevation_arcs = traj.find_arcs(&device).unwrap();
-                for arc in elevation_arcs {
-                    info!("Working on {arc}");
-                    let strand_start = arc.rise.state.epoch();
-                    let strand_end = arc.fall.state.epoch();
-
-                    let mut strand_range = EpochRanges {
-                        start: strand_start,
-                        end: strand_end,
-                    };
-
-                    if let Cadence::Intermittent { on, off } = scheduler.cadence {
-                        // Check that the next start time is within the allocated time
-                        if let Some(prev_strand) = built_cfg[name].strands.as_ref().unwrap().last()
-                        {
-                            if prev_strand.end + off > strand_range.start {
-                                // We're turning on the tracking sooner than the schedule allows, so let's fix that.
-                                strand_range.start = prev_strand.end + off;
-                                // Check that we didn't eat into the whole tracking opportunity
-                                if strand_range.start > strand_end {
-                                    // Lost this whole opportunity.
-                                    info!("Discarding {name} opportunity from {} to {} due to cadence {:?}", strand_start, strand_end, scheduler.cadence);
+                match traj.find_arcs(&device) {
+                    Err(_) => info!("No measurements from {name}"),
+                    Ok(elevation_arcs) => {
+                        for arc in elevation_arcs {
+                            info!("Working on {arc}");
+                            let strand_start = arc.rise.state.epoch();
+                            let strand_end = arc.fall.state.epoch();
+
+                            let mut strand_range = EpochRanges {
+                                start: strand_start,
+                                end: strand_end,
+                            };
+
+                            if let Cadence::Intermittent { on, off } = scheduler.cadence {
+                                // Check that the next start time is within the allocated time
+                                if let Some(prev_strand) =
+                                    built_cfg[name].strands.as_ref().unwrap().last()
+                                {
+                                    if prev_strand.end + off > strand_range.start {
+                                        // We're turning on the tracking sooner than the schedule allows, so let's fix that.
+                                        strand_range.start = prev_strand.end + off;
+                                        // Check that we didn't eat into the whole tracking opportunity
+                                        if strand_range.start > strand_end {
+                                            // Lost this whole opportunity.
+                                            info!("Discarding {name} opportunity from {} to {} due to cadence {:?}", strand_start, strand_end, scheduler.cadence);
+                                        }
+                                    }
+                                }
+                                // Check that we aren't tracking for longer than configured
+                                if strand_range.end - strand_range.start > on {
+                                    strand_range.end = strand_range.start + on;
                                 }
                             }
+
+                            // We've found when the spacecraft is below the horizon, so this is a new strand.
+                            built_cfg
+                                .get_mut(name)
+                                .unwrap()
+                                .strands
+                                .as_mut()
+                                .unwrap()
+                                .push(strand_range);
                         }
-                        // Check that we aren't tracking for longer than configured
-                        if strand_range.end - strand_range.start > on {
-                            strand_range.end = strand_range.start + on;
-                        }
-                    }
 
-                    // We've found when the spacecraft is below the horizon, so this is a new strand.
-                    built_cfg
-                        .get_mut(name)
-                        .unwrap()
-                        .strands
-                        .as_mut()
-                        .unwrap()
-                        .push(strand_range);
+                        info!(
+                            "Built {} tracking strands for {name}",
+                            built_cfg[name].strands.as_ref().unwrap().len()
+                        );
+                    }
                 }
-
-                info!(
-                    "Built {} tracking strands for {name}",
-                    built_cfg[name].strands.as_ref().unwrap().len()
-                );
             }
         }
         // todo!("remove overlaps")

tests/orbit_determination/xhat_dev.rs

@@ -71,7 +71,7 @@ fn xhat_dev_test_ekf_two_body() {
 
     let (err_p, err_v) = rss_orbit_errors(&initial_state_dev, &initial_state);
     println!(
-        "Initial state dev: {:.3} m\t{:.3} m/s\n{}",
+        "Initial state dev: {:.3} m\t{:.3} m/s\nDelta: {}",
         err_p * 1e3,
         err_v * 1e3,
         initial_state - initial_state_dev
@@ -85,6 +85,7 @@ fn xhat_dev_test_ekf_two_body() {
         .unwrap();
 
     // Simulate tracking data
+    println!("{traj}");
     let mut arc_sim = TrackingArcSim::with_seed(all_stations, traj.clone(), configs, 0).unwrap();
     arc_sim.build_schedule(cosm.clone()).unwrap();