1use crate::observation::PierSide;
2use crate::session::Session;
3
4pub struct ObsResidual {
5 pub ha_deg: f64,
6 pub dec_deg: f64,
7 pub dh: f64,
8 pub dd: f64,
9 pub dx: f64,
10 pub dr: f64,
11 pub index: usize,
12 pub pier_east: bool,
13}
14
15pub fn compute_residuals(session: &Session) -> Vec<ObsResidual> {
16 let lat = session.latitude();
17 session
18 .observations
19 .iter()
20 .enumerate()
21 .filter(|(_, obs)| !obs.masked)
22 .map(|(i, obs)| build_residual(i, obs, &session.model, lat))
23 .collect()
24}
25
26fn build_residual(
27 index: usize,
28 obs: &crate::observation::Observation,
29 model: &crate::model::PointingModel,
30 lat: f64,
31) -> ObsResidual {
32 let h = obs.commanded_ha.radians();
33 let dec = obs.catalog_dec.radians();
34 let pier = obs.pier_side.sign();
35 let (model_dh, model_dd) = model.apply_equatorial(h, dec, lat, pier);
36 let raw_dh = (obs.actual_ha - obs.commanded_ha).arcseconds();
37 let raw_dd = (obs.observed_dec - obs.catalog_dec).arcseconds();
38 let dh = raw_dh - model_dh;
39 let dd = raw_dd - model_dd;
40 let dx = dh * libm::cos(dec);
41 let dr = libm::sqrt(dx * dx + dd * dd);
42 ObsResidual {
43 ha_deg: obs.commanded_ha.degrees(),
44 dec_deg: obs.catalog_dec.degrees(),
45 dh,
46 dd,
47 dx,
48 dr,
49 index,
50 pier_east: obs.pier_side == PierSide::East,
51 }
52}
53
54pub fn require_fit(session: &Session) -> crate::error::Result<()> {
55 if session.last_fit.is_none() {
56 return Err(crate::error::Error::Fit(
57 "no fit results - run FIT first".into(),
58 ));
59 }
60 Ok(())
61}
62
63#[cfg(test)]
64mod tests {
65 use super::*;
66 use crate::observation::{Observation, PierSide};
67 use celestial_core::Angle;
68
69 fn make_obs(
70 cmd_ha_arcsec: f64,
71 act_ha_arcsec: f64,
72 cat_dec_deg: f64,
73 obs_dec_deg: f64,
74 pier: PierSide,
75 masked: bool,
76 ) -> Observation {
77 Observation {
78 catalog_ra: Angle::from_hours(0.0),
79 catalog_dec: Angle::from_degrees(cat_dec_deg),
80 observed_ra: Angle::from_hours(0.0),
81 observed_dec: Angle::from_degrees(obs_dec_deg),
82 lst: Angle::from_hours(0.0),
83 commanded_ha: Angle::from_arcseconds(cmd_ha_arcsec),
84 actual_ha: Angle::from_arcseconds(act_ha_arcsec),
85 pier_side: pier,
86 masked,
87 }
88 }
89
90 #[test]
91 fn empty_session_returns_empty() {
92 let session = Session::new();
93 let residuals = compute_residuals(&session);
94 assert!(residuals.is_empty());
95 }
96
97 #[test]
98 fn masked_observations_excluded() {
99 let mut session = Session::new();
100 session
101 .observations
102 .push(make_obs(0.0, 100.0, 45.0, 45.0, PierSide::East, false));
103 session
104 .observations
105 .push(make_obs(0.0, 200.0, 30.0, 30.0, PierSide::East, true));
106 let residuals = compute_residuals(&session);
107 assert_eq!(residuals.len(), 1);
108 assert_eq!(residuals[0].index, 0);
109 }
110
111 #[test]
112 fn residual_no_model_equals_raw() {
113 let mut session = Session::new();
114 session
115 .observations
116 .push(make_obs(0.0, 3600.0, 0.0, 2.0, PierSide::East, false));
117 let residuals = compute_residuals(&session);
118 assert_eq!(residuals.len(), 1);
119 let r = &residuals[0];
120 assert_eq!(r.dh, 3600.0);
121 assert_eq!(r.dd, 7200.0);
122 assert_eq!(r.dx, 3600.0 * libm::cos(0.0_f64));
123 assert_eq!(r.dr, libm::sqrt(3600.0_f64.powi(2) + 7200.0_f64.powi(2)));
124 }
125
126 #[test]
127 fn residual_with_ih_model() {
128 let mut session = Session::new();
129 session
130 .observations
131 .push(make_obs(0.0, 100.0, 45.0, 45.0, PierSide::East, false));
132 session.model.add_term("IH").unwrap();
133 session.model.set_coefficients(&[-100.0]).unwrap();
134 let residuals = compute_residuals(&session);
135 let r = &residuals[0];
136 let dec_rad = 45.0_f64.to_radians();
137 let model_dh = 100.0;
138 let expected_dh = 100.0 - model_dh;
139 assert_eq!(r.dh, expected_dh);
140 assert_eq!(r.dx, expected_dh * libm::cos(dec_rad));
141 }
142
143 #[test]
144 fn pier_side_recorded() {
145 let mut session = Session::new();
146 session
147 .observations
148 .push(make_obs(0.0, 0.0, 0.0, 0.0, PierSide::East, false));
149 session
150 .observations
151 .push(make_obs(0.0, 0.0, 0.0, 0.0, PierSide::West, false));
152 let residuals = compute_residuals(&session);
153 assert!(residuals[0].pier_east);
154 assert!(!residuals[1].pier_east);
155 }
156
157 #[test]
158 fn index_tracks_original_position() {
159 let mut session = Session::new();
160 session
161 .observations
162 .push(make_obs(0.0, 0.0, 0.0, 0.0, PierSide::East, true));
163 session
164 .observations
165 .push(make_obs(0.0, 0.0, 0.0, 0.0, PierSide::East, false));
166 session
167 .observations
168 .push(make_obs(0.0, 0.0, 0.0, 0.0, PierSide::East, false));
169 let residuals = compute_residuals(&session);
170 assert_eq!(residuals.len(), 2);
171 assert_eq!(residuals[0].index, 1);
172 assert_eq!(residuals[1].index, 2);
173 }
174
175 #[test]
176 fn require_fit_no_fit() {
177 let session = Session::new();
178 assert!(require_fit(&session).is_err());
179 }
180
181 #[test]
182 fn require_fit_with_fit() {
183 let mut session = Session::new();
184 session.last_fit = Some(crate::solver::FitResult {
185 coefficients: vec![1.0],
186 sigma: vec![0.1],
187 sky_rms: 5.0,
188 term_names: vec!["IH".to_string()],
189 });
190 assert!(require_fit(&session).is_ok());
191 }
192
193 #[test]
194 fn ha_and_dec_degrees_populated() {
195 let mut session = Session::new();
196 let cmd_ha_arcsec = 3600.0 * 15.0;
197 session.observations.push(make_obs(
198 cmd_ha_arcsec,
199 cmd_ha_arcsec,
200 45.0,
201 45.0,
202 PierSide::East,
203 false,
204 ));
205 let residuals = compute_residuals(&session);
206 let r = &residuals[0];
207 assert_eq!(r.ha_deg, Angle::from_arcseconds(cmd_ha_arcsec).degrees());
208 assert_eq!(r.dec_deg, 45.0);
209 }
210
211 #[test]
212 fn dr_is_sqrt_dx2_dd2() {
213 let mut session = Session::new();
214 session
215 .observations
216 .push(make_obs(0.0, 3.0, 0.0, 4.0 / 3600.0, PierSide::East, false));
217 let residuals = compute_residuals(&session);
218 let r = &residuals[0];
219 let expected_dr = libm::sqrt(r.dx * r.dx + r.dd * r.dd);
220 assert_eq!(r.dr, expected_dr);
221 }
222}