Skip to main content

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, Scale, 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/// Earth/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::Custom`].
55///
56/// Allows you to specify exactly which 0-based column contains each value
57/// when your input file does not match a standard IERS layout.
58#[derive(Debug, Clone, Copy)]
59pub struct CustomEopCols {
60    pub mjd: usize,
61    pub offset: usize,
62    pub pm_x: Option<usize>,
63    pub pm_y: Option<usize>,
64}
65
66/// A single parsed row of Earth Orientation Parameters.
67///
68/// - `mjd` — Modified Julian Date
69/// - `offset` — UT1 − UTC (or equivalent) in **seconds**
70/// - `pm_x`, `pm_y` — Polar motion in **arcseconds**
71#[derive(Debug, Clone, Copy)]
72pub struct EopDataRow {
73    pub mjd: Real,
74    /// e.g. UT1-UTC(s)
75    pub offset: Real,
76    /// polar motion x (arcsec)
77    pub pm_x: Real,
78    /// polar motion y (arcsec)
79    pub pm_y: Real,
80}
81
82/// Container for Earth/Body Orientation Parameters data.
83///
84/// - On Earth this would enable time scale conversions to and from
85///   the **UT1 time scale**.
86/// - Earth Orientation Parameters data is available from: <https://maia.usno.navy.mil/ser7/finals2000A.all>
87#[derive(Debug, Clone)]
88pub struct EopData {
89    rows: Vec<EopDataRow>,
90}
91
92#[cfg(feature = "std")]
93impl EopData {
94    /// Parse EOP data from any `std::io::BufRead` (file, network stream, etc.).
95    ///
96    /// Lines starting with `#` or longer than [`EopData::MAX_LINE_LEN`] are skipped.
97    /// The returned vector is always sorted by MJD.
98    pub fn data_from_reader<R: std::io::BufRead>(
99        mut reader: R,
100        format: EopFormat,
101        separator: Separator,
102    ) -> Result<Vec<EopDataRow>, DtErr> {
103        let mut line_buf = String::with_capacity(256);
104        let mut rows = Vec::new();
105
106        loop {
107            line_buf.clear();
108
109            let bytes_read = match reader.read_line(&mut line_buf) {
110                Ok(0) => break,
111                Ok(n) => n,
112                Err(e) => {
113                    return Err(an_err!(DtErrKind::IOErr, "{}", e));
114                }
115            };
116
117            if bytes_read > Self::MAX_LINE_LEN {
118                continue;
119            }
120
121            let trimmed = line_buf.trim();
122            if trimmed.is_empty() || trimmed.starts_with('#') {
123                continue;
124            }
125
126            if let Some(row) = Self::try_parse_row(trimmed, format, separator) {
127                rows.push(row);
128            }
129        }
130
131        if rows.is_empty() {
132            return Err(an_err!(DtErrKind::Empty));
133        }
134
135        rows.sort_by(|a, b| a.mjd.partial_cmp(&b.mjd).unwrap_or(Ordering::Equal));
136        Ok(rows)
137    }
138
139    /// Returns a [`Vec`] of [`EopDataRow`] from a text file on disk.
140    ///
141    /// ## Examples
142    ///
143    /// ```rust
144    /// # #[cfg(feature = "eop-tests")]
145    /// # {
146    /// use deep_time::eop::{EopData, EopFormat, Separator};
147    ///
148    /// let path = "finals.all.iau2000.txt";
149    /// let rows = EopData::data_from_text_file(path, EopFormat::Finals2000A, Separator::Whitespace).unwrap();
150    /// # }
151    /// ```
152    ///
153    /// ## See also
154    ///
155    /// - [`EopData::from_text_file`](../eop/struct.EopData.html#method.from_text_file)
156    pub fn data_from_text_file<P: AsRef<std::path::Path>>(
157        path: P,
158        format: EopFormat,
159        separator: Separator,
160    ) -> Result<Vec<EopDataRow>, DtErr> {
161        use std::fs::File;
162        use std::io::BufReader;
163
164        let path = path.as_ref();
165        let file = File::open(path).map_err(|e| an_err!(DtErrKind::IOErr, "{}", e))?;
166
167        let reader = BufReader::new(file);
168        Self::data_from_reader(reader, format, separator)
169    }
170
171    /// Create an [`EopData`] by loading from a text file on disk.
172    ///
173    /// ## Examples
174    ///
175    /// ```rust
176    /// # #[cfg(feature = "eop-tests")]
177    /// # {
178    /// use deep_time::eop::{EopData, EopFormat, Separator};
179    ///
180    /// let path = "finals.all.iau2000.txt";
181    /// let provider = EopData::from_text_file(path, EopFormat::Finals2000A, Separator::Whitespace).unwrap();
182    /// # }
183    /// ```
184    pub fn from_text_file<P: AsRef<std::path::Path>>(
185        path: P,
186        format: EopFormat,
187        separator: Separator,
188    ) -> Result<Self, DtErr> {
189        let rows = Self::data_from_text_file(path, format, separator)?;
190        Ok(Self { rows })
191    }
192}
193
194impl EopData {
195    pub const MAX_LINE_LEN: usize = 8192;
196
197    // Small helper — parses ONE row (shared by all paths)
198    fn try_parse_row(trimmed: &str, format: EopFormat, separator: Separator) -> Option<EopDataRow> {
199        let parts: Vec<&str> = match separator {
200            Separator::Whitespace => trimmed.split_whitespace().collect(),
201            Separator::Comma => trimmed.split(',').map(|s| s.trim()).collect(),
202            Separator::Tab => trimmed.split('\t').map(|s| s.trim()).collect(),
203            Separator::Pipe => trimmed.split('|').map(|s| s.trim()).collect(),
204            Separator::Semicolon => trimmed.split(';').map(|s| s.trim()).collect(),
205        };
206
207        if parts.len() < 2 {
208            return None;
209        }
210
211        let (mjd, offset, pm_x, pm_y) = match format {
212            EopFormat::Finals2000A => {
213                let mjd_idx = parts.iter().position(|p| {
214                    p.contains('.') && p.parse::<Real>().is_ok_and(|v| v > 30000.0)
215                })?;
216
217                let mut flag_count = 0;
218                let mut offset_value: Option<Real> = None;
219                let mut pm_x_val: Real = 0.0;
220                let mut pm_y_val: Real = 0.0;
221
222                for i in (mjd_idx + 1)..parts.len() {
223                    let token = parts[i];
224
225                    let is_flag = token == "I"
226                        || token == "P"
227                        || token.starts_with("I-")
228                        || token.starts_with("P-");
229
230                    if is_flag {
231                        flag_count += 1;
232
233                        if flag_count == 1 {
234                            // After first flag: x_p and y_p
235                            if let (Some(x_str), Some(y_str)) = (parts.get(i + 1), parts.get(i + 3))
236                            {
237                                pm_x_val = x_str.parse::<Real>().unwrap_or(0.0);
238                                pm_y_val = y_str.parse::<Real>().unwrap_or(0.0);
239                            }
240                        }
241
242                        if flag_count == 2 {
243                            let value_str = if token.starts_with("I-") || token.starts_with("P-") {
244                                &token[1..]
245                            } else if i + 1 < parts.len() {
246                                parts[i + 1]
247                            } else {
248                                break;
249                            };
250                            if let Ok(val) = value_str.parse::<Real>() {
251                                offset_value = Some(val);
252                            }
253                            break;
254                        }
255                    }
256                }
257
258                let offset = offset_value?;
259                let mjd_val = parts[mjd_idx].parse::<Real>().ok()?;
260
261                (mjd_val, offset, pm_x_val, pm_y_val)
262            }
263
264            EopFormat::C04 => {
265                let mjd = parts.get(4)?.parse::<Real>().ok()?;
266                let pm_x = parts
267                    .get(5)
268                    .unwrap_or(&"0.0")
269                    .parse::<Real>()
270                    .unwrap_or(0.0);
271                let pm_y = parts
272                    .get(6)
273                    .unwrap_or(&"0.0")
274                    .parse::<Real>()
275                    .unwrap_or(0.0);
276                let offset = parts.get(7)?.parse::<Real>().ok()?;
277                (mjd, offset, pm_x, pm_y)
278            }
279
280            EopFormat::Custom(cols) => {
281                let mjd = parts.get(cols.mjd)?.parse::<Real>().ok()?;
282                let offset = parts.get(cols.offset)?.parse::<Real>().ok()?;
283                let pm_x = if let Some(pm_x_col) = cols.pm_x {
284                    parts
285                        .get(pm_x_col)
286                        .unwrap_or(&"0.0")
287                        .parse::<Real>()
288                        .ok()
289                        .unwrap_or(0.0)
290                } else {
291                    0.0
292                };
293                let pm_y = if let Some(pm_y_col) = cols.pm_y {
294                    parts
295                        .get(pm_y_col)
296                        .unwrap_or(&"0.0")
297                        .parse::<Real>()
298                        .ok()
299                        .unwrap_or(0.0)
300                } else {
301                    0.0
302                };
303                (mjd, offset, pm_x, pm_y)
304            }
305        };
306
307        Some(EopDataRow {
308            mjd,
309            offset,
310            pm_x,
311            pm_y,
312        })
313    }
314
315    fn parse_lines<'a>(
316        lines: impl Iterator<Item = &'a str>,
317        format: EopFormat,
318        separator: Separator,
319    ) -> Result<Vec<EopDataRow>, DtErr> {
320        let mut rows = Vec::new();
321
322        for line in lines {
323            let trimmed = line.trim();
324            if trimmed.is_empty() || trimmed.starts_with('#') || trimmed.len() > Self::MAX_LINE_LEN
325            {
326                continue;
327            }
328
329            if let Some(row) = Self::try_parse_row(trimmed, format, separator) {
330                rows.push(row);
331            }
332        }
333
334        if rows.is_empty() {
335            return Err(an_err!(DtErrKind::Empty));
336        }
337
338        rows.sort_by(|a, b| a.mjd.partial_cmp(&b.mjd).unwrap_or(Ordering::Equal));
339        Ok(rows)
340    }
341
342    /// Parse EOP data from a `&str`.
343    ///
344    /// Useful when the data is already in memory (embedded resource,
345    /// downloaded string, etc.).
346    pub fn data_from_str(
347        s: &str,
348        format: EopFormat,
349        separator: Separator,
350    ) -> Result<Vec<EopDataRow>, DtErr> {
351        Self::parse_lines(s.lines(), format, separator)
352    }
353
354    /// Parse EOP data from raw bytes.
355    ///
356    /// The bytes are interpreted as UTF-8. Invalid UTF-8 sequences
357    /// result in an empty string (and therefore an error).
358    pub fn data_from_bytes(
359        bytes: &[u8],
360        format: EopFormat,
361        separator: Separator,
362    ) -> Result<Vec<EopDataRow>, DtErr> {
363        let s = core::str::from_utf8(bytes).unwrap_or("");
364        Self::data_from_str(s, format, separator)
365    }
366
367    /// Create an [`EopData`] from a string slice.
368    pub fn from_str(s: &str, format: EopFormat, separator: Separator) -> Result<Self, DtErr> {
369        let rows = Self::data_from_str(s, format, separator)?;
370        Ok(Self { rows })
371    }
372
373    /// Create an [`EopData`] from raw bytes.
374    pub fn from_bytes(
375        bytes: &[u8],
376        format: EopFormat,
377        separator: Separator,
378    ) -> Result<Self, DtErr> {
379        let rows = Self::data_from_bytes(bytes, format, separator)?;
380        Ok(Self { rows })
381    }
382
383    /// Returns all interpolated orientation parameters (offset + polar motion)
384    /// at the given MJD.
385    ///
386    /// Returns `None` if the table is empty or the MJD is completely outside
387    /// the loaded data.
388    pub fn eop_offset(&self, mjd: Real) -> Option<EopOffset> {
389        if self.rows.is_empty() {
390            return None;
391        }
392
393        let idx = match self
394            .rows
395            .binary_search_by(|probe| probe.mjd.partial_cmp(&mjd).unwrap_or(Ordering::Equal))
396        {
397            Ok(i) => i,
398            Err(i) => {
399                if i == 0 {
400                    let row = &self.rows[0];
401                    return Some(EopOffset {
402                        offset: row.offset,
403                        pm_x: row.pm_x,
404                        pm_y: row.pm_y,
405                    });
406                }
407                if i >= self.rows.len() {
408                    let row = &self.rows[self.rows.len() - 1];
409                    return Some(EopOffset {
410                        offset: row.offset,
411                        pm_x: row.pm_x,
412                        pm_y: row.pm_y,
413                    });
414                }
415                i - 1
416            }
417        };
418
419        if idx + 1 < self.rows.len() {
420            let e0 = &self.rows[idx];
421            let e1 = &self.rows[idx + 1];
422
423            let t = (mjd - e0.mjd) / (e1.mjd - e0.mjd);
424
425            let offset = e0.offset + t * (e1.offset - e0.offset);
426            let pm_x = e0.pm_x + t * (e1.pm_x - e0.pm_x);
427            let pm_y = e0.pm_y + t * (e1.pm_y - e0.pm_y);
428
429            Some(EopOffset { offset, pm_x, pm_y })
430        } else {
431            let row = &self.rows[idx];
432            Some(EopOffset {
433                offset: row.offset,
434                pm_x: row.pm_x,
435                pm_y: row.pm_y,
436            })
437        }
438    }
439}
440
441/// Interpolated Body Orientation Parameters at a specific MJD.
442///
443/// Contains everything needed for high-precision sidereal time
444/// and polar-motion corrections.
445#[derive(Debug, Clone, Copy, Default)]
446pub struct EopOffset {
447    /// Value in **seconds** e.g. UT1 − UTC offset
448    pub offset: Real,
449    /// Polar motion x-coordinate in **arcseconds**
450    pub pm_x: Real,
451    /// Polar motion y-coordinate in **arcseconds**
452    pub pm_y: Real,
453}
454
455impl Dt {
456    /// Get an orientation parameters offset in seconds inside a struct: ([`EopOffset`])
457    /// for a particular Modified Julian Date.
458    ///
459    /// - On Earth this would be the UT1 time scale.
460    /// - Earth Orientation Parameters data is available from: <https://maia.usno.navy.mil/ser7/finals2000A.all>
461    pub fn mjd_to_eop_offset(mjd: Real, op_data: &EopData) -> Result<EopOffset, DtErr> {
462        let offset = op_data
463            .eop_offset(mjd)
464            .ok_or_else(|| an_err!(DtErrKind::MjdOutOfRange, "{mjd}"))?;
465        Ok(offset)
466    }
467
468    /// Get an orientation parameters offset in seconds for a particular Modified Julian Date.
469    ///
470    /// - On Earth this would be the UT1 time scale.
471    /// - Earth Orientation Parameters data is available from: <https://maia.usno.navy.mil/ser7/finals2000A.all>
472    #[inline]
473    pub fn mjd_to_eop_offset_f(mjd: Real, op_data: &EopData) -> Result<Real, DtErr> {
474        Self::mjd_to_eop_offset(mjd, op_data).map(|res| res.offset)
475    }
476
477    /// Offsets a [`Dt`] using orientation parameters data.
478    ///
479    /// - On Earth this would be the UT1 time scale.
480    /// - Earth Orientation Parameters data is available from: <https://maia.usno.navy.mil/ser7/finals2000A.all>
481    #[inline]
482    pub fn to_eop(&self, op_data: &EopData) -> Result<Self, DtErr> {
483        Ok(self.add(Dt::from_sec_f(
484            Self::mjd_to_eop_offset_f(self.to_mjd_f_raw(), op_data)?,
485            Scale::TAI,
486        )))
487    }
488
489    /// Convert a [`Dt`] already offset using orientation parameters data back to whatever
490    /// it was before.
491    ///
492    /// - On Earth this would be the UT1 time scale.
493    /// - Earth Orientation Parameters data is available from: <https://maia.usno.navy.mil/ser7/finals2000A.all>
494    pub fn from_eop(&self, op_data: &EopData) -> Result<Self, DtErr> {
495        if op_data.rows.is_empty() {
496            return Err(an_err!(DtErrKind::Empty));
497        }
498        let mut guess = *self;
499
500        for _ in 0..8 {
501            let mjd = guess.to_mjd_f_raw();
502            let offset = op_data
503                .eop_offset(mjd)
504                .ok_or_else(|| an_err!(DtErrKind::MjdOutOfRange, "{mjd}"))?
505                .offset;
506
507            guess = self.sub(Dt::from_sec_f(offset, Scale::TAI)); // TODO: guess or self?
508        }
509
510        Ok(guess)
511    }
512}