sidereon_core/tides/
mod.rs

1//! Solid-earth tide station displacement (IERS Conventions, Chapter 7).
2//!
3//! [`solid_earth_tide`] computes the tidal displacement of an Earth-fixed (ITRF)
4//! GNSS station caused by the lunar and solar gravitational attraction. It is a
5//! derived work of the IERS Conventions (2010) reference routine
6//! `DEHANTTIDEINEL.F` (and its companion routines `ST1IDIU`, `ST1ISEM`,
7//! `ST1L1`, `STEP2DIU`, `STEP2LON`, `CAL2JD`, `DAT`), reproduced here in Rust.
8//!
9//! How this derived work is based upon and differs from the original Software:
10//!
11//! * It is a line-for-line Rust translation of the in-phase degree-2/degree-3
12//!   displacement, the out-of-phase corrections (`ST1IDIU`, `ST1ISEM`), the
13//!   latitude-dependence correction (`ST1L1`), and the frequency-dependent
14//!   step-2 diurnal/long-period band corrections (`STEP2DIU`, `STEP2LON`),
15//!   evaluating the identical Love/Shida numbers, Doodson/argument tables, and
16//!   leap-second table.
17//! * It keeps the permanent (mean) tide deformation: the original routine's
18//!   commented-out "Step 3" permanent-tide removal is left disabled, matching
19//!   the ITRF/IGS conform-to-mean-tide convention.
20//! * The routine names are changed from the IERS originals (per the IERS
21//!   Conventions Software License), and the Fortran subroutine structure is
22//!   inlined into private helpers.
23//!
24//! The Sun and Moon geocentric positions are inputs (metres, ECEF/ITRF); the
25//! caller supplies them, e.g. from [`crate::astro::bodies::sun_moon_ecef`].
26//!
27//! IERS Conventions Software License: permission is granted to use this software
28//! for any purpose, including commercial applications, free of charge, and to
29//! distribute derived works subject to the conditions reproduced above. Results
30//! obtained with this software acknowledge use of the IERS Conventions software.
31
32#[cfg(all(test, sidereon_repo_tests))]
33mod tests;
34
35mod ocean;
36mod pole;
37pub use ocean::{ocean_tide_loading, OceanLoadingBlq, NUM_OCEAN_CONSTITUENTS};
38pub use pole::solid_earth_pole_tide;
39
40use crate::astro::constants::models::iers::SOLID_TIDE_EARTH_RADIUS_M;
41use crate::astro::constants::time::{
42    DAYS_PER_JULIAN_CENTURY, J2000_JD, SECONDS_PER_DAY, TT_MINUS_TAI_S,
43};
44use crate::astro::constants::units::{DEG_TO_RAD, KM_TO_M};
45use crate::astro::math::vec3::{dot3_ref as dot, norm3_ref as norm8};
46use crate::validate::{self, FieldError};
47
48#[derive(Debug, Clone, Copy, PartialEq, Eq)]
49pub enum TideInputErrorKind {
50    Missing,
51    NonFinite,
52    NotPositive,
53    Negative,
54    OutOfRange,
55    FloatParse,
56    IntParse,
57    InvalidCivilDate,
58    InvalidCivilTime,
59}
60
61impl core::fmt::Display for TideInputErrorKind {
62    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
63        f.write_str(match self {
64            Self::Missing => "missing",
65            Self::NonFinite => "not finite",
66            Self::NotPositive => "not positive",
67            Self::Negative => "negative",
68            Self::OutOfRange => "out of range",
69            Self::FloatParse => "invalid float",
70            Self::IntParse => "invalid integer",
71            Self::InvalidCivilDate => "invalid civil date",
72            Self::InvalidCivilTime => "invalid civil time",
73        })
74    }
75}
76
77impl From<&FieldError> for TideInputErrorKind {
78    fn from(error: &FieldError) -> Self {
79        match error {
80            FieldError::Missing { .. } => Self::Missing,
81            FieldError::NonFinite { .. } => Self::NonFinite,
82            FieldError::NotPositive { .. } => Self::NotPositive,
83            FieldError::Negative { .. } => Self::Negative,
84            FieldError::OutOfRange { .. } => Self::OutOfRange,
85            FieldError::FloatParse { .. } => Self::FloatParse,
86            FieldError::IntParse { .. } => Self::IntParse,
87            FieldError::InvalidCivilDate { .. } => Self::InvalidCivilDate,
88            FieldError::InvalidCivilTime { .. } => Self::InvalidCivilTime,
89        }
90    }
91}
92
93#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
94pub enum TideError {
95    #[error("invalid solid-earth tide input {field}: {kind}")]
96    InvalidInput {
97        field: &'static str,
98        kind: TideInputErrorKind,
99    },
100}
101
102fn invalid_tide_input(error: FieldError) -> TideError {
103    TideError::InvalidInput {
104        field: error.field(),
105        kind: (&error).into(),
106    }
107}
108
109/// Solid-earth tide displacement of an ITRF station, in metres (ECEF).
110///
111/// Arguments mirror the IERS reference routine:
112/// * `xsta` - geocentric station position (m, ITRF).
113/// * `year`, `month`, `day` - UTC calendar date.
114/// * `fhr` - UTC fractional hour of the day (hour + min/60 + sec/3600).
115/// * `xsun` - geocentric Sun position (m, ECEF).
116/// * `xmon` - geocentric Moon position (m, ECEF).
117///
118/// Returns the displacement vector `dxtide` (m, geocentric ITRF). The permanent
119/// (mean) tide deformation is retained (ITRF/IGS convention).
120///
121/// Returns [`TideError`] when inputs are non-finite or geometrically
122/// degenerate: the station vector must be non-zero and non-polar, and Sun/Moon
123/// vectors must be non-zero.
124pub fn solid_earth_tide(
125    xsta: &[f64; 3],
126    year: i32,
127    month: i32,
128    day: i32,
129    fhr: f64,
130    xsun: &[f64; 3],
131    xmon: &[f64; 3],
132) -> Result<[f64; 3], TideError> {
133    validate_tide_domain(xsta, year, month, day, fhr, xsun, xmon)?;
134    Ok(solid_earth_tide_unchecked(
135        xsta, year, month, day, fhr, xsun, xmon,
136    ))
137}
138
139fn validate_tide_domain(
140    xsta: &[f64; 3],
141    year: i32,
142    month: i32,
143    day: i32,
144    fhr: f64,
145    xsun: &[f64; 3],
146    xmon: &[f64; 3],
147) -> Result<(), TideError> {
148    validate::finite_vec3(*xsta, "station position").map_err(invalid_tide_input)?;
149    validate::civil_datetime_with_second_policy(
150        i64::from(year),
151        i64::from(month),
152        i64::from(day),
153        0,
154        0,
155        0.0,
156        validate::CivilSecondPolicy::Continuous,
157    )
158    .map_err(invalid_tide_input)?;
159    validate::finite_in_range_exclusive_upper(fhr, 0.0, 24.0, "fractional hour")
160        .map_err(invalid_tide_input)?;
161    validate::finite_vec3(*xsun, "sun position").map_err(invalid_tide_input)?;
162    validate::finite_vec3(*xmon, "moon position").map_err(invalid_tide_input)?;
163
164    validate::finite_positive(norm8(xsta), "station radius").map_err(invalid_tide_input)?;
165    let station_horizontal_radius = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt();
166    validate::finite_positive(station_horizontal_radius, "station horizontal radius")
167        .map_err(invalid_tide_input)?;
168    validate::finite_positive(norm8(xsun), "sun radius").map_err(invalid_tide_input)?;
169    validate::finite_positive(norm8(xmon), "moon radius").map_err(invalid_tide_input)?;
170
171    Ok(())
172}
173
174fn solid_earth_tide_unchecked(
175    xsta: &[f64; 3],
176    year: i32,
177    month: i32,
178    day: i32,
179    fhr: f64,
180    xsun: &[f64; 3],
181    xmon: &[f64; 3],
182) -> [f64; 3] {
183    // Nominal second- and third-degree Love and Shida numbers.
184    const H20: f64 = 0.6078;
185    const L20: f64 = 0.0847;
186    const H3: f64 = 0.292;
187    const L3: f64 = 0.015;
188
189    // Scalar product of station vector with Sun/Moon vector (SPROD).
190    let rsta = norm8(xsta);
191    let rsun = norm8(xsun);
192    let rmon = norm8(xmon);
193    let scs = dot(xsta, xsun);
194    let scm = dot(xsta, xmon);
195    let scsun = scs / rsta / rsun;
196    let scmon = scm / rsta / rmon;
197
198    // Latitude-corrected H2 and L2.
199    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
200    let h2 = H20 - 0.0006 * (1.0 - 3.0 / 2.0 * cosphi * cosphi);
201    let l2 = L20 + 0.0002 * (1.0 - 3.0 / 2.0 * cosphi * cosphi);
202
203    // P2 term.
204    let p2sun = 3.0 * (h2 / 2.0 - l2) * scsun * scsun - h2 / 2.0;
205    let p2mon = 3.0 * (h2 / 2.0 - l2) * scmon * scmon - h2 / 2.0;
206
207    // P3 term.
208    let p3sun = 5.0 / 2.0 * (H3 - 3.0 * L3) * scsun.powi(3) + 3.0 / 2.0 * (L3 - H3) * scsun;
209    let p3mon = 5.0 / 2.0 * (H3 - 3.0 * L3) * scmon.powi(3) + 3.0 / 2.0 * (L3 - H3) * scmon;
210
211    // Term in direction of Sun/Moon vector.
212    let x2sun = 3.0 * l2 * scsun;
213    let x2mon = 3.0 * l2 * scmon;
214    let x3sun = 3.0 * L3 / 2.0 * (5.0 * scsun * scsun - 1.0);
215    let x3mon = 3.0 * L3 / 2.0 * (5.0 * scmon * scmon - 1.0);
216
217    // Factors for Sun/Moon (IAU current best estimates).
218    const MASS_RATIO_SUN: f64 = 332946.0482;
219    const MASS_RATIO_MOON: f64 = 0.0123000371;
220    const RE: f64 = SOLID_TIDE_EARTH_RADIUS_M;
221    let fac2sun = MASS_RATIO_SUN * RE * (RE / rsun).powi(3);
222    let fac2mon = MASS_RATIO_MOON * RE * (RE / rmon).powi(3);
223    let fac3sun = fac2sun * (RE / rsun);
224    let fac3mon = fac2mon * (RE / rmon);
225
226    // Total in-phase degree-2/degree-3 displacement.
227    let mut dxtide = [0.0_f64; 3];
228    for i in 0..3 {
229        dxtide[i] = fac2sun * (x2sun * xsun[i] / rsun + p2sun * xsta[i] / rsta)
230            + fac2mon * (x2mon * xmon[i] / rmon + p2mon * xsta[i] / rsta)
231            + fac3sun * (x3sun * xsun[i] / rsun + p3sun * xsta[i] / rsta)
232            + fac3mon * (x3mon * xmon[i] / rmon + p3mon * xsta[i] / rsta);
233    }
234
235    // Out-of-phase corrections (diurnal, semi-diurnal) and latitude dependence.
236    let c = st1idiu(xsta, xsun, xmon, fac2sun, fac2mon);
237    for i in 0..3 {
238        dxtide[i] += c[i];
239    }
240    let c = st1isem(xsta, xsun, xmon, fac2sun, fac2mon);
241    for i in 0..3 {
242        dxtide[i] += c[i];
243    }
244    let c = st1l1(xsta, xsun, xmon, fac2sun, fac2mon);
245    for i in 0..3 {
246        dxtide[i] += c[i];
247    }
248
249    // Step 2 corrections need the date in Julian centuries (TT).
250    let (jjm0, jjm1) = cal2jd(year, month, day);
251    let fhrd = fhr / 24.0;
252    let mut t = ((jjm0 - J2000_JD) + jjm1 + fhrd) / DAYS_PER_JULIAN_CENTURY;
253    let dtt = dat(year, month, day) + TT_MINUS_TAI_S;
254    t += dtt / (SECONDS_PER_DAY * DAYS_PER_JULIAN_CENTURY);
255
256    let c = step2diu(xsta, fhr, t);
257    for i in 0..3 {
258        dxtide[i] += c[i];
259    }
260    let c = step2lon(xsta, t);
261    for i in 0..3 {
262        dxtide[i] += c[i];
263    }
264
265    // Step 3 of the IERS routine, the permanent (zero-frequency) tide removal,
266    // is intentionally not applied, so the permanent (mean) tide deformation is
267    // retained (the ITRF/IGS conform-to-mean-tide convention; see module docs).
268    dxtide
269}
270
271/// Out-of-phase part of the Love numbers, diurnal band (ST1IDIU).
272fn st1idiu(
273    xsta: &[f64; 3],
274    xsun: &[f64; 3],
275    xmon: &[f64; 3],
276    fac2sun: f64,
277    fac2mon: f64,
278) -> [f64; 3] {
279    const DHI: f64 = -0.0025;
280    const DLI: f64 = -0.0007;
281    let rsta = norm8(xsta);
282    let sinphi = xsta[2] / rsta;
283    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
284    let cos2phi = cosphi * cosphi - sinphi * sinphi;
285    let sinla = xsta[1] / cosphi / rsta;
286    let cosla = xsta[0] / cosphi / rsta;
287    let rmon = norm8(xmon);
288    let rsun = norm8(xsun);
289
290    let drsun =
291        -3.0 * DHI * sinphi * cosphi * fac2sun * xsun[2] * (xsun[0] * sinla - xsun[1] * cosla)
292            / (rsun * rsun);
293    let drmon =
294        -3.0 * DHI * sinphi * cosphi * fac2mon * xmon[2] * (xmon[0] * sinla - xmon[1] * cosla)
295            / (rmon * rmon);
296    let dnsun = -3.0 * DLI * cos2phi * fac2sun * xsun[2] * (xsun[0] * sinla - xsun[1] * cosla)
297        / (rsun * rsun);
298    let dnmon = -3.0 * DLI * cos2phi * fac2mon * xmon[2] * (xmon[0] * sinla - xmon[1] * cosla)
299        / (rmon * rmon);
300    let desun = -3.0 * DLI * sinphi * fac2sun * xsun[2] * (xsun[0] * cosla + xsun[1] * sinla)
301        / (rsun * rsun);
302    let demon = -3.0 * DLI * sinphi * fac2mon * xmon[2] * (xmon[0] * cosla + xmon[1] * sinla)
303        / (rmon * rmon);
304
305    let dr = drsun + drmon;
306    let dn = dnsun + dnmon;
307    let de = desun + demon;
308
309    [
310        dr * cosla * cosphi - de * sinla - dn * sinphi * cosla,
311        dr * sinla * cosphi + de * cosla - dn * sinphi * sinla,
312        dr * sinphi + dn * cosphi,
313    ]
314}
315
316/// Out-of-phase part of the Love numbers, semi-diurnal band (ST1ISEM).
317fn st1isem(
318    xsta: &[f64; 3],
319    xsun: &[f64; 3],
320    xmon: &[f64; 3],
321    fac2sun: f64,
322    fac2mon: f64,
323) -> [f64; 3] {
324    const DHI: f64 = -0.0022;
325    const DLI: f64 = -0.0007;
326    let rsta = norm8(xsta);
327    let sinphi = xsta[2] / rsta;
328    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
329    let sinla = xsta[1] / cosphi / rsta;
330    let cosla = xsta[0] / cosphi / rsta;
331    let costwola = cosla * cosla - sinla * sinla;
332    let sintwola = 2.0 * cosla * sinla;
333    let rmon = norm8(xmon);
334    let rsun = norm8(xsun);
335
336    let drsun = -3.0 / 4.0
337        * DHI
338        * cosphi
339        * cosphi
340        * fac2sun
341        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * sintwola - 2.0 * xsun[0] * xsun[1] * costwola)
342        / (rsun * rsun);
343    let drmon = -3.0 / 4.0
344        * DHI
345        * cosphi
346        * cosphi
347        * fac2mon
348        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * sintwola - 2.0 * xmon[0] * xmon[1] * costwola)
349        / (rmon * rmon);
350    let dnsun = 3.0 / 2.0
351        * DLI
352        * sinphi
353        * cosphi
354        * fac2sun
355        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * sintwola - 2.0 * xsun[0] * xsun[1] * costwola)
356        / (rsun * rsun);
357    let dnmon = 3.0 / 2.0
358        * DLI
359        * sinphi
360        * cosphi
361        * fac2mon
362        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * sintwola - 2.0 * xmon[0] * xmon[1] * costwola)
363        / (rmon * rmon);
364    let desun = -3.0 / 2.0
365        * DLI
366        * cosphi
367        * fac2sun
368        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * costwola + 2.0 * xsun[0] * xsun[1] * sintwola)
369        / (rsun * rsun);
370    let demon = -3.0 / 2.0
371        * DLI
372        * cosphi
373        * fac2mon
374        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * costwola + 2.0 * xmon[0] * xmon[1] * sintwola)
375        / (rmon * rmon);
376
377    let dr = drsun + drmon;
378    let dn = dnsun + dnmon;
379    let de = desun + demon;
380
381    [
382        dr * cosla * cosphi - de * sinla - dn * sinphi * cosla,
383        dr * sinla * cosphi + de * cosla - dn * sinphi * sinla,
384        dr * sinphi + dn * cosphi,
385    ]
386}
387
388/// Latitude dependence of the Love numbers, part L^(1) (ST1L1).
389fn st1l1(
390    xsta: &[f64; 3],
391    xsun: &[f64; 3],
392    xmon: &[f64; 3],
393    fac2sun: f64,
394    fac2mon: f64,
395) -> [f64; 3] {
396    const L1D: f64 = 0.0012;
397    const L1SD: f64 = 0.0024;
398    let rsta = norm8(xsta);
399    let sinphi = xsta[2] / rsta;
400    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
401    let sinla = xsta[1] / cosphi / rsta;
402    let cosla = xsta[0] / cosphi / rsta;
403    let rmon = norm8(xmon);
404    let rsun = norm8(xsun);
405
406    // Diurnal band.
407    let mut l1 = L1D;
408    let dnsun = -l1 * sinphi * sinphi * fac2sun * xsun[2] * (xsun[0] * cosla + xsun[1] * sinla)
409        / (rsun * rsun);
410    let dnmon = -l1 * sinphi * sinphi * fac2mon * xmon[2] * (xmon[0] * cosla + xmon[1] * sinla)
411        / (rmon * rmon);
412    let desun = l1
413        * sinphi
414        * (cosphi * cosphi - sinphi * sinphi)
415        * fac2sun
416        * xsun[2]
417        * (xsun[0] * sinla - xsun[1] * cosla)
418        / (rsun * rsun);
419    let demon = l1
420        * sinphi
421        * (cosphi * cosphi - sinphi * sinphi)
422        * fac2mon
423        * xmon[2]
424        * (xmon[0] * sinla - xmon[1] * cosla)
425        / (rmon * rmon);
426
427    let de = 3.0 * (desun + demon);
428    let dn = 3.0 * (dnsun + dnmon);
429
430    let mut xcorsta = [
431        -de * sinla - dn * sinphi * cosla,
432        de * cosla - dn * sinphi * sinla,
433        dn * cosphi,
434    ];
435
436    // Semi-diurnal band.
437    l1 = L1SD;
438    let costwola = cosla * cosla - sinla * sinla;
439    let sintwola = 2.0 * cosla * sinla;
440
441    let dnsun = -l1 / 2.0
442        * sinphi
443        * cosphi
444        * fac2sun
445        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * costwola + 2.0 * xsun[0] * xsun[1] * sintwola)
446        / (rsun * rsun);
447    let dnmon = -l1 / 2.0
448        * sinphi
449        * cosphi
450        * fac2mon
451        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * costwola + 2.0 * xmon[0] * xmon[1] * sintwola)
452        / (rmon * rmon);
453    let desun = -l1 / 2.0
454        * sinphi
455        * sinphi
456        * cosphi
457        * fac2sun
458        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * sintwola - 2.0 * xsun[0] * xsun[1] * costwola)
459        / (rsun * rsun);
460    let demon = -l1 / 2.0
461        * sinphi
462        * sinphi
463        * cosphi
464        * fac2mon
465        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * sintwola - 2.0 * xmon[0] * xmon[1] * costwola)
466        / (rmon * rmon);
467
468    let de = 3.0 * (desun + demon);
469    let dn = 3.0 * (dnsun + dnmon);
470
471    xcorsta[0] += -de * sinla - dn * sinphi * cosla;
472    xcorsta[1] += de * cosla - dn * sinphi * sinla;
473    xcorsta[2] += dn * cosphi;
474    xcorsta
475}
476
477/// In-phase / out-of-phase frequency-dependent corrections, diurnal band
478/// (STEP2DIU). `fhr` is UTC fractional hour, `t` is Julian centuries (TT).
479fn step2diu(xsta: &[f64; 3], fhr: f64, t: f64) -> [f64; 3] {
480    // DATDI(9,31): {l, l', F, D, Omega(Ps), Adr, Adi, Anr, Ani} per wave.
481    #[rustfmt::skip]
482    const DATDI: [[f64; 9]; 31] = [
483        [-3.0, 0.0, 2.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
484        [-3.0, 2.0, 0.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
485        [-2.0, 0.0, 1.0, -1.0, 0.0, -0.02, 0.0, 0.0, 0.0],
486        [-2.0, 0.0, 1.0, 0.0, 0.0, -0.08, 0.0, -0.01, 0.01],
487        [-2.0, 2.0, -1.0, 0.0, 0.0, -0.02, 0.0, 0.0, 0.0],
488        [-1.0, 0.0, 0.0, -1.0, 0.0, -0.10, 0.0, 0.0, 0.0],
489        [-1.0, 0.0, 0.0, 0.0, 0.0, -0.51, 0.0, -0.02, 0.03],
490        [-1.0, 2.0, 0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0],
491        [0.0, -2.0, 1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0],
492        [0.0, 0.0, -1.0, 0.0, 0.0, 0.02, 0.0, 0.0, 0.0],
493        [0.0, 0.0, 1.0, 0.0, 0.0, 0.06, 0.0, 0.0, 0.0],
494        [0.0, 0.0, 1.0, 1.0, 0.0, 0.01, 0.0, 0.0, 0.0],
495        [0.0, 2.0, -1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0],
496        [1.0, -3.0, 0.0, 0.0, 1.0, -0.06, 0.0, 0.0, 0.0],
497        [1.0, -2.0, 0.0, -1.0, 0.0, 0.01, 0.0, 0.0, 0.0],
498        [1.0, -2.0, 0.0, 0.0, 0.0, -1.23, -0.07, 0.06, 0.01],
499        [1.0, -1.0, 0.0, 0.0, -1.0, 0.02, 0.0, 0.0, 0.0],
500        [1.0, -1.0, 0.0, 0.0, 1.0, 0.04, 0.0, 0.0, 0.0],
501        [1.0, 0.0, 0.0, -1.0, 0.0, -0.22, 0.01, 0.01, 0.0],
502        [1.0, 0.0, 0.0, 0.0, 0.0, 12.00, -0.80, -0.67, -0.03],
503        [1.0, 0.0, 0.0, 1.0, 0.0, 1.73, -0.12, -0.10, 0.0],
504        [1.0, 0.0, 0.0, 2.0, 0.0, -0.04, 0.0, 0.0, 0.0],
505        [1.0, 1.0, 0.0, 0.0, -1.0, -0.50, -0.01, 0.03, 0.0],
506        [1.0, 1.0, 0.0, 0.0, 1.0, 0.01, 0.0, 0.0, 0.0],
507        [0.0, 1.0, 0.0, 1.0, -1.0, -0.01, 0.0, 0.0, 0.0],
508        [1.0, 2.0, -2.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
509        [1.0, 2.0, 0.0, 0.0, 0.0, -0.11, 0.01, 0.01, 0.0],
510        [2.0, -2.0, 1.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
511        [2.0, 0.0, -1.0, 0.0, 0.0, -0.02, 0.0, 0.0, 0.0],
512        [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
513        [3.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
514    ];
515    let mut s = 218.31664563 + (481267.88194 + (-0.0014663889 + 0.00000185139 * t) * t) * t;
516    let mut tau = fhr * 15.0
517        + 280.4606184
518        + (36000.7700536 + (0.00038793 + -0.0000000258 * t) * t) * t
519        + (-s);
520    let pr = (1.396971278 + (0.000308889 + (0.000000021 + 0.000000007 * t) * t) * t) * t;
521    s += pr;
522    let mut h = 280.46645
523        + (36000.7697489 + (0.00030322222 + (0.000000020 + -0.00000000654 * t) * t) * t) * t;
524    let mut p = 83.35324312
525        + (4069.01363525 + (-0.01032172222 + (-0.0000124991 + 0.00000005263 * t) * t) * t) * t;
526    let mut zns = 234.95544499
527        + (1934.13626197 + (-0.00207561111 + (-0.00000213944 + 0.00000001650 * t) * t) * t) * t;
528    let mut ps = 282.93734098
529        + (1.71945766667 + (0.00045688889 + (-0.00000001778 + -0.00000000334 * t) * t) * t) * t;
530
531    s = s.rem_euclid(360.0);
532    tau = tau.rem_euclid(360.0);
533    h = h.rem_euclid(360.0);
534    p = p.rem_euclid(360.0);
535    zns = zns.rem_euclid(360.0);
536    ps = ps.rem_euclid(360.0);
537
538    let rsta = (xsta[0] * xsta[0] + xsta[1] * xsta[1] + xsta[2] * xsta[2]).sqrt();
539    let sinphi = xsta[2] / rsta;
540    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
541    let cosla = xsta[0] / cosphi / rsta;
542    let sinla = xsta[1] / cosphi / rsta;
543    let zla = xsta[1].atan2(xsta[0]);
544
545    let mut xcorsta = [0.0_f64; 3];
546    for w in &DATDI {
547        let thetaf = (tau + w[0] * s + w[1] * h + w[2] * p + w[3] * zns + w[4] * ps) * DEG_TO_RAD;
548        let dr = w[5] * 2.0 * sinphi * cosphi * (thetaf + zla).sin()
549            + w[6] * 2.0 * sinphi * cosphi * (thetaf + zla).cos();
550        let dn = w[7] * (cosphi * cosphi - sinphi * sinphi) * (thetaf + zla).sin()
551            + w[8] * (cosphi * cosphi - sinphi * sinphi) * (thetaf + zla).cos();
552        let de = w[7] * sinphi * (thetaf + zla).cos() - w[8] * sinphi * (thetaf + zla).sin();
553
554        xcorsta[0] += dr * cosla * cosphi - de * sinla - dn * sinphi * cosla;
555        xcorsta[1] += dr * sinla * cosphi + de * cosla - dn * sinphi * sinla;
556        xcorsta[2] += dr * sinphi + dn * cosphi;
557    }
558    for v in &mut xcorsta {
559        *v /= KM_TO_M;
560    }
561    xcorsta
562}
563
564/// In-phase / out-of-phase frequency-dependent corrections, long-period band
565/// (STEP2LON). `t` is Julian centuries (TT).
566fn step2lon(xsta: &[f64; 3], t: f64) -> [f64; 3] {
567    #[rustfmt::skip]
568    const DATDI: [[f64; 9]; 5] = [
569        [0.0, 0.0, 0.0, 1.0, 0.0, 0.47, 0.23, 0.16, 0.07],
570        [0.0, 2.0, 0.0, 0.0, 0.0, -0.20, -0.12, -0.11, -0.05],
571        [1.0, 0.0, -1.0, 0.0, 0.0, -0.11, -0.08, -0.09, -0.04],
572        [2.0, 0.0, 0.0, 0.0, 0.0, -0.13, -0.11, -0.15, -0.07],
573        [2.0, 0.0, 0.0, 1.0, 0.0, -0.05, -0.05, -0.06, -0.03],
574    ];
575    let mut s = 218.31664563 + (481267.88194 + (-0.0014663889 + 0.00000185139 * t) * t) * t;
576    let pr = (1.396971278 + (0.000308889 + (0.000000021 + 0.000000007 * t) * t) * t) * t;
577    s += pr;
578    let mut h = 280.46645
579        + (36000.7697489 + (0.00030322222 + (0.000000020 + -0.00000000654 * t) * t) * t) * t;
580    let mut p = 83.35324312
581        + (4069.01363525 + (-0.01032172222 + (-0.0000124991 + 0.00000005263 * t) * t) * t) * t;
582    let mut zns = 234.95544499
583        + (1934.13626197 + (-0.00207561111 + (-0.00000213944 + 0.00000001650 * t) * t) * t) * t;
584    let mut ps = 282.93734098
585        + (1.71945766667 + (0.00045688889 + (-0.00000001778 + -0.00000000334 * t) * t) * t) * t;
586
587    let rsta = (xsta[0] * xsta[0] + xsta[1] * xsta[1] + xsta[2] * xsta[2]).sqrt();
588    let sinphi = xsta[2] / rsta;
589    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
590    let cosla = xsta[0] / cosphi / rsta;
591    let sinla = xsta[1] / cosphi / rsta;
592
593    s = s.rem_euclid(360.0);
594    h = h.rem_euclid(360.0);
595    p = p.rem_euclid(360.0);
596    zns = zns.rem_euclid(360.0);
597    ps = ps.rem_euclid(360.0);
598
599    let mut xcorsta = [0.0_f64; 3];
600    for w in &DATDI {
601        let thetaf = (w[0] * s + w[1] * h + w[2] * p + w[3] * zns + w[4] * ps) * DEG_TO_RAD;
602        let dr = w[5] * (3.0 * sinphi * sinphi - 1.0) / 2.0 * thetaf.cos()
603            + w[7] * (3.0 * sinphi * sinphi - 1.0) / 2.0 * thetaf.sin();
604        let dn = w[6] * (cosphi * sinphi * 2.0) * thetaf.cos()
605            + w[8] * (cosphi * sinphi * 2.0) * thetaf.sin();
606        let de = 0.0;
607
608        xcorsta[0] += dr * cosla * cosphi - de * sinla - dn * sinphi * cosla;
609        xcorsta[1] += dr * sinla * cosphi + de * cosla - dn * sinphi * sinla;
610        xcorsta[2] += dr * sinphi + dn * cosphi;
611    }
612    for v in &mut xcorsta {
613        *v /= KM_TO_M;
614    }
615    xcorsta
616}
617
618/// Gregorian calendar date -> (MJD epoch 2400000.5, MJD) (SOFA CAL2JD).
619///
620/// This is a SOFA parity adapter, deliberately NOT routed through
621/// [`crate::astro::time::civil`]: the solid-Earth/ocean/pole tide models are
622/// validated bit-for-bit against the IERS/SOFA reference (the
623/// `ocean_loading_oracle` test), so the calendar-to-MJD step must reproduce
624/// SOFA's `iauCal2jd` exactly. It is kept local under this tides-specific name
625/// so it is not mistaken for a duplicate of the canonical civil conversions and
626/// is never consolidated into them.
627fn cal2jd(iy: i32, im: i32, id: i32) -> (f64, f64) {
628    let my = (im - 14) / 12;
629    let iypmy = iy + my;
630    let djm0 = 2400000.5;
631    let djm = ((1461 * (iypmy + 4800)) / 4 + (367 * (im - 2 - 12 * my)) / 12
632        - (3 * ((iypmy + 4900) / 100)) / 4
633        + id
634        - 2432076) as f64;
635    (djm0, djm)
636}
637
638/// TAI-UTC (Delta(AT)) in seconds for the given date (SOFA DAT, leap-second
639/// table only; the four golden dates are all post-1972 so the pre-1972 drift
640/// branch is not exercised, but it is retained for completeness).
641fn dat(iy: i32, im: i32, _id: i32) -> f64 {
642    // Post-1972 leap-second table: (year, month, Delta(AT) seconds).
643    const IDAT: [(i32, i32, f64); 28] = [
644        (1972, 1, 10.0),
645        (1972, 7, 11.0),
646        (1973, 1, 12.0),
647        (1974, 1, 13.0),
648        (1975, 1, 14.0),
649        (1976, 1, 15.0),
650        (1977, 1, 16.0),
651        (1978, 1, 17.0),
652        (1979, 1, 18.0),
653        (1980, 1, 19.0),
654        (1981, 7, 20.0),
655        (1982, 7, 21.0),
656        (1983, 7, 22.0),
657        (1985, 7, 23.0),
658        (1988, 1, 24.0),
659        (1990, 1, 25.0),
660        (1991, 1, 26.0),
661        (1992, 7, 27.0),
662        (1993, 7, 28.0),
663        (1994, 7, 29.0),
664        (1996, 1, 30.0),
665        (1997, 7, 31.0),
666        (1999, 1, 32.0),
667        (2006, 1, 33.0),
668        (2009, 1, 34.0),
669        (2012, 7, 35.0),
670        (2015, 7, 36.0),
671        (2017, 1, 37.0),
672    ];
673    let m = 12 * iy + im;
674    let mut da = IDAT[0].2;
675    for &(y, mo, d) in &IDAT {
676        if m >= 12 * y + mo {
677            da = d;
678        }
679    }
680    da
681}
sidereon_core/tides/mod.rs

sidereon_core/tides/
mod.rs
