deep_time/eop/
mod.rs

1//! Earth Orientation Parameters (EOP) data parser and interpolator.
2//!
3//! Loads UT1–UTC offsets and polar motion (x, y) from standard IERS formats
4//! (`Finals2000A`, `C04`) or custom column layouts.
5//!
6//! Provides interpolation at any MJD and integrates with [`Dt`]
7//! for UT1 time scale conversions.
8//!
9//! Tries to be planet/body agnostic, such that custom data for any world
10//! might be able to be loaded and used.
11
12#![allow(clippy::indexing_slicing)]
13#![allow(clippy::excessive_precision)]
14#![allow(clippy::approx_constant)]
15#![allow(clippy::eq_op)]
16
17use crate::{Dt, DtErr, DtErrKind, Real, an_err};
18use alloc::string::String;
19use alloc::vec::Vec;
20use core::cmp::Ordering;
21
22#[derive(Debug, Clone, Copy, Default)]
23pub enum Separator {
24    #[default]
25    Whitespace,
26    Comma,
27    Tab,
28    Pipe,
29    Semicolon,
30}
31
32/// Body Orientation Parameters Format.
33///
34/// Formats to provide to the parser, including a
35/// custom one to allow specific column indices.
36///
37/// - `Finals2000A` such as is available from
38///   https://maia.usno.navy.mil/ser7/finals2000A.all
39/// - `C04` such as is available from
40///   https://datacenter.iers.org/data/latestVersion/EOP_20u24_C04_one_file_1962-now.txt
41/// - `Custom` so you can provide your own specific column indices
42///   using [`CustomEopCols`].
43#[derive(Debug, Clone, Copy, Default)]
44pub enum EopFormat {
45    /// finals2000A.all / finals.all.iau2000.txt style files
46    #[default]
47    Finals2000A,
48    /// C04 long-term series
49    C04,
50    /// User-defined column indices (0-based)
51    Custom(CustomEopCols),
52}
53
54/// For use with [`EopFormat`].
55#[derive(Debug, Clone, Copy)]
56pub struct CustomEopCols {
57    pub mjd: usize,
58    pub offset: usize,
59    pub pm_x: Option<usize>,
60    pub pm_y: Option<usize>,
61}
62
63#[derive(Debug, Clone, Copy)]
64pub struct EopDataRow {
65    pub mjd: Real,
66    /// e.g. UT1-UTC(s)
67    pub offset: Real,
68    /// polar motion x (arcsec)
69    pub pm_x: Real,
70    /// polar motion y (arcsec)
71    pub pm_y: Real,
72}
73
74/// Container for Body Orientation Parameters data.
75///
76/// - On Earth this would enable time scale conversions to and from
77///   the **UT1 time scale**.
78/// - Earth Orientation Parameters data is available from: https://maia.usno.navy.mil/ser7/finals2000A.all
79#[derive(Debug, Clone)]
80pub struct EopData {
81    rows: Vec<EopDataRow>,
82}
83
84#[cfg(feature = "std")]
85impl EopData {
86    pub fn data_from_reader<R: std::io::BufRead>(
87        mut reader: R,
88        format: EopFormat,
89        separator: Separator,
90    ) -> Result<Vec<EopDataRow>, DtErr> {
91        let mut line_buf = String::with_capacity(256);
92        let mut rows = Vec::new();
93
94        loop {
95            line_buf.clear();
96
97            let bytes_read = match reader.read_line(&mut line_buf) {
98                Ok(0) => break,
99                Ok(n) => n,
100                Err(e) => {
101                    return Err(an_err!(DtErrKind::IOErr, "read line: {}", e));
102                }
103            };
104
105            if bytes_read > Self::MAX_LINE_LEN {
106                continue;
107            }
108
109            let trimmed = line_buf.trim();
110            if trimmed.is_empty() || trimmed.starts_with('#') {
111                continue;
112            }
113
114            if let Some(row) = Self::try_parse_row(trimmed, format, separator) {
115                rows.push(row);
116            }
117        }
118
119        if rows.is_empty() {
120            return Err(an_err!(DtErrKind::Incomplete, "no valid rows"));
121        }
122
123        rows.sort_by(|a, b| a.mjd.partial_cmp(&b.mjd).unwrap_or(Ordering::Equal));
124        Ok(rows)
125    }
126
127    pub fn data_from_text_file<P: AsRef<std::path::Path>>(
128        path: P,
129        format: EopFormat,
130        separator: Separator,
131    ) -> Result<Vec<EopDataRow>, DtErr> {
132        use std::fs::File;
133        use std::io::BufReader;
134
135        let path = path.as_ref();
136        let file = File::open(path)
137            .map_err(|e| an_err!(DtErrKind::IOErr, "open file: '{}': {}", path.display(), e))?;
138
139        let reader = BufReader::new(file);
140        Self::data_from_reader(reader, format, separator)
141    }
142
143    pub fn from_text_file<P: AsRef<std::path::Path>>(
144        path: P,
145        format: EopFormat,
146        separator: Separator,
147    ) -> Result<Self, DtErr> {
148        let rows = Self::data_from_text_file(path, format, separator)?;
149        Ok(Self { rows })
150    }
151}
152
153impl EopData {
154    pub const MAX_LINE_LEN: usize = 8192;
155
156    // Small helper — parses ONE row (shared by all paths)
157    fn try_parse_row(trimmed: &str, format: EopFormat, separator: Separator) -> Option<EopDataRow> {
158        let parts: Vec<&str> = match separator {
159            Separator::Whitespace => trimmed.split_whitespace().collect(),
160            Separator::Comma => trimmed.split(',').map(|s| s.trim()).collect(),
161            Separator::Tab => trimmed.split('\t').map(|s| s.trim()).collect(),
162            Separator::Pipe => trimmed.split('|').map(|s| s.trim()).collect(),
163            Separator::Semicolon => trimmed.split(';').map(|s| s.trim()).collect(),
164        };
165
166        if parts.len() < 2 {
167            return None;
168        }
169
170        let (mjd, offset, pm_x, pm_y) = match format {
171            EopFormat::Finals2000A => {
172                let mjd_idx = parts.iter().position(|p| {
173                    p.contains('.') && p.parse::<Real>().is_ok_and(|v| v > 30000.0)
174                })?;
175
176                let mut flag_count = 0;
177                let mut offset_value: Option<Real> = None;
178                let mut pm_x_val: Real = 0.0;
179                let mut pm_y_val: Real = 0.0;
180
181                for i in (mjd_idx + 1)..parts.len() {
182                    let token = parts[i];
183
184                    let is_flag = token == "I"
185                        || token == "P"
186                        || token.starts_with("I-")
187                        || token.starts_with("P-");
188
189                    if is_flag {
190                        flag_count += 1;
191
192                        if flag_count == 1 {
193                            // After first flag: x_p and y_p
194                            if let (Some(x_str), Some(y_str)) = (parts.get(i + 1), parts.get(i + 3))
195                            {
196                                pm_x_val = x_str.parse::<Real>().unwrap_or(0.0);
197                                pm_y_val = y_str.parse::<Real>().unwrap_or(0.0);
198                            }
199                        }
200
201                        if flag_count == 2 {
202                            let value_str = if token.starts_with("I-") || token.starts_with("P-") {
203                                &token[1..]
204                            } else if i + 1 < parts.len() {
205                                parts[i + 1]
206                            } else {
207                                break;
208                            };
209                            if let Ok(val) = value_str.parse::<Real>() {
210                                offset_value = Some(val);
211                            }
212                            break;
213                        }
214                    }
215                }
216
217                let offset = offset_value?;
218                let mjd_val = parts[mjd_idx].parse::<Real>().ok()?;
219
220                (mjd_val, offset, pm_x_val, pm_y_val)
221            }
222
223            EopFormat::C04 => {
224                let mjd = parts.get(4)?.parse::<Real>().ok()?;
225                let pm_x = parts
226                    .get(5)
227                    .unwrap_or(&"0.0")
228                    .parse::<Real>()
229                    .unwrap_or(0.0);
230                let pm_y = parts
231                    .get(6)
232                    .unwrap_or(&"0.0")
233                    .parse::<Real>()
234                    .unwrap_or(0.0);
235                let offset = parts.get(7)?.parse::<Real>().ok()?;
236                (mjd, offset, pm_x, pm_y)
237            }
238
239            EopFormat::Custom(cols) => {
240                let mjd = parts.get(cols.mjd)?.parse::<Real>().ok()?;
241                let offset = parts.get(cols.offset)?.parse::<Real>().ok()?;
242                let pm_x = if let Some(pm_x_col) = cols.pm_x {
243                    parts
244                        .get(pm_x_col)
245                        .unwrap_or(&"0.0")
246                        .parse::<Real>()
247                        .ok()
248                        .unwrap_or(0.0)
249                } else {
250                    0.0
251                };
252                let pm_y = if let Some(pm_y_col) = cols.pm_y {
253                    parts
254                        .get(pm_y_col)
255                        .unwrap_or(&"0.0")
256                        .parse::<Real>()
257                        .ok()
258                        .unwrap_or(0.0)
259                } else {
260                    0.0
261                };
262                (mjd, offset, pm_x, pm_y)
263            }
264        };
265
266        Some(EopDataRow {
267            mjd,
268            offset,
269            pm_x,
270            pm_y,
271        })
272    }
273
274    fn parse_lines<'a>(
275        lines: impl Iterator<Item = &'a str>,
276        format: EopFormat,
277        separator: Separator,
278    ) -> Result<Vec<EopDataRow>, DtErr> {
279        let mut rows = Vec::new();
280
281        for line in lines {
282            let trimmed = line.trim();
283            if trimmed.is_empty() || trimmed.starts_with('#') || trimmed.len() > Self::MAX_LINE_LEN
284            {
285                continue;
286            }
287
288            if let Some(row) = Self::try_parse_row(trimmed, format, separator) {
289                rows.push(row);
290            }
291        }
292
293        if rows.is_empty() {
294            return Err(an_err!(DtErrKind::Incomplete, "no valid rows"));
295        }
296
297        rows.sort_by(|a, b| a.mjd.partial_cmp(&b.mjd).unwrap_or(Ordering::Equal));
298        Ok(rows)
299    }
300
301    pub fn data_from_str(
302        s: &str,
303        format: EopFormat,
304        separator: Separator,
305    ) -> Result<Vec<EopDataRow>, DtErr> {
306        Self::parse_lines(s.lines(), format, separator)
307    }
308
309    pub fn data_from_bytes(
310        bytes: &[u8],
311        format: EopFormat,
312        separator: Separator,
313    ) -> Result<Vec<EopDataRow>, DtErr> {
314        let s = core::str::from_utf8(bytes).unwrap_or("");
315        Self::data_from_str(s, format, separator)
316    }
317
318    pub fn from_str(s: &str, format: EopFormat, separator: Separator) -> Result<Self, DtErr> {
319        let rows = Self::data_from_str(s, format, separator)?;
320        Ok(Self { rows })
321    }
322
323    pub fn from_bytes(
324        bytes: &[u8],
325        format: EopFormat,
326        separator: Separator,
327    ) -> Result<Self, DtErr> {
328        let rows = Self::data_from_bytes(bytes, format, separator)?;
329        Ok(Self { rows })
330    }
331
332    /// Returns all interpolated orientation parameters (offset + polar motion)
333    /// at the given MJD.
334    ///
335    /// Returns `None` if the table is empty or the MJD is completely outside
336    /// the loaded data.
337    pub fn eop_offset(&self, mjd: Real) -> Option<EopOffset> {
338        if self.rows.is_empty() {
339            return None;
340        }
341
342        let idx = match self
343            .rows
344            .binary_search_by(|probe| probe.mjd.partial_cmp(&mjd).unwrap_or(Ordering::Equal))
345        {
346            Ok(i) => i,
347            Err(i) => {
348                if i == 0 {
349                    let row = &self.rows[0];
350                    return Some(EopOffset {
351                        offset: row.offset,
352                        pm_x: row.pm_x,
353                        pm_y: row.pm_y,
354                    });
355                }
356                if i >= self.rows.len() {
357                    let row = &self.rows[self.rows.len() - 1];
358                    return Some(EopOffset {
359                        offset: row.offset,
360                        pm_x: row.pm_x,
361                        pm_y: row.pm_y,
362                    });
363                }
364                i - 1
365            }
366        };
367
368        if idx + 1 < self.rows.len() {
369            let e0 = &self.rows[idx];
370            let e1 = &self.rows[idx + 1];
371
372            let t = (mjd - e0.mjd) / (e1.mjd - e0.mjd);
373
374            let offset = e0.offset + t * (e1.offset - e0.offset);
375            let pm_x = e0.pm_x + t * (e1.pm_x - e0.pm_x);
376            let pm_y = e0.pm_y + t * (e1.pm_y - e0.pm_y);
377
378            Some(EopOffset { offset, pm_x, pm_y })
379        } else {
380            let row = &self.rows[idx];
381            Some(EopOffset {
382                offset: row.offset,
383                pm_x: row.pm_x,
384                pm_y: row.pm_y,
385            })
386        }
387    }
388}
389
390/// Interpolated Body Orientation Parameters at a specific MJD.
391///
392/// Contains everything needed for high-precision sidereal time
393/// and polar-motion corrections.
394#[derive(Debug, Clone, Copy, Default)]
395pub struct EopOffset {
396    /// Value in **seconds** e.g. UT1 − UTC offset
397    pub offset: Real,
398    /// Polar motion x-coordinate in **arcseconds**
399    pub pm_x: Real,
400    /// Polar motion y-coordinate in **arcseconds**
401    pub pm_y: Real,
402}
403
404impl Dt {
405    /// Get an orientation parameters offset in seconds inside a struct: ([`EopOffset`])
406    /// for a particular Modified Julian Date.
407    ///
408    /// - On Earth this would be the UT1 time scale.
409    /// - Earth Orientation Parameters data is available from: https://maia.usno.navy.mil/ser7/finals2000A.all
410    pub fn mjd_to_eop_offset(mjd: Real, op_data: &EopData) -> Result<EopOffset, DtErr> {
411        let offset = op_data
412            .eop_offset(mjd)
413            .ok_or_else(|| an_err!(DtErrKind::OutOfRange, "mjd: {mjd}"))?;
414        Ok(offset)
415    }
416
417    /// Get an orientation parameters offset in seconds for a particular Modified Julian Date.
418    ///
419    /// - On Earth this would be the UT1 time scale.
420    /// - Earth Orientation Parameters data is available from: https://maia.usno.navy.mil/ser7/finals2000A.all
421    #[inline]
422    pub fn mjd_to_eop_offset_f(mjd: Real, op_data: &EopData) -> Result<Real, DtErr> {
423        Self::mjd_to_eop_offset(mjd, op_data).map(|res| res.offset)
424    }
425
426    /// Offsets a [`Dt`] using orientation parameters data.
427    ///
428    /// - On Earth this would be the UT1 time scale.
429    /// - Earth Orientation Parameters data is available from: https://maia.usno.navy.mil/ser7/finals2000A.all
430    #[inline]
431    pub fn to_eop(&self, op_data: &EopData) -> Result<Self, DtErr> {
432        Ok(self.add(Dt::from_sec_f(Self::mjd_to_eop_offset_f(
433            self.to_mjd_f(),
434            op_data,
435        )?)))
436    }
437
438    /// Convert a [`Dt`] already offset using orientation parameters data back to whatever
439    /// it was before.
440    ///
441    /// - On Earth this would be the UT1 time scale.
442    /// - Earth Orientation Parameters data is available from: https://maia.usno.navy.mil/ser7/finals2000A.all
443    pub fn from_eop(&self, op_data: &EopData) -> Result<Self, DtErr> {
444        if op_data.rows.is_empty() {
445            return Err(an_err!(DtErrKind::InternalErr, "contains no data"));
446        }
447        let mut guess = *self;
448
449        for _ in 0..8 {
450            let mjd = guess.to_mjd_f();
451            let offset = op_data
452                .eop_offset(mjd)
453                .ok_or_else(|| an_err!(DtErrKind::OutOfRange, "mjd: {mjd}"))?
454                .offset;
455
456            guess = self.sub(Dt::from_sec_f(offset)); // TODO: guess or self?
457        }
458
459        Ok(guess)
460    }
461}
deep_time/eop/mod.rs

deep_time/eop/
mod.rs
