celestial_pointing/commands/
outl.rs

1use super::{Command, CommandOutput};
2use crate::error::{Error, Result};
3use crate::session::Session;
4
5pub struct Outl;
6
7impl Command for Outl {
8    fn name(&self) -> &str {
9        "OUTL"
10    }
11    fn description(&self) -> &str {
12        "Identify outlier observations"
13    }
14
15    fn execute(&self, session: &mut Session, args: &[&str]) -> Result<CommandOutput> {
16        if args.is_empty() {
17            return Err(Error::Parse("OUTL requires a sigma threshold".into()));
18        }
19        let threshold: f64 = args[0]
20            .parse()
21            .map_err(|e| Error::Parse(format!("invalid threshold: {}", e)))?;
22        let do_mask = args.get(1).is_some_and(|a| a.eq_ignore_ascii_case("M"));
23
24        let fit = session
25            .last_fit
26            .as_ref()
27            .ok_or_else(|| Error::Fit("no fit results available (run FIT first)".into()))?;
28        let rms = fit.sky_rms;
29        let cutoff = threshold * rms;
30
31        let lat = session.latitude();
32        let mut outliers: Vec<(usize, f64)> = Vec::new();
33
34        for (i, obs) in session.observations.iter().enumerate() {
35            if obs.masked {
36                continue;
37            }
38            let h = obs.commanded_ha.radians();
39            let dec = obs.catalog_dec.radians();
40            let pier = obs.pier_side.sign();
41            let (model_dh, model_dd) = session.model.apply_equatorial(h, dec, lat, pier);
42            let raw_dh = (obs.actual_ha - obs.commanded_ha).arcseconds();
43            let raw_dd = (obs.observed_dec - obs.catalog_dec).arcseconds();
44            let dh = raw_dh - model_dh;
45            let dd = raw_dd - model_dd;
46            let dx = dh * libm::cos(dec);
47            let dr = libm::sqrt(dx * dx + dd * dd);
48            if dr > cutoff {
49                outliers.push((i, dr));
50            }
51        }
52
53        if outliers.is_empty() {
54            return Ok(CommandOutput::Text(format!(
55                "No outliers (threshold {:.1} * {:.2}\" = {:.2}\")",
56                threshold, rms, cutoff
57            )));
58        }
59
60        let mut output = format!(
61            "Outliers (residual > {:.1} * {:.2}\" = {:.2}\"):\n",
62            threshold, rms, cutoff
63        );
64        for &(idx, dr) in &outliers {
65            output += &format!("  obs {:>4}: {:.1}\"\n", idx + 1, dr);
66        }
67
68        if do_mask {
69            for &(idx, _) in &outliers {
70                session.observations[idx].masked = true;
71            }
72            output += &format!("\nMasked {} observations", outliers.len());
73        } else {
74            output += &format!("\nUse OUTL {:.1} M to mask", threshold);
75        }
76
77        Ok(CommandOutput::Text(output))
78    }
79}
celestial_pointing/commands/outl.rs

celestial_pointing/commands/
outl.rs
