lox_earth/
eop.rs

1// SPDX-FileCopyrightText: 2025 Helge Eichhorn <git@helgeeichhorn.de>
2//
3// SPDX-License-Identifier: MPL-2.0
4
5use std::{
6    path::{Path, PathBuf},
7    sync::Arc,
8};
9
10use csv::{ByteRecord, ReaderBuilder};
11use lox_core::f64::consts::{SECONDS_BETWEEN_MJD_AND_J2000, SECONDS_PER_DAY};
12use lox_core::units::Angle;
13use lox_frames::iers::{Corrections, polar_motion::PoleCoords};
14use lox_io::spice::lsk::LeapSecondsKernel;
15use lox_math::series::{InterpolationType, Series, SeriesError};
16use lox_time::{
17    deltas::TimeDelta,
18    julian_dates::JulianDate,
19    utc::{
20        Utc, UtcError,
21        leap_seconds::{DefaultLeapSecondsProvider, LeapSecondsProvider},
22        transformations::ToUtc,
23    },
24};
25use serde::Deserialize;
26use thiserror::Error;
27
28#[derive(Debug, Deserialize)]
29struct EopRecord {
30    #[serde(rename = "MJD")]
31    modified_julian_date: f64,
32    #[serde(rename = "Year")]
33    year: i64,
34    #[serde(rename = "Month")]
35    month: u8,
36    #[serde(rename = "Day")]
37    day: u8,
38    x_pole: Option<f64>,
39    y_pole: Option<f64>,
40    #[serde(rename = "UT1-UTC")]
41    delta_ut1_utc: Option<f64>,
42    #[serde(rename = "dPsi")]
43    dpsi: Option<f64>,
44    #[serde(rename = "dEpsilon")]
45    deps: Option<f64>,
46    #[serde(rename = "dX")]
47    dx: Option<f64>,
48    #[serde(rename = "dY")]
49    dy: Option<f64>,
50}
51
52impl EopRecord {
53    fn merge(mut self, other: &Self) -> Self {
54        self.dpsi = self.dpsi.or(other.dpsi);
55        self.deps = self.deps.or(other.deps);
56        self.dx = self.dx.or(other.dx);
57        self.dy = self.dy.or(other.dy);
58        self
59    }
60}
61
62/// Errors that can occur when parsing EOP data files.
63#[derive(Debug, Error)]
64pub enum EopParserError {
65    /// CSV parsing error.
66    #[error("{0}")]
67    Csv(String),
68    /// No input files were provided.
69    #[error("either a 'finals.all.csv' or a 'finals2000A.all.csv' file needs to be provided")]
70    NoFiles,
71    /// Column lengths do not match.
72    #[error("mismatched dimensions for columns '{0} (n={1})' and '{2} (n={3})'")]
73    DimensionMismatch(String, usize, String, usize),
74    /// UTC conversion error.
75    #[error(transparent)]
76    Utc(#[from] UtcError),
77    /// Interpolation series error.
78    #[error(transparent)]
79    Series(#[from] SeriesError),
80}
81
82impl From<csv::Error> for EopParserError {
83    fn from(err: csv::Error) -> Self {
84        EopParserError::Csv(err.to_string())
85    }
86}
87
88/// Builder for parsing IERS EOP CSV files into an [`EopProvider`].
89#[derive(Default)]
90pub struct EopParser {
91    paths: (Option<PathBuf>, Option<PathBuf>),
92    lsk: Option<LeapSecondsKernel>,
93}
94
95impl EopParser {
96    /// Creates a new parser with no files configured.
97    pub fn new() -> Self {
98        Self::default()
99    }
100
101    /// Sets a single EOP CSV file path.
102    pub fn from_path(mut self, path: impl AsRef<Path>) -> Self {
103        self.paths = (Some(path.as_ref().to_owned()), None);
104        self
105    }
106
107    /// Sets two EOP CSV file paths (IAU 1980 and IAU 2000).
108    pub fn from_paths(mut self, path1: impl AsRef<Path>, path2: impl AsRef<Path>) -> Self {
109        self.paths = (
110            Some(path1.as_ref().to_owned()),
111            Some(path2.as_ref().to_owned()),
112        );
113        self
114    }
115
116    /// Configures a leap seconds kernel for UTC conversion.
117    pub fn with_leap_seconds_kernel(mut self, lsp: LeapSecondsKernel) -> Self {
118        self.lsk = Some(lsp);
119        self
120    }
121
122    /// Parses the configured files and returns an [`EopProvider`].
123    pub fn parse(self) -> Result<EopProvider, EopParserError> {
124        let n = if let Some(iau1980) = self.paths.0.as_ref() {
125            let mut reader = ReaderBuilder::new().delimiter(b';').from_path(iau1980)?;
126            reader.records().count()
127        } else {
128            return Err(EopParserError::NoFiles);
129        };
130
131        let mut j2000: Vec<f64> = Vec::with_capacity(n);
132        let mut delta_ut1_tai: Vec<f64> = Vec::with_capacity(n);
133        let mut x_pole: Vec<f64> = Vec::with_capacity(n);
134        let mut y_pole: Vec<f64> = Vec::with_capacity(n);
135        let mut dpsi: Vec<f64> = Vec::with_capacity(n);
136        let mut deps: Vec<f64> = Vec::with_capacity(n);
137        let mut dx: Vec<f64> = Vec::with_capacity(n);
138        let mut dy: Vec<f64> = Vec::with_capacity(n);
139
140        // The EOP from the IERS are split into two different files based on the two IAU conventions for some reason.
141        // Both files have the same header and identical data except that the columns for the IAU1980 parameters are
142        // only populated in `finals.all.csv` and the columns for the IAU2000 parameters are only populated in
143        // `finals2000A.all.csv`.
144        // We do not want to force the user to provide both files if they do not need them. At the same time, we want
145        // to parse both files in one pass if they are present to avoid looping over tens of thousands of lines
146        // multiple time.
147        // I have not been able to fulfil both requirements and make the compiler happy without pretending that there
148        // are always two files. In case the user provided two files, we parse both files and merge the records.
149        // If we have only one file, we still parse it twice and record merging is a no-op.
150        let mut rdr1 = ReaderBuilder::new()
151            .delimiter(b';')
152            .from_path(self.paths.0.as_ref().unwrap())?;
153        let mut rdr2 = ReaderBuilder::new()
154            .delimiter(b';')
155            .from_path(self.paths.1.as_ref().or(self.paths.0.as_ref()).unwrap())?;
156
157        let mut raw1 = ByteRecord::new();
158        let mut raw2 = ByteRecord::new();
159        let headers = rdr1.byte_headers()?.clone();
160
161        while rdr1.read_byte_record(&mut raw1)? && rdr2.read_byte_record(&mut raw2)? {
162            let r1: EopRecord = raw1.deserialize(Some(&headers))?;
163            let r2: EopRecord = raw2.deserialize(Some(&headers))?;
164            let r = r1.merge(&r2);
165
166            j2000.push(r.modified_julian_date * SECONDS_PER_DAY - SECONDS_BETWEEN_MJD_AND_J2000);
167
168            if let Some(delta_ut1_utc) = r.delta_ut1_utc {
169                let utc = Utc::builder().with_ymd(r.year, r.month, r.day).build()?;
170                let delta_tai_utc = self
171                    .lsk
172                    .as_ref()
173                    .map(|lsp| lsp.delta_utc_tai(utc))
174                    .unwrap_or_else(|| DefaultLeapSecondsProvider.delta_utc_tai(utc));
175                delta_ut1_tai.push(delta_ut1_utc + delta_tai_utc.to_seconds().to_f64())
176            }
177
178            if let (Some(xp), Some(yp)) = (r.x_pole, r.y_pole) {
179                x_pole.push(xp);
180                y_pole.push(yp);
181            }
182
183            if let (Some(d_psi), Some(d_eps)) = (r.dpsi, r.deps) {
184                dpsi.push(d_psi);
185                deps.push(d_eps);
186            }
187
188            if let (Some(d_x), Some(d_y)) = (r.dx, r.dy) {
189                dx.push(d_x);
190                dy.push(d_y);
191            }
192        }
193
194        let n = x_pole.len();
195        let npn = dpsi.len().max(dx.len());
196
197        let index: Arc<[f64]> = Arc::from(&j2000[0..n]);
198        let np_index: Arc<[f64]> = Arc::from(&j2000[0..npn]);
199
200        Ok(EopProvider {
201            polar_motion: (
202                Series::try_new(index.clone(), x_pole, InterpolationType::CubicSpline)?,
203                Series::try_new(index.clone(), y_pole, InterpolationType::CubicSpline)?,
204            ),
205            delta_ut1_tai: Series::try_new(
206                index.clone(),
207                delta_ut1_tai,
208                InterpolationType::CubicSpline,
209            )?,
210            nut_prec: NutPrecCorrections {
211                iau1980: if !dpsi.is_empty() {
212                    Some((
213                        Series::try_new(np_index.clone(), dpsi, InterpolationType::CubicSpline)?,
214                        Series::try_new(np_index.clone(), deps, InterpolationType::CubicSpline)?,
215                    ))
216                } else {
217                    None
218                },
219                iau2000: if !dx.is_empty() {
220                    Some((
221                        Series::try_new(np_index.clone(), dx, InterpolationType::CubicSpline)?,
222                        Series::try_new(np_index.clone(), dy, InterpolationType::CubicSpline)?,
223                    ))
224                } else {
225                    None
226                },
227            },
228            lsk: self.lsk.clone(),
229        })
230    }
231}
232
233/// Errors returned by [`EopProvider`] queries.
234#[derive(Debug, Error)]
235pub enum EopProviderError {
236    /// Time-scale offset conversion failed.
237    #[error("offset error: {0}")]
238    Offset(String),
239    /// The queried epoch lies outside the data range (single value).
240    #[error("value was extrapolated")]
241    ExtrapolatedValue(f64),
242    /// The queried epoch lies outside the data range (pair of values).
243    #[error("values were extrapolated")]
244    ExtrapolatedValues(f64, f64),
245    /// No IAU 1980 nutation/precession data was loaded.
246    #[error("no 'finals.all.csv' file was loaded")]
247    MissingIau1980,
248    /// No IAU 2000 nutation/precession data was loaded.
249    #[error("no 'finals2000A.all.csv' file was loaded")]
250    MissingIau2000,
251}
252
253#[derive(Debug, Clone)]
254struct NutPrecCorrections {
255    iau1980: Option<(Series, Series)>,
256    iau2000: Option<(Series, Series)>,
257}
258
259/// Interpolates Earth orientation parameters from parsed IERS data.
260#[derive(Debug)]
261pub struct EopProvider {
262    polar_motion: (Series, Series),
263    delta_ut1_tai: Series,
264    nut_prec: NutPrecCorrections,
265    lsk: Option<LeapSecondsKernel>,
266}
267
268impl EopProvider {
269    /// Returns the pole coordinates at the given epoch.
270    pub fn polar_motion<T: ToUtc>(&self, t: T) -> Result<PoleCoords, EopProviderError> {
271        let t = t.to_utc().seconds_since_j2000();
272        let xp = self.polar_motion.0.interpolate(t);
273        let yp = self.polar_motion.1.interpolate(t);
274        let (min, _) = self.polar_motion.0.first();
275        let (max, _) = self.polar_motion.0.last();
276        if t < min || t > max {
277            return Err(EopProviderError::ExtrapolatedValues(xp, yp));
278        }
279        Ok(PoleCoords {
280            xp: Angle::arcseconds(xp),
281            yp: Angle::arcseconds(yp),
282        })
283    }
284
285    /// Returns the IAU 1980 nutation/precession corrections at the given epoch.
286    pub fn nutation_precession_iau1980<T: ToUtc>(
287        &self,
288        t: T,
289    ) -> Result<Corrections, EopProviderError> {
290        let Some(nut_prec) = &self.nut_prec.iau1980 else {
291            return Err(EopProviderError::MissingIau1980);
292        };
293        let t = t.to_utc().seconds_since_j2000();
294        let dpsi = nut_prec.0.interpolate(t);
295        let deps = nut_prec.1.interpolate(t);
296        let (min, _) = nut_prec.0.first();
297        let (max, _) = nut_prec.0.last();
298        if t < min || t > max {
299            return Err(EopProviderError::ExtrapolatedValues(dpsi, deps));
300        }
301        Ok(Corrections(
302            Angle::arcseconds(dpsi * 1e-3),
303            Angle::arcseconds(deps * 1e-3),
304        ))
305    }
306
307    /// Returns the IAU 2000 nutation/precession corrections at the given epoch.
308    pub fn nutation_precession_iau2000<T: ToUtc>(
309        &self,
310        t: T,
311    ) -> Result<Corrections, EopProviderError> {
312        let Some(nut_prec) = &self.nut_prec.iau2000 else {
313            return Err(EopProviderError::MissingIau2000);
314        };
315        let t = t.to_utc().seconds_since_j2000();
316        let dx = nut_prec.0.interpolate(t);
317        let dy = nut_prec.1.interpolate(t);
318        let (min, _) = nut_prec.0.first();
319        let (max, _) = nut_prec.0.last();
320        if t < min || t > max {
321            return Err(EopProviderError::ExtrapolatedValues(dx, dy));
322        }
323        Ok(Corrections(
324            Angle::arcseconds(dx * 1e-3),
325            Angle::arcseconds(dy * 1e-3),
326        ))
327    }
328
329    /// Returns the UT1-TAI offset for the given TAI epoch.
330    pub fn delta_ut1_tai(&self, tai: TimeDelta) -> Result<TimeDelta, EopProviderError> {
331        let seconds = tai.seconds_since_j2000();
332        let (t0, _) = self.delta_ut1_tai.first();
333        let (tn, _) = self.delta_ut1_tai.last();
334        let val = self.delta_ut1_tai.interpolate(seconds);
335        if seconds < t0 || seconds > tn {
336            return Err(EopProviderError::ExtrapolatedValue(val));
337        }
338        Ok(TimeDelta::from_seconds_f64(val))
339    }
340
341    /// Returns the TAI-UT1 offset for the given UT1 epoch.
342    pub fn delta_tai_ut1(&self, ut1: TimeDelta) -> Result<TimeDelta, EopProviderError> {
343        let seconds = ut1.seconds_since_j2000();
344        let (t0, _) = self.delta_ut1_tai.first();
345        let (tn, _) = self.delta_ut1_tai.last();
346        // Use the UT1 offset as an initial guess even though the table is based on TAI
347        let mut val = self.delta_ut1_tai.interpolate(seconds);
348        // Interpolate again with the adjusted offsets
349        for _ in 0..2 {
350            val = self.delta_ut1_tai.interpolate(seconds - val);
351        }
352        if seconds < t0 || seconds > tn {
353            return Err(EopProviderError::ExtrapolatedValue(val));
354        }
355        Ok(-TimeDelta::from_seconds_f64(val))
356    }
357
358    /// Returns the leap seconds kernel, if one was configured.
359    pub fn get_lsk(&self) -> Option<&LeapSecondsKernel> {
360        self.lsk.as_ref()
361    }
362}
363
364#[cfg(test)]
365mod tests {
366    use lox_test_utils::data_file;
367    use lox_time::{Time, time_scales::Tai};
368
369    use super::*;
370
371    #[test]
372    #[should_panic(expected = "NoFiles")]
373    fn test_eop_parser_no_files() {
374        EopParser::new().parse().unwrap();
375    }
376
377    #[test]
378    #[should_panic(expected = "MissingIau2000")]
379    fn test_eop_provider_missing_iau2000() {
380        let t: Time<Tai> = Time::default();
381        let eop = EopParser::new()
382            .from_path(data_file("iers/finals.all.csv"))
383            .parse()
384            .unwrap();
385        eop.nutation_precession_iau2000(t).unwrap();
386    }
387
388    #[test]
389    #[should_panic(expected = "MissingIau1980")]
390    fn test_eop_provider_missing_iau1980() {
391        let t: Time<Tai> = Time::default();
392        let eop = EopParser::new()
393            .from_path(data_file("iers/finals2000A.all.csv"))
394            .parse()
395            .unwrap();
396        eop.nutation_precession_iau1980(t).unwrap();
397    }
398}
lox_earth/eop.rs

lox_earth/
eop.rs
