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

sidereon_core/tides/
mod.rs
