deep_time/eop/
mod.rs

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

deep_time/eop/
mod.rs
