diff --git a/.idea/copyright/Lox.xml b/.idea/copyright/Lox.xml
index 3ca717b3..23676d50 100644
--- a/.idea/copyright/Lox.xml
+++ b/.idea/copyright/Lox.xml
@@ -1,6 +1,6 @@
-
+
\ No newline at end of file
diff --git a/.idea/lox-space.iml b/.idea/lox-space.iml
index 7e2df88b..85dd8d82 100644
--- a/.idea/lox-space.iml
+++ b/.idea/lox-space.iml
@@ -12,6 +12,8 @@
+
+
diff --git a/Cargo.lock b/Cargo.lock
index 2fa525d1..41c9add2 100644
--- a/Cargo.lock
+++ b/Cargo.lock
@@ -50,6 +50,27 @@ version = "1.0.0"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "baf1de4339761588bc0619e3cbc0120ee582ebb74b53b4efbf79117bd2da40fd"
+[[package]]
+name = "csv"
+version = "1.3.0"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "ac574ff4d437a7b5ad237ef331c17ccca63c46479e5b5453eb8e10bb99a759fe"
+dependencies = [
+ "csv-core",
+ "itoa",
+ "ryu",
+ "serde",
+]
+
+[[package]]
+name = "csv-core"
+version = "0.1.11"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "5efa2b3d7902f4b634a20cae3c9c4e6209dc4779feb6863329607560143efa70"
+dependencies = [
+ "memchr",
+]
+
[[package]]
name = "errno"
version = "0.3.5"
@@ -208,6 +229,12 @@ version = "2.0.4"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "1e186cfbae8084e513daff4240b4797e342f988cecda4fb6c939150f96315fd8"
+[[package]]
+name = "itoa"
+version = "1.0.9"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "af150ab688ff2122fcef229be89cb50dd66af9e01a4ff320cc137eecc9bacc38"
+
[[package]]
name = "lazy_static"
version = "1.4.0"
@@ -242,6 +269,17 @@ dependencies = [
"scopeguard",
]
+[[package]]
+name = "lox-eop"
+version = "0.1.0"
+dependencies = [
+ "csv",
+ "float_eq",
+ "num-traits",
+ "serde",
+ "thiserror",
+]
+
[[package]]
name = "lox-space"
version = "0.0.0"
@@ -689,6 +727,12 @@ dependencies = [
"wait-timeout",
]
+[[package]]
+name = "ryu"
+version = "1.0.15"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "1ad4cc8da4ef723ed60bced201181d83791ad433213d8c24efffda1eec85d741"
+
[[package]]
name = "scopeguard"
version = "1.2.0"
@@ -701,6 +745,26 @@ version = "1.0.20"
source = "registry+https://github.com/rust-lang/crates.io-index"
checksum = "836fa6a3e1e547f9a2c4040802ec865b5d85f4014efe00555d7090a3dcaa1090"
+[[package]]
+name = "serde"
+version = "1.0.190"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "91d3c334ca1ee894a2c6f6ad698fe8c435b76d504b13d436f0685d648d6d96f7"
+dependencies = [
+ "serde_derive",
+]
+
+[[package]]
+name = "serde_derive"
+version = "1.0.190"
+source = "registry+https://github.com/rust-lang/crates.io-index"
+checksum = "67c5609f394e5c2bd7fc51efda478004ea80ef42fee983d5c67a65e34f32c0e3"
+dependencies = [
+ "proc-macro2",
+ "quote",
+ "syn",
+]
+
[[package]]
name = "slab"
version = "0.4.9"
diff --git a/crates/lox-eop/Cargo.toml b/crates/lox-eop/Cargo.toml
new file mode 100644
index 00000000..0687137d
--- /dev/null
+++ b/crates/lox-eop/Cargo.toml
@@ -0,0 +1,18 @@
+[package]
+name = "lox-eop"
+version = "0.1.0"
+rust-version.workspace = true
+edition.workspace = true
+license.workspace = true
+authors.workspace = true
+
+# See more keys and their definitions at https://doc.rust-lang.org/cargo/reference/manifest.html
+
+[dependencies]
+csv = "1.3.0"
+num-traits = "0.2.17"
+serde = { version = "1.0.190", features = ["derive"] }
+thiserror.workspace = true
+
+[dev-dependencies]
+float_eq = "1.0.1"
diff --git a/crates/lox-eop/src/akima.rs b/crates/lox-eop/src/akima.rs
new file mode 100644
index 00000000..e59251da
--- /dev/null
+++ b/crates/lox-eop/src/akima.rs
@@ -0,0 +1,231 @@
+/*
+ * Copyright (c) 2023. Helge Eichhorn and the LOX contributors
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+use std::f64;
+use std::iter::zip;
+use std::ops::Sub;
+
+use thiserror::Error;
+
+fn diff + Copy>(x: &Vec) -> Vec {
+ let n = x.len();
+ let mut dx = vec![];
+ for (i1, i2) in zip(0..=n - 2, 1..=n - 1) {
+ let dxi = x[i2] - x[i1];
+ dx.push(dxi);
+ }
+ dx
+}
+
+#[derive(Error, Debug)]
+pub enum AkimaError {
+ #[error("size of `x` and `y` must match")]
+ SizeMismatch,
+}
+
+pub struct Akima {
+ n: usize,
+ pub x: Vec,
+ pub y: Vec,
+ b: Vec,
+ c: Vec,
+ d: Vec,
+}
+
+impl Akima {
+ pub fn new(x: Vec, y: Vec) -> Result {
+ let n = x.len();
+ if n != y.len() {
+ return Err(AkimaError::SizeMismatch);
+ }
+
+ let mut b = vec![];
+ let mut c = vec![];
+ let mut d = vec![];
+
+ let dx = diff(&x);
+ let dy = diff(&y);
+ let mut m = vec![];
+ for (dxi, dyi) in zip(dx.clone(), dy) {
+ m.push(dyi / f64::from(dxi));
+ }
+ m.insert(0, 2.0 * m[0] - m[1]);
+ m.insert(0, 2.0 * m[0] - m[1]);
+ m.push(2.0 * m[m.len() - 1] - m[m.len() - 2]);
+ m.push(2.0 * m[m.len() - 1] - m[m.len() - 2]);
+
+ for (i1, i2) in zip(3..m.len(), 0..m.len() - 3) {
+ b.push(0.5 * (m[i1] + m[i2]));
+ }
+
+ let dm: Vec = diff(&m).iter().map(|x| x.abs()).collect();
+ let f1 = &dm[2..n + 2];
+ let f2 = &dm[0..n];
+ let f12: Vec = zip(f1, f2).map(|(x1, x2)| x1 + x2).collect();
+ let f12_max = f12.iter().cloned().fold(-1. / 0. /* inf */, f64::max);
+ let ind: Vec = f12
+ .iter()
+ .enumerate()
+ .filter(|(_, x)| **x > 1e-9 * f12_max)
+ .map(|(idx, _)| idx)
+ .collect();
+ for i in ind {
+ b[i] = (f1[i] * m[i + 1] + f2[i] * m[i + 2]) / f12[i]
+ }
+
+ for i in 0..n - 1 {
+ c.push((3.0 * m[i + 2] - 2.0 * b[i] - b[i + 1]) / f64::from(dx[i]))
+ }
+
+ for i in 0..n - 1 {
+ d.push((b[i] + b[i + 1] - 2.0 * m[i + 2]) / f64::from(dx[i].pow(2)))
+ }
+ Ok(Self {
+ n: x.len(),
+ x,
+ y,
+ b,
+ c,
+ d,
+ })
+ }
+
+ pub fn interpolate(&self, xi: f64) -> f64 {
+ let x = &self.x;
+ let x0 = f64::from(x[0]);
+ let xn = f64::from(x[x.len() - 1]);
+ if xi <= x0 {
+ return self.y[0];
+ }
+ if xi >= xn {
+ return self.y[self.n - 1];
+ }
+
+ let idx = x.binary_search(&(xi.trunc() as i32)).unwrap();
+
+ let wj = xi - f64::from(self.x[idx]);
+ let y = self.y[idx];
+ let b = self.b[idx];
+ let c = self.c[idx];
+ let d = self.d[idx];
+ d.mul_add(wj, c).mul_add(wj, b).mul_add(wj, y)
+ }
+}
+
+#[cfg(test)]
+mod tests {
+ use float_eq::assert_float_eq;
+
+ use super::*;
+
+ #[test]
+ fn test_diff() {
+ let x = vec![
+ -0.00561583609169947,
+ 0.29513581230551944,
+ -0.944145132062222,
+ 0.8539804645085572,
+ -0.6630410427468136,
+ -0.33045519762661285,
+ -0.5237166946868412,
+ -1.1435794359757951,
+ -0.5221715292393267,
+ 0.4762176135879527,
+ ];
+
+ let dx_exp = vec![
+ 0.3007516483972189,
+ -1.2392809443677413,
+ 1.7981255965707792,
+ -1.5170215072553708,
+ 0.3325858451202008,
+ -0.19326149706022838,
+ -0.6198627412889539,
+ 0.6214079067364684,
+ 0.9983891428272795,
+ ];
+
+ let dx_act = diff(&x);
+ assert_eq!(dx_act, dx_exp);
+ }
+
+ #[test]
+ fn test() {
+ let x = vec![0, 1, 2, 3, 4, 5, 6, 7, 8, 9, 10];
+ let y = vec![0.0, 2.0, 1.0, 3.0, 2.0, 6.0, 5.5, 5.5, 2.7, 5.1, 3.0];
+
+ let b = vec![
+ 3.5,
+ 0.5,
+ 0.5,
+ 0.875,
+ 1.0,
+ -0.09090909090909091,
+ -0.1917808219178082,
+ -0.2456140350877193,
+ -0.8054794520547944,
+ -0.01237113402061866,
+ -4.349999999999999,
+ ];
+
+ let c = vec![
+ -1.5,
+ -4.5,
+ 4.125,
+ -5.75,
+ 10.090909090909092,
+ -1.1264009962640098,
+ 0.6291756789233357,
+ -7.103292477769766,
+ 8.823330038130205,
+ -1.9252577319587632,
+ ];
+
+ let d = vec![
+ 0.0,
+ 3.0,
+ -2.625,
+ 3.875,
+ -7.090909090909091,
+ 0.7173100871731009,
+ -0.4373948570055275,
+ 4.548906512857486,
+ -5.617850586075412,
+ -0.16237113402061798,
+ ];
+
+ let xi = vec![
+ 0.0, 0.5, 1., 1.5, 2.5, 3.5, 4.5, 5.1, 6.5, 7.2, 8.6, 9.9, 10.,
+ ];
+ let yi = vec![
+ 0.,
+ 1.375,
+ 2.,
+ 1.5,
+ 1.953125,
+ 2.484375,
+ 4.136_363_636_363_637,
+ 5.980_362_391_033_624,
+ 5.506_729_151_646_239,
+ 5.203_136_745_974_525,
+ 4.179_655_415_901_708,
+ 3.411_038_659_793_813,
+ 3.0,
+ ];
+
+ let akima = Akima::new(x, y).expect("vectors lengths should match");
+ assert_eq!(akima.b, b);
+ assert_eq!(akima.c, c);
+ assert_eq!(akima.d, d);
+
+ for (xi, yi_exp) in zip(xi, yi) {
+ let yi_act = akima.interpolate(xi);
+ assert_float_eq!(yi_act, yi_exp, rel <= 1e-8);
+ }
+ }
+}
diff --git a/crates/lox-eop/src/lib.rs b/crates/lox-eop/src/lib.rs
new file mode 100644
index 00000000..93fb740a
--- /dev/null
+++ b/crates/lox-eop/src/lib.rs
@@ -0,0 +1,658 @@
+/*
+ * Copyright (c) 2023. Helge Eichhorn and the LOX contributors
+ *
+ * This Source Code Form is subject to the terms of the Mozilla Public
+ * License, v. 2.0. If a copy of the MPL was not distributed with this
+ * file, You can obtain one at http://mozilla.org/MPL/2.0/.
+ */
+
+use std::collections::HashMap;
+use std::fs;
+use std::path::Path;
+
+use serde::Deserialize;
+use thiserror::Error;
+
+use crate::akima::{Akima, AkimaError};
+
+mod akima;
+
+#[derive(Error, Debug)]
+enum LoxEopError {
+ #[error(transparent)]
+ Csv(#[from] csv::Error),
+ #[error(transparent)]
+ Io(#[from] std::io::Error),
+ #[error(transparent)]
+ Interpolation(#[from] AkimaError),
+}
+
+#[derive(Debug, Deserialize)]
+enum ValueType {
+ #[serde(rename = "prediction")]
+ Prediction,
+ #[serde(rename = "final")]
+ Final,
+}
+
+#[derive(Debug, Deserialize)]
+struct Record {
+ #[serde(rename = "MJD")]
+ modified_julian_date: i32,
+ // #[serde(rename = "Year")]
+ // year: u64,
+ // #[serde(rename = "Month")]
+ // month: u64,
+ // #[serde(rename = "Day")]
+ // day: u64,
+ // type_polar_motion: ValueType,
+ x_pole: Option,
+ sigma_x_pole: Option,
+ y_pole: Option,
+ sigma_y_pole: Option,
+ // type_rotation: String,
+ #[serde(rename = "UT1-UTC")]
+ delta_ut1_utc: Option,
+ #[serde(rename = "sigma_UT1-UTC")]
+ sigma_delta_ut1_utc: Option,
+ #[serde(rename = "LOD")]
+ lod: Option,
+ #[serde(rename = "sigma_LOD")]
+ sigma_lod: Option,
+ // type_equator: String,
+ #[serde(rename = "dPsi")]
+ delta_psi: Option,
+ #[serde(rename = "sigma_dPsi")]
+ sigma_delta_psi: Option,
+ #[serde(rename = "dEpsilon")]
+ delta_epsilon: Option,
+ #[serde(rename = "sigma_dEpsilon")]
+ sigma_delta_epsilon: Option,
+ #[serde(rename = "dX")]
+ delta_x: Option,
+ #[serde(rename = "sigma_dX")]
+ sigma_delta_x: Option,
+ #[serde(rename = "dY")]
+ delta_y: Option,
+ #[serde(rename = "sigma_dY")]
+ sigma_delta_y: Option,
+ // type_polar_motion_b: String,
+ // #[serde(rename = "bulB/x_pole")]
+ // x_pole_b: Option,
+ // #[serde(rename = "bulB/y_pole")]
+ // y_pole_b: Option,
+ // type_rotation_b: String,
+ // #[serde(rename = "bulB/UT-UTC")]
+ // delta_ut1_utc_b: Option,
+ // type_equator_b: String,
+ // #[serde(rename = "bulB/dPsi")]
+ // delta_psi_b: Option,
+ // #[serde(rename = "bulB/dEpsilon")]
+ // delta_epsilon_b: Option,
+ // #[serde(rename = "bulB/dX")]
+ // delta_x_b: Option,
+ // #[serde(rename = "bulB/dY")]
+ // delta_y_b: Option,
+}
+
+fn read_records>(path: P) -> Result, LoxEopError> {
+ let contents = fs::read_to_string(path)?
+ // Replace duplicate `Type` headers
+ .replacen("Type", "type_polar_motion", 1)
+ .replacen("Type", "type_rotation", 1)
+ .replacen("Type", "type_equator", 1)
+ .replacen("Type", "type_polar_motion_b", 1)
+ .replacen("Type", "type_rotation_b", 1)
+ .replacen("Type", "type_equator_b", 1);
+ let mut rdr = csv::ReaderBuilder::new()
+ .delimiter(b';')
+ .from_reader(contents.as_bytes());
+ let mut records: Vec = vec![];
+ for result in rdr.deserialize() {
+ records.push(result?);
+ }
+ Ok(records)
+}
+
+struct Records {
+ modified_julian_date: Vec,
+ x_pole: Vec