1use super::record::{EopParameters, EopRecord};
2use crate::{CoordError, CoordResult};
3
4#[derive(Debug, Clone, Copy, PartialEq, Eq)]
5pub enum InterpolationMethod {
6 Linear,
7
8 Lagrange5,
9}
10
11pub struct EopInterpolator {
12 records: Vec<EopRecord>,
13
14 method: InterpolationMethod,
15
16 max_gap_days: f64,
17}
18
19impl EopInterpolator {
20 pub fn new(mut records: Vec<EopRecord>) -> Self {
21 records.sort_by(|a, b| a.mjd.partial_cmp(&b.mjd).unwrap());
22
23 Self {
24 records,
25 method: InterpolationMethod::Linear,
26 max_gap_days: 5.0,
27 }
28 }
29
30 pub fn with_method(mut self, method: InterpolationMethod) -> Self {
31 self.method = method;
32 self
33 }
34
35 pub fn with_max_gap(mut self, max_gap_days: f64) -> Self {
36 self.max_gap_days = max_gap_days;
37 self
38 }
39
40 pub fn get(&self, mjd: f64) -> CoordResult<EopParameters> {
41 if self.records.is_empty() {
42 return Err(CoordError::data_unavailable(
43 "No EOP records available for interpolation",
44 ));
45 }
46
47 if let Ok(idx) = self
48 .records
49 .binary_search_by(|r| r.mjd.partial_cmp(&mjd).unwrap())
50 {
51 let mut params = self.records[idx].to_parameters();
52 params.compute_s_prime();
53 return Ok(params);
54 }
55
56 let (before_idx, after_idx) = self.find_interpolation_interval(mjd)?;
57
58 let gap = self.records[after_idx].mjd - self.records[before_idx].mjd;
59 if gap > self.max_gap_days {
60 return Err(CoordError::data_unavailable(format!(
61 "Gap of {:.1} days exceeds maximum interpolation gap of {:.1} days",
62 gap, self.max_gap_days
63 )));
64 }
65
66 match self.method {
67 InterpolationMethod::Linear => self.linear_interpolate(mjd, before_idx, after_idx),
68 InterpolationMethod::Lagrange5 => self.lagrange_interpolate(mjd, 5),
69 }
70 }
71
72 fn find_interpolation_interval(&self, mjd: f64) -> CoordResult<(usize, usize)> {
73 if mjd < self.records[0].mjd {
74 return Err(CoordError::data_unavailable(format!(
75 "MJD {:.1} is before first available record (MJD {:.1})",
76 mjd, self.records[0].mjd
77 )));
78 }
79
80 if mjd > self.records.last().unwrap().mjd {
81 return Err(CoordError::data_unavailable(format!(
82 "MJD {:.1} is after last available record (MJD {:.1})",
83 mjd,
84 self.records.last().unwrap().mjd
85 )));
86 }
87
88 let mut left = 0;
89 let mut right = self.records.len() - 1;
90
91 while right - left > 1 {
92 let mid = (left + right) / 2;
93 if self.records[mid].mjd <= mjd {
94 left = mid;
95 } else {
96 right = mid;
97 }
98 }
99
100 Ok((left, right))
101 }
102
103 fn linear_interpolate(
104 &self,
105 mjd: f64,
106 before_idx: usize,
107 after_idx: usize,
108 ) -> CoordResult<EopParameters> {
109 let r1 = &self.records[before_idx];
110 let r2 = &self.records[after_idx];
111
112 let t = (mjd - r1.mjd) / (r2.mjd - r1.mjd);
113
114 let p1 = r1.to_parameters();
115 let p2 = r2.to_parameters();
116
117 let x_p = p1.x_p + t * (p2.x_p - p1.x_p);
118 let y_p = p1.y_p + t * (p2.y_p - p1.y_p);
119 let ut1_utc = p1.ut1_utc + t * (p2.ut1_utc - p1.ut1_utc);
120 let lod = p1.lod + t * (p2.lod - p1.lod);
121
122 let dx = match (p1.dx, p2.dx) {
123 (Some(dx1), Some(dx2)) => Some(dx1 + t * (dx2 - dx1)),
124 _ => None,
125 };
126
127 let dy = match (p1.dy, p2.dy) {
128 (Some(dy1), Some(dy2)) => Some(dy1 + t * (dy2 - dy1)),
129 _ => None,
130 };
131
132 let xrt = match (p1.xrt, p2.xrt) {
133 (Some(v1), Some(v2)) => Some(v1 + t * (v2 - v1)),
134 _ => None,
135 };
136
137 let yrt = match (p1.yrt, p2.yrt) {
138 (Some(v1), Some(v2)) => Some(v1 + t * (v2 - v1)),
139 _ => None,
140 };
141
142 let mut params = EopParameters {
143 mjd,
144 x_p,
145 y_p,
146 ut1_utc,
147 lod,
148 dx,
149 dy,
150 xrt,
151 yrt,
152 s_prime: 0.0,
153 flags: p1.flags,
154 };
155
156 params.compute_s_prime();
157
158 Ok(params)
159 }
160
161 fn lagrange_interpolate(&self, mjd: f64, n: usize) -> CoordResult<EopParameters> {
162 if n > self.records.len() {
163 return Err(CoordError::invalid_coordinate(format!(
164 "Not enough records for {}-point Lagrange interpolation",
165 n
166 )));
167 }
168
169 let center_idx = self.find_center_index(mjd)?;
170 let half_n = n / 2;
171
172 let start_idx = center_idx.saturating_sub(half_n);
173 let end_idx = (start_idx + n).min(self.records.len());
174 let actual_start = end_idx.saturating_sub(n);
175
176 if end_idx - actual_start < n {
177 return Err(CoordError::data_unavailable(
178 "Insufficient records for Lagrange interpolation",
179 ));
180 }
181
182 let points: Vec<_> = self.records[actual_start..actual_start + n]
183 .iter()
184 .map(|r| r.to_parameters())
185 .collect();
186
187 let x_p = self.lagrange_interpolate_value(mjd, &points, |p| p.x_p);
188 let y_p = self.lagrange_interpolate_value(mjd, &points, |p| p.y_p);
189 let ut1_utc = self.lagrange_interpolate_value(mjd, &points, |p| p.ut1_utc);
190 let lod = self.lagrange_interpolate_value(mjd, &points, |p| p.lod);
191
192 let dx = if points.iter().all(|p| p.dx.is_some()) {
193 Some(self.lagrange_interpolate_value(mjd, &points, |p| p.dx.unwrap()))
194 } else {
195 None
196 };
197
198 let dy = if points.iter().all(|p| p.dy.is_some()) {
199 Some(self.lagrange_interpolate_value(mjd, &points, |p| p.dy.unwrap()))
200 } else {
201 None
202 };
203
204 let xrt = if points.iter().all(|p| p.xrt.is_some()) {
205 Some(self.lagrange_interpolate_value(mjd, &points, |p| p.xrt.unwrap()))
206 } else {
207 None
208 };
209
210 let yrt = if points.iter().all(|p| p.yrt.is_some()) {
211 Some(self.lagrange_interpolate_value(mjd, &points, |p| p.yrt.unwrap()))
212 } else {
213 None
214 };
215
216 let mut params = EopParameters {
217 mjd,
218 x_p,
219 y_p,
220 ut1_utc,
221 lod,
222 dx,
223 dy,
224 xrt,
225 yrt,
226 s_prime: 0.0,
227 flags: points[0].flags,
228 };
229
230 params.compute_s_prime();
231
232 Ok(params)
233 }
234
235 fn find_center_index(&self, mjd: f64) -> CoordResult<usize> {
236 let idx = self
237 .records
238 .binary_search_by(|r| r.mjd.partial_cmp(&mjd).unwrap());
239
240 match idx {
241 Ok(i) => Ok(i),
242 Err(i) => {
243 if i == 0 {
244 Ok(0)
245 } else if i >= self.records.len() {
246 Ok(self.records.len() - 1)
247 } else {
248 let before = (self.records[i - 1].mjd - mjd).abs();
249 let after = (self.records[i].mjd - mjd).abs();
250 if before <= after {
251 Ok(i - 1)
252 } else {
253 Ok(i)
254 }
255 }
256 }
257 }
258 }
259
260 fn lagrange_interpolate_value<F>(&self, mjd: f64, points: &[EopParameters], extract: F) -> f64
261 where
262 F: Fn(&EopParameters) -> f64,
263 {
264 let n = points.len();
265 let mut result = 0.0;
266
267 for i in 0..n {
268 let yi = extract(&points[i]);
269 let xi = points[i].mjd;
270
271 let mut li = 1.0;
272 for (j, point) in points.iter().enumerate().take(n) {
273 if i != j {
274 let xj = point.mjd;
275 li *= (mjd - xj) / (xi - xj);
276 }
277 }
278
279 result += yi * li;
280 }
281
282 result
283 }
284
285 pub fn time_span(&self) -> Option<(f64, f64)> {
286 if self.records.is_empty() {
287 None
288 } else {
289 Some((self.records[0].mjd, self.records.last().unwrap().mjd))
290 }
291 }
292
293 pub fn record_count(&self) -> usize {
294 self.records.len()
295 }
296
297 pub fn extend(&mut self, records: Vec<EopRecord>) {
298 self.records.extend(records);
299 self.records
300 .sort_by(|a, b| a.mjd.partial_cmp(&b.mjd).unwrap());
301 self.records.dedup_by(|a, b| a.mjd == b.mjd);
302 }
303}
304
305#[cfg(test)]
306mod tests {
307 use super::*;
308
309 fn create_test_records() -> Vec<EopRecord> {
310 let mut records = Vec::new();
311
312 for i in 0..5 {
313 let mjd = 59945.0 + i as f64;
314 let x_p = 0.1 + 0.001 * i as f64;
315 let y_p = 0.2 + 0.002 * i as f64;
316 let ut1_utc = 0.01 + 0.0001 * i as f64;
317 let lod = 0.001 + 0.00001 * i as f64;
318
319 let record = EopRecord::new(mjd, x_p, y_p, ut1_utc, lod).unwrap();
320 records.push(record);
321 }
322
323 records
324 }
325
326 #[test]
327 fn test_linear_interpolation() {
328 let records = create_test_records();
329 let interpolator = EopInterpolator::new(records);
330
331 let mjd = 59946.5;
332 let params = interpolator.get(mjd).unwrap();
333
334 let expected_x_p = (0.101 + 0.102) / 2.0;
335 let expected_y_p = (0.202 + 0.204) / 2.0;
336
337 assert!((params.x_p - expected_x_p).abs() < 1e-10);
338 assert!((params.y_p - expected_y_p).abs() < 1e-10);
339 assert_eq!(params.mjd, mjd);
340 }
341
342 #[test]
343 fn test_exact_match() {
344 let records = create_test_records();
345 let interpolator = EopInterpolator::new(records);
346
347 let mjd = 59947.0;
348 let params = interpolator.get(mjd).unwrap();
349
350 assert_eq!(params.mjd, mjd);
351 assert!((params.x_p - 0.102).abs() < 1e-10);
352 assert!((params.y_p - 0.204).abs() < 1e-10);
353 }
354
355 #[test]
356 fn test_linear_interpolation_cip_offsets() {
357 let mut records = Vec::new();
358
359 for i in 0..=1 {
360 let mjd = 59945.0 + i as f64;
361 let mut record = EopRecord::new(mjd, 0.1 + 0.001 * i as f64, 0.2, 0.01, 0.001).unwrap();
362 let dx = 1.0 + i as f64;
363 let dy = -0.2 - 0.1 * i as f64;
364 record = record.with_cip_offsets(dx, dy).unwrap();
365 records.push(record);
366 }
367
368 let interpolator = EopInterpolator::new(records);
369 let params = interpolator.get(59945.5).unwrap();
370
371 assert!(params.dx.is_some());
372 assert!(params.dy.is_some());
373 assert!((params.dx.unwrap() - 1.5).abs() < 1e-10);
374 assert!((params.dy.unwrap() - (-0.25)).abs() < 1e-10);
375 assert!(params.flags.has_cip_offsets);
376 }
377
378 #[test]
379 fn test_lagrange_interpolation() {
380 let records = create_test_records();
381 let interpolator =
382 EopInterpolator::new(records).with_method(InterpolationMethod::Lagrange5);
383
384 let mjd = 59947.0;
385 let params = interpolator.get(mjd).unwrap();
386
387 assert!((params.x_p - 0.102).abs() < 1e-10);
388 assert!((params.y_p - 0.204).abs() < 1e-10);
389 }
390
391 #[test]
392 fn test_lagrange_interpolation_cip_offsets() {
393 let mut records = Vec::new();
394 for i in 0..6 {
395 let mjd = 59945.0 + i as f64;
396 let mut record = EopRecord::new(mjd, 0.1 + 0.001 * i as f64, 0.2, 0.01, 0.001).unwrap();
397 let dx = 1.0 + 0.5 * i as f64;
398 let dy = -0.2 + 0.05 * i as f64;
399 record = record.with_cip_offsets(dx, dy).unwrap();
400 records.push(record);
401 }
402
403 let interpolator =
404 EopInterpolator::new(records).with_method(InterpolationMethod::Lagrange5);
405
406 let target_mjd = 59947.5;
407 let params = interpolator.get(target_mjd).unwrap();
408
409 let expected_dx = 1.0 + 0.5 * (target_mjd - 59945.0);
410 let expected_dy = -0.2 + 0.05 * (target_mjd - 59945.0);
411
412 assert!(params.dx.is_some());
413 assert!(params.dy.is_some());
414 assert!((params.dx.unwrap() - expected_dx).abs() < 1e-10);
415 assert!((params.dy.unwrap() - expected_dy).abs() < 1e-10);
416 assert!(params.flags.has_cip_offsets);
417 }
418
419 #[test]
420 fn test_out_of_range() {
421 let records = create_test_records();
422 let interpolator = EopInterpolator::new(records);
423
424 let result = interpolator.get(59944.0);
425 assert!(result.is_err());
426
427 let result = interpolator.get(59950.0);
428 assert!(result.is_err());
429 }
430
431 #[test]
432 fn test_max_gap_enforcement() {
433 let mut records = create_test_records();
434
435 records[3].mjd = 59955.0;
436 records[4].mjd = 59956.0;
437
438 let interpolator = EopInterpolator::new(records).with_max_gap(3.0);
439
440 let result = interpolator.get(59950.0);
441 assert!(result.is_err());
442 }
443
444 #[test]
445 fn test_time_span() {
446 let records = create_test_records();
447 let interpolator = EopInterpolator::new(records);
448
449 let (start, end) = interpolator.time_span().unwrap();
450 assert_eq!(start, 59945.0);
451 assert_eq!(end, 59949.0);
452 }
453
454 #[test]
455 fn test_empty_records() {
456 let interpolator = EopInterpolator::new(vec![]);
457 let result = interpolator.get(59945.0);
458 assert!(result.is_err());
459 let err = result.unwrap_err();
460 assert!(err.to_string().contains("No EOP records available"));
461 }
462
463 #[test]
464 fn test_empty_records_time_span() {
465 let interpolator = EopInterpolator::new(vec![]);
466 assert_eq!(interpolator.time_span(), None);
467 }
468
469 #[test]
470 fn test_lagrange_insufficient_points() {
471 let records = vec![
472 EopRecord::new(59945.0, 0.1, 0.2, 0.01, 0.001).unwrap(),
473 EopRecord::new(59946.0, 0.101, 0.202, 0.0101, 0.0011).unwrap(),
474 ];
475
476 let interpolator =
477 EopInterpolator::new(records).with_method(InterpolationMethod::Lagrange5);
478
479 let result = interpolator.get(59945.5);
480 assert!(result.is_err());
481 }
482
483 #[test]
484 fn test_record_count() {
485 let records = create_test_records();
486 let interpolator = EopInterpolator::new(records);
487 assert_eq!(interpolator.record_count(), 5);
488
489 let empty = EopInterpolator::new(vec![]);
490 assert_eq!(empty.record_count(), 0);
491 }
492
493 #[test]
494 fn test_find_center_index() {
495 let records = create_test_records();
496 let interpolator = EopInterpolator::new(records);
497
498 let center = interpolator.find_center_index(59947.0).unwrap();
499 assert_eq!(center, 2);
500
501 let edge = interpolator.find_center_index(59945.1).unwrap();
502 assert_eq!(edge, 0);
503 }
504
505 #[test]
506 fn test_lagrange_edge_cases() {
507 let mut records = Vec::new();
508 for i in 0..10 {
509 let mjd = 59945.0 + i as f64;
510 records.push(EopRecord::new(mjd, 0.1 + 0.001 * i as f64, 0.2, 0.01, 0.001).unwrap());
511 }
512
513 let interpolator =
514 EopInterpolator::new(records).with_method(InterpolationMethod::Lagrange5);
515
516 let near_start = interpolator.get(59946.0).unwrap();
517 assert!(near_start.x_p > 0.1);
518
519 let near_end = interpolator.get(59953.0).unwrap();
520 assert!(near_end.x_p > 0.1);
521 }
522}