Reenable obstruction in tracking data

examples/04_lro_od/dsn-network.yaml

@@ -15,7 +15,7 @@ DSS-65 Madrid:
       white_noise:
         mean: 0.0
         sigma: 50.0e-6 # 5 cm/s
-  light_time_correction: false
+  light_time_correction: true
   latitude_deg: 40.427222
   longitude_deg: 4.250556
   height_km: 0.834939
@@ -43,7 +43,7 @@ DSS-34 Canberra:
       white_noise:
         mean: 0.0
         sigma: 50.0e-6 # 5 cm/s
-  light_time_correction: false
+  light_time_correction: true
   measurement_types:
     - range_km
     - doppler_km_s
@@ -68,7 +68,7 @@ DSS-13 Goldstone:
       white_noise:
         mean: 0.0
         sigma: 50.0e-6 # 5 cm/s
-  light_time_correction: false
+  light_time_correction: true
   measurement_types:
     - range_km
     - doppler_km_s

examples/04_lro_od/main.rs

@@ -253,7 +253,7 @@ fn main() -> Result<(), Box<dyn Error>> {
         // 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-11, 1e-11, 1e-11]),
+        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 mover than 4 sigmas away from what we expect.

examples/04_lro_od/plot_od_rslt.py

@@ -4,9 +4,13 @@
 import numpy as np
 import plotly.graph_objects as go
 import plotly.express as px
+import click
 
-if __name__ == "__main__":
-    df = pl.read_parquet("output_data/ekf_rng_dpl_az_el_odp.parquet")
+
+@click.command
+@click.option("-p", "--path", type=str, default="./04_lro_od_results.parquet")
+def main(path: str):
+    df = pl.read_parquet(path)
 
     df = (
         df.with_columns(pl.col("Epoch (UTC)").str.to_datetime("%Y-%m-%dT%H:%M:%S%.f"))
@@ -49,68 +53,25 @@
     # Add scatter points
     fig_qq.add_trace(
         go.Scatter(
-            x=qq[0][0],
-            y=qq[0][1],
-            mode='markers',
-            name='Sample Data',
-            marker=dict(color='blue')
+            x=qq[0][0], y=qq[0][1], mode="markers", name="Residuals ratios (QQ)", marker=dict(color="blue")
         )
     )
 
     # Add the theoretical line
     fig_qq.add_trace(
-        go.Scatter(
-            x=x_qq,
-            y=y_qq,
-            mode='lines',
-            name='Theoretical Normal',
-            line=dict(color='red')
-        )
+        go.Scatter(x=x_qq, y=y_qq, mode="lines", name="Theoretical Normal", line=dict(color="red"))
     )
 
     # Update layout
     fig_qq.update_layout(
-        title='Normal Q-Q Plot',
-        xaxis_title='Theoretical Quantiles',
-        yaxis_title='Sample Quantiles',
+        title="Normal Q-Q Plot",
+        xaxis_title="Theoretical Quantiles",
+        yaxis_title="Sample Quantiles",
     )
 
     # Show QQ plot
     fig_qq.show()
 
-    # Create histogram with normal distribution overlay
-    hist_fig = px.histogram(
-        df,
-        x="Residual ratio",
-        color="Tracker",
-        marginal="rug",
-        hover_data=df.columns,
-    )
-
-    # Calculate normal distribution parameters
-    mean = residual_ratio.mean()
-    std = residual_ratio.std()
-    x_range = np.linspace(residual_ratio.min(), residual_ratio.max(), 100)
-    y_normal = stats.norm.pdf(x_range, mean, std)
-
-    # Scale the normal distribution to match histogram height
-    max_hist_height = 100
-    scaling_factor = max_hist_height / max(y_normal)
-    y_normal_scaled = y_normal * scaling_factor
-
-    # Add normal distribution curve
-    hist_fig.add_trace(
-        go.Scatter(
-            x=x_range,
-            y=y_normal_scaled,
-            name='Normal Distribution',
-            line=dict(color='red', width=2),
-        )
-    )
-
-    # Show histogram with normal overlay
-    hist_fig.show()
-
     px.histogram(
         df,
         x="Residual ratio",
@@ -163,7 +124,6 @@
         y=["Sigma Vx (RIC) (km/s)", "Sigma Vy (RIC) (km/s)", "Sigma Vz (RIC) (km/s)"],
     ).show()
 
-    raise AssertionError("stop")
     # Load the RIC diff.
     for fname, errname in [
         ("04_lro_od_truth_error", "OD vs Flown"),
@@ -219,3 +179,7 @@
             ],
             title=f"Velocity error with {errname} ({fname})",
         ).show()
+
+
+if __name__ == "__main__":
+    main()

examples/04_lro_od/requirements.txt

@@ -0,0 +1,38 @@
+asttokens==2.4.1
+attrs==23.2.0
+click==8.1.7
+decorator==5.1.1
+executing==2.0.1
+iniconfig==2.0.0
+ipython==8.25.0
+jedi==0.19.1
+matplotlib-inline==0.1.7
+maturin==1.6.0
+numpy==1.26.4
+packaging==24.1
+pandas==2.1.4
+parso==0.8.4
+pexpect==4.9.0
+pip==24.0
+plotly==5.16.1
+pluggy==1.5.0
+polars==0.20.31
+prompt-toolkit==3.0.47
+ptyprocess==0.7.0
+pure-eval==0.2.2
+pyarrow==13.0.0
+pygments==2.18.0
+pytest==7.2.2
+python-dateutil==2.9.0.post0
+python-slugify==8.0.4
+pytz==2024.1
+ruff==0.5.0
+scipy==1.11.4
+six==1.16.0
+stack-data==0.6.3
+tenacity==8.3.0
+text-unidecode==1.3
+traitlets==5.14.3
+typing-extensions==4.12.2
+tzdata==2024.1
+wcwidth==0.2.13

src/od/ground_station/trk_device.rs

@@ -84,6 +84,14 @@ impl TrackingDevice<Spacecraft> for GroundStation {
                         self.name, self.elevation_mask_deg, aer_t0.elevation_deg, aer_t1.elevation_deg
                     );
                     return Ok(None);
+                } else if aer_t0.is_obstructed() || aer_t1.is_obstructed() {
+                    debug!(
+                        "{} obstruction at t0={}, t1={} -- no measurement",
+                        self.name,
+                        aer_t0.is_obstructed(),
+                        aer_t1.is_obstructed()
+                    );
+                    return Ok(None);
                 }
 
                 // Noises are computed at the midpoint of the integration time.
@@ -128,7 +136,7 @@ impl TrackingDevice<Spacecraft> for GroundStation {
                 action: "computing AER",
             })?;
 
-        if aer.elevation_deg >= self.elevation_mask_deg {
+        if aer.elevation_deg >= self.elevation_mask_deg && !aer.is_obstructed() {
             // Only update the noises if the measurement is valid.
             let noises = self.noises(rx.orbit.epoch, rng)?;