1use super::{MountTypeFlags, Term};
2use crate::error::{Error, Result};
3
4#[derive(Debug, Clone, Copy, PartialEq, Eq)]
5pub enum TrigFunc {
6 Sin,
7 Cos,
8}
9
10#[derive(Debug, Clone, Copy, PartialEq, Eq)]
11pub enum Coordinate {
12 H,
13 D,
14 A,
15 E,
16 Z,
17}
18
19#[derive(Debug, Clone, Copy, PartialEq, Eq)]
20pub enum ResultType {
21 H,
22 D,
23 X,
24 P,
25 A,
26 E,
27 Z,
28 S,
29}
30
31#[derive(Debug, Clone)]
32pub struct HarmonicComponent {
33 pub func: TrigFunc,
34 pub coord: Coordinate,
35 pub frequency: u8,
36}
37
38#[derive(Debug, Clone)]
39pub struct HarmonicTerm {
40 pub result: ResultType,
41 pub components: Vec<HarmonicComponent>,
42 pub original_name: String,
43}
44
45pub fn parse_harmonic(name: &str) -> Result<HarmonicTerm> {
46 let chars: Vec<char> = name.chars().collect();
47 if chars.len() < 4 || chars[0] != 'H' {
48 return Err(Error::InvalidHarmonic(format!(
49 "too short or missing H prefix: {}",
50 name
51 )));
52 }
53 let result = parse_result_type(chars[1], name)?;
54 let components = parse_components(&chars[2..], name)?;
55 if components.is_empty() {
56 return Err(Error::InvalidHarmonic(format!(
57 "no components found: {}",
58 name
59 )));
60 }
61 Ok(HarmonicTerm {
62 result,
63 components,
64 original_name: name.to_string(),
65 })
66}
67
68fn parse_result_type(ch: char, name: &str) -> Result<ResultType> {
69 match ch {
70 'H' => Ok(ResultType::H),
71 'D' => Ok(ResultType::D),
72 'X' => Ok(ResultType::X),
73 'P' => Ok(ResultType::P),
74 'A' => Ok(ResultType::A),
75 'E' => Ok(ResultType::E),
76 'Z' => Ok(ResultType::Z),
77 'S' => Ok(ResultType::S),
78 _ => Err(Error::InvalidHarmonic(format!(
79 "unknown result type '{}' in {}",
80 ch, name
81 ))),
82 }
83}
84
85fn parse_func(ch: char, name: &str) -> Result<TrigFunc> {
86 match ch {
87 'S' => Ok(TrigFunc::Sin),
88 'C' => Ok(TrigFunc::Cos),
89 _ => Err(Error::InvalidHarmonic(format!(
90 "unknown function '{}' in {}",
91 ch, name
92 ))),
93 }
94}
95
96fn parse_coord(ch: char, name: &str) -> Result<Coordinate> {
97 match ch {
98 'H' => Ok(Coordinate::H),
99 'D' => Ok(Coordinate::D),
100 'A' => Ok(Coordinate::A),
101 'E' => Ok(Coordinate::E),
102 'Z' => Ok(Coordinate::Z),
103 _ => Err(Error::InvalidHarmonic(format!(
104 "unknown coordinate '{}' in {}",
105 ch, name
106 ))),
107 }
108}
109
110fn parse_components(chars: &[char], name: &str) -> Result<Vec<HarmonicComponent>> {
111 let mut components = Vec::new();
112 let mut i = 0;
113 while i < chars.len() {
114 if i + 1 >= chars.len() {
115 return Err(Error::InvalidHarmonic(format!(
116 "incomplete component in {}",
117 name
118 )));
119 }
120 let func = parse_func(chars[i], name)?;
121 let coord = parse_coord(chars[i + 1], name)?;
122 i += 2;
123 let frequency = if i < chars.len() && chars[i].is_ascii_digit() {
124 let d = chars[i].to_digit(10).unwrap() as u8;
125 i += 1;
126 d
127 } else {
128 1
129 };
130 if frequency == 0 {
131 return Err(Error::InvalidHarmonic(format!(
132 "frequency 0 not allowed in {}",
133 name
134 )));
135 }
136 components.push(HarmonicComponent {
137 func,
138 coord,
139 frequency,
140 });
141 }
142 Ok(components)
143}
144
145impl HarmonicTerm {
146 fn evaluate_equatorial(&self, h: f64, dec: f64) -> f64 {
147 self.components.iter().fold(1.0, |acc, comp| {
148 let coord_val = match comp.coord {
149 Coordinate::H => h,
150 Coordinate::D => dec,
151 Coordinate::A | Coordinate::E | Coordinate::Z => 0.0,
152 };
153 acc * trig_value(comp.func, coord_val, comp.frequency)
154 })
155 }
156
157 fn evaluate_altaz(&self, az: f64, el: f64) -> f64 {
158 self.components.iter().fold(1.0, |acc, comp| {
159 let coord_val = match comp.coord {
160 Coordinate::A => az,
161 Coordinate::E => el,
162 Coordinate::Z => std::f64::consts::FRAC_PI_2 - el,
163 Coordinate::H | Coordinate::D => 0.0,
164 };
165 acc * trig_value(comp.func, coord_val, comp.frequency)
166 })
167 }
168}
169
170fn trig_value(func: TrigFunc, coord_val: f64, frequency: u8) -> f64 {
171 let arg = coord_val * frequency as f64;
172 match func {
173 TrigFunc::Sin => libm::sin(arg),
174 TrigFunc::Cos => libm::cos(arg),
175 }
176}
177
178impl Term for HarmonicTerm {
179 fn name(&self) -> &str {
180 &self.original_name
181 }
182 fn description(&self) -> &str {
183 "Harmonic correction term"
184 }
185
186 fn pier_sensitive(&self) -> bool {
187 matches!(self.result, ResultType::P)
188 }
189
190 fn applicable_mounts(&self) -> MountTypeFlags {
191 match self.result {
192 ResultType::H | ResultType::D | ResultType::X | ResultType::P => {
193 MountTypeFlags::EQUATORIAL
194 }
195 ResultType::A | ResultType::E | ResultType::Z | ResultType::S => MountTypeFlags::ALTAZ,
196 }
197 }
198
199 fn jacobian_equatorial(&self, h: f64, dec: f64, _lat: f64, pier: f64) -> (f64, f64) {
200 let val = self.evaluate_equatorial(h, dec);
201 match self.result {
202 ResultType::H => (val, 0.0),
203 ResultType::D => (0.0, val),
204 ResultType::X => (val / libm::cos(dec), 0.0),
205 ResultType::P => (-pier * libm::tan(dec) * val, 0.0),
206 _ => (0.0, 0.0),
207 }
208 }
209
210 fn jacobian_altaz(&self, az: f64, el: f64, _lat: f64) -> (f64, f64) {
211 let val = self.evaluate_altaz(az, el);
212 match self.result {
213 ResultType::A => (val, 0.0),
214 ResultType::E => (0.0, val),
215 ResultType::Z => (0.0, -val),
216 ResultType::S => (val / libm::cos(el), 0.0),
217 _ => (0.0, 0.0),
218 }
219 }
220}
221
222#[cfg(test)]
223mod tests {
224 use super::*;
225 use std::f64::consts::{FRAC_PI_2, FRAC_PI_4};
226
227 #[test]
228 fn parse_hdsh() {
229 let term = parse_harmonic("HDSH").unwrap();
230 assert_eq!(term.result, ResultType::D);
231 assert_eq!(term.components.len(), 1);
232 assert_eq!(term.components[0].func, TrigFunc::Sin);
233 assert_eq!(term.components[0].coord, Coordinate::H);
234 assert_eq!(term.components[0].frequency, 1);
235 }
236
237 #[test]
238 fn parse_hxch3() {
239 let term = parse_harmonic("HXCH3").unwrap();
240 assert_eq!(term.result, ResultType::X);
241 assert_eq!(term.components.len(), 1);
242 assert_eq!(term.components[0].func, TrigFunc::Cos);
243 assert_eq!(term.components[0].coord, Coordinate::H);
244 assert_eq!(term.components[0].frequency, 3);
245 }
246
247 #[test]
248 fn parse_hdch2cd4() {
249 let term = parse_harmonic("HDCH2CD4").unwrap();
250 assert_eq!(term.result, ResultType::D);
251 assert_eq!(term.components.len(), 2);
252 assert_eq!(term.components[0].func, TrigFunc::Cos);
253 assert_eq!(term.components[0].coord, Coordinate::H);
254 assert_eq!(term.components[0].frequency, 2);
255 assert_eq!(term.components[1].func, TrigFunc::Cos);
256 assert_eq!(term.components[1].coord, Coordinate::D);
257 assert_eq!(term.components[1].frequency, 4);
258 }
259
260 #[test]
261 fn parse_too_short() {
262 assert!(parse_harmonic("HD").is_err());
263 }
264
265 #[test]
266 fn parse_bad_prefix() {
267 assert!(parse_harmonic("XDSH").is_err());
268 }
269
270 #[test]
271 fn parse_bad_result_type() {
272 assert!(parse_harmonic("HQSH").is_err());
273 }
274
275 #[test]
276 fn parse_bad_function() {
277 assert!(parse_harmonic("HDXH").is_err());
278 }
279
280 #[test]
281 fn parse_bad_coordinate() {
282 assert!(parse_harmonic("HDSQ").is_err());
283 }
284
285 #[test]
286 fn parse_zero_frequency() {
287 assert!(parse_harmonic("HDSH0").is_err());
288 }
289
290 #[test]
291 fn hdsh_jacobian_at_pi_over_2() {
292 let term = parse_harmonic("HDSH").unwrap();
293 let (dh, dd) = term.jacobian_equatorial(FRAC_PI_2, 0.0, 0.0, 1.0);
294 assert_eq!(dh, 0.0);
295 assert_eq!(dd, 1.0);
296 }
297
298 #[test]
299 fn hdsh_jacobian_at_zero() {
300 let term = parse_harmonic("HDSH").unwrap();
301 let (dh, dd) = term.jacobian_equatorial(0.0, 0.0, 0.0, 1.0);
302 assert_eq!(dh, 0.0);
303 assert_eq!(dd, 0.0);
304 }
305
306 #[test]
307 fn hhch_jacobian_at_zero() {
308 let term = parse_harmonic("HHCH").unwrap();
309 let (dh, dd) = term.jacobian_equatorial(0.0, 0.0, 0.0, 1.0);
310 assert_eq!(dh, 1.0);
311 assert_eq!(dd, 0.0);
312 }
313
314 #[test]
315 fn hxch_at_zero_ha_zero_dec() {
316 let term = parse_harmonic("HXCH").unwrap();
317 let (dh, dd) = term.jacobian_equatorial(0.0, 0.0, 0.0, 1.0);
318 assert_eq!(dh, 1.0);
319 assert_eq!(dd, 0.0);
320 }
321
322 #[test]
323 fn hp_result_pier_sensitive() {
324 let term = parse_harmonic("HPSH").unwrap();
325 assert!(term.pier_sensitive());
326 }
327
328 #[test]
329 fn hd_result_not_pier_sensitive() {
330 let term = parse_harmonic("HDSH").unwrap();
331 assert!(!term.pier_sensitive());
332 }
333
334 #[test]
335 fn equatorial_result_mount_flags() {
336 for prefix in &["HH", "HD", "HX", "HP"] {
337 let name = format!("{}SH", prefix);
338 let term = parse_harmonic(&name).unwrap();
339 assert_eq!(term.applicable_mounts(), MountTypeFlags::EQUATORIAL);
340 }
341 }
342
343 #[test]
344 fn altaz_result_mount_flags() {
345 for prefix in &["HA", "HE", "HZ", "HS"] {
346 let name = format!("{}SA", prefix);
347 let term = parse_harmonic(&name).unwrap();
348 assert_eq!(term.applicable_mounts(), MountTypeFlags::ALTAZ);
349 }
350 }
351
352 #[test]
353 fn altaz_harmonic_jacobian() {
354 let term = parse_harmonic("HASA").unwrap();
355 let az = FRAC_PI_2;
356 let (da, de) = term.jacobian_altaz(az, 0.0, 0.0);
357 assert_eq!(da, 1.0);
358 assert_eq!(de, 0.0);
359 }
360
361 #[test]
362 fn altaz_elevation_result() {
363 let term = parse_harmonic("HESE").unwrap();
364 let el = FRAC_PI_4;
365 let (da, de) = term.jacobian_altaz(0.0, el, 0.0);
366 assert_eq!(da, 0.0);
367 assert_eq!(de, libm::sin(el));
368 }
369
370 #[test]
371 fn altaz_zenith_result() {
372 let term = parse_harmonic("HZSE").unwrap();
373 let el = FRAC_PI_4;
374 let (da, de) = term.jacobian_altaz(0.0, el, 0.0);
375 assert_eq!(da, 0.0);
376 assert_eq!(de, -libm::sin(el));
377 }
378
379 #[test]
380 fn two_component_product() {
381 let term = parse_harmonic("HDCH2CD4").unwrap();
382 let h = 0.0;
383 let dec = 0.0;
384 let (dh, dd) = term.jacobian_equatorial(h, dec, 0.0, 1.0);
385 assert_eq!(dh, 0.0);
386 assert_eq!(dd, 1.0);
387 }
388
389 #[test]
390 fn equatorial_harmonic_returns_zero_for_altaz() {
391 let term = parse_harmonic("HDSH").unwrap();
392 let (da, de) = term.jacobian_altaz(1.0, 0.5, 0.7);
393 assert_eq!((da, de), (0.0, 0.0));
394 }
395
396 #[test]
397 fn altaz_harmonic_returns_zero_for_equatorial() {
398 let term = parse_harmonic("HASA").unwrap();
399 let (dh, dd) = term.jacobian_equatorial(1.0, 0.5, 0.7, 1.0);
400 assert_eq!((dh, dd), (0.0, 0.0));
401 }
402
403 #[test]
404 fn original_name_preserved() {
405 let term = parse_harmonic("HDSH").unwrap();
406 assert_eq!(term.name(), "HDSH");
407 }
408
409 #[test]
410 fn frequency_3_multiplies_argument() {
411 let term = parse_harmonic("HHSH3").unwrap();
412 let h = FRAC_PI_2;
413 let (dh, _dd) = term.jacobian_equatorial(h, 0.0, 0.0, 1.0);
414 assert_eq!(dh, libm::sin(3.0 * h));
415 }
416}