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::{
38    ocean_tide_loading, parse_ocean_loading_blq_block, parse_ocean_loading_blq_blocks,
39    OceanLoadingBlq, OceanLoadingBlqBlock, OceanTideConstituent, NUM_OCEAN_CONSTITUENTS,
40    OCEAN_LOADING_CONSTITUENTS,
41};
42pub use pole::solid_earth_pole_tide;
43
44use crate::astro::bodies::{sun_moon_ecef_with_polar_motion, SunMoonError};
45use crate::astro::constants::models::iers::SOLID_TIDE_EARTH_RADIUS_M;
46use crate::astro::constants::time::{
47    DAYS_PER_JULIAN_CENTURY, J2000_JD, SECONDS_PER_DAY, TT_MINUS_TAI_S,
48};
49use crate::astro::constants::units::{ARCSEC_TO_RAD, DEG_TO_RAD, KM_TO_M};
50use crate::astro::frames::transforms::{FrameTransformError, PolarMotion};
51use crate::astro::math::vec3::{dot3_ref as dot, norm3_ref as norm8};
52use crate::astro::time::{CoverageError, TimeScaleInputErrorKind, TimeScales};
53use crate::frame::{geodetic_to_itrf, ItrfPositionM, Wgs84Geodetic};
54use crate::validate::{self, FieldError};
55
56#[derive(Debug, Clone, Copy, PartialEq, Eq)]
57pub enum TideInputErrorKind {
58    Missing,
59    NonFinite,
60    NotPositive,
61    Negative,
62    OutOfRange,
63    FloatParse,
64    IntParse,
65    InvalidCivilDate,
66    InvalidCivilTime,
67}
68
69#[derive(Debug, Clone, PartialEq, Eq)]
70pub enum BlqParseErrorKind {
71    Empty,
72    MissingStation,
73    MissingCoefficientRows {
74        station: String,
75        expected: usize,
76        found: usize,
77    },
78    TooManyCoefficientRows {
79        station: String,
80    },
81    WrongColumnCount {
82        expected: usize,
83        found: usize,
84    },
85    InvalidNumber {
86        token: String,
87    },
88    NonFiniteNumber {
89        token: String,
90    },
91    UnsupportedConstituent {
92        constituent: String,
93    },
94    DuplicateConstituent {
95        constituent: String,
96    },
97    MultipleBlocks {
98        found: usize,
99    },
100}
101
102impl core::fmt::Display for BlqParseErrorKind {
103    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
104        match self {
105            Self::Empty => f.write_str("empty BLQ block"),
106            Self::MissingStation => f.write_str("missing station identifier"),
107            Self::MissingCoefficientRows {
108                station,
109                expected,
110                found,
111            } => write!(
112                f,
113                "station {station} has {found} coefficient rows, expected {expected}"
114            ),
115            Self::TooManyCoefficientRows { station } => {
116                write!(f, "station {station} has more than 6 coefficient rows")
117            }
118            Self::WrongColumnCount { expected, found } => {
119                write!(
120                    f,
121                    "coefficient row has {found} columns, expected {expected}"
122                )
123            }
124            Self::InvalidNumber { token } => write!(f, "invalid number {token:?}"),
125            Self::NonFiniteNumber { token } => write!(f, "non-finite number {token:?}"),
126            Self::UnsupportedConstituent { constituent } => {
127                write!(f, "unsupported constituent {constituent}")
128            }
129            Self::DuplicateConstituent { constituent } => {
130                write!(f, "duplicate constituent {constituent}")
131            }
132            Self::MultipleBlocks { found } => {
133                write!(f, "expected one BLQ station block, found {found}")
134            }
135        }
136    }
137}
138
139impl core::fmt::Display for TideInputErrorKind {
140    fn fmt(&self, f: &mut core::fmt::Formatter<'_>) -> core::fmt::Result {
141        f.write_str(match self {
142            Self::Missing => "missing",
143            Self::NonFinite => "not finite",
144            Self::NotPositive => "not positive",
145            Self::Negative => "negative",
146            Self::OutOfRange => "out of range",
147            Self::FloatParse => "invalid float",
148            Self::IntParse => "invalid integer",
149            Self::InvalidCivilDate => "invalid civil date",
150            Self::InvalidCivilTime => "invalid civil time",
151        })
152    }
153}
154
155impl From<&FieldError> for TideInputErrorKind {
156    fn from(error: &FieldError) -> Self {
157        match error {
158            FieldError::Missing { .. } => Self::Missing,
159            FieldError::NonFinite { .. } => Self::NonFinite,
160            FieldError::NotPositive { .. } => Self::NotPositive,
161            FieldError::Negative { .. } => Self::Negative,
162            FieldError::OutOfRange { .. } => Self::OutOfRange,
163            FieldError::FloatParse { .. } => Self::FloatParse,
164            FieldError::IntParse { .. } => Self::IntParse,
165            FieldError::InvalidCivilDate { .. } => Self::InvalidCivilDate,
166            FieldError::InvalidCivilTime { .. } => Self::InvalidCivilTime,
167        }
168    }
169}
170
171#[derive(Debug, Clone, PartialEq, Eq, thiserror::Error)]
172pub enum TideError {
173    #[error("invalid solid-earth tide input {field}: {kind}")]
174    InvalidInput {
175        field: &'static str,
176        kind: TideInputErrorKind,
177    },
178    #[error("station displacement time-scale conversion failed: {0}")]
179    TimeScale(#[from] CoverageError),
180    #[error("station displacement frame transform failed: {0}")]
181    FrameTransform(#[from] FrameTransformError),
182    #[error("station displacement Sun/Moon evaluation failed: {0}")]
183    SunMoon(#[from] SunMoonError),
184    #[error("missing station displacement input {field}")]
185    MissingInput { field: &'static str },
186    #[error("invalid BLQ block at line {line}: {kind}")]
187    BlqParse {
188        line: usize,
189        kind: BlqParseErrorKind,
190    },
191}
192
193fn invalid_tide_input(error: FieldError) -> TideError {
194    TideError::InvalidInput {
195        field: error.field(),
196        kind: (&error).into(),
197    }
198}
199
200fn map_time_input(error: CoverageError) -> TideError {
201    match error {
202        CoverageError::InvalidInput { field, kind } => TideError::InvalidInput {
203            field,
204            kind: tide_kind_from_time_kind(kind),
205        },
206        other => TideError::TimeScale(other),
207    }
208}
209
210fn tide_kind_from_time_kind(kind: TimeScaleInputErrorKind) -> TideInputErrorKind {
211    match kind {
212        TimeScaleInputErrorKind::Missing => TideInputErrorKind::Missing,
213        TimeScaleInputErrorKind::NonFinite => TideInputErrorKind::NonFinite,
214        TimeScaleInputErrorKind::NotPositive => TideInputErrorKind::NotPositive,
215        TimeScaleInputErrorKind::Negative => TideInputErrorKind::Negative,
216        TimeScaleInputErrorKind::OutOfRange => TideInputErrorKind::OutOfRange,
217        TimeScaleInputErrorKind::FloatParse => TideInputErrorKind::FloatParse,
218        TimeScaleInputErrorKind::IntParse => TideInputErrorKind::IntParse,
219        TimeScaleInputErrorKind::InvalidCivilDate => TideInputErrorKind::InvalidCivilDate,
220        TimeScaleInputErrorKind::InvalidCivilTime => TideInputErrorKind::InvalidCivilTime,
221    }
222}
223
224/// Station position accepted by the high-level displacement API.
225#[derive(Debug, Clone, Copy, PartialEq)]
226pub enum StationDisplacementPosition {
227    /// ITRF/ECEF metres.
228    Ecef(ItrfPositionM),
229    /// WGS84 geodetic radians/metres. Converted to ITRF through the public
230    /// frame conversion path before any tide model is evaluated.
231    Geodetic(Wgs84Geodetic),
232}
233
234impl From<ItrfPositionM> for StationDisplacementPosition {
235    fn from(value: ItrfPositionM) -> Self {
236        Self::Ecef(value)
237    }
238}
239
240impl From<Wgs84Geodetic> for StationDisplacementPosition {
241    fn from(value: Wgs84Geodetic) -> Self {
242        Self::Geodetic(value)
243    }
244}
245
246impl StationDisplacementPosition {
247    /// Construct from raw ITRF/ECEF metre components.
248    pub fn from_ecef_m(position_m: [f64; 3]) -> Result<Self, TideError> {
249        let position =
250            ItrfPositionM::new(position_m[0], position_m[1], position_m[2]).map_err(|error| {
251                match error {
252                    crate::frame::FrameValueError::InvalidInput { field, reason: _ } => {
253                        TideError::InvalidInput {
254                            field,
255                            kind: TideInputErrorKind::NonFinite,
256                        }
257                    }
258                }
259            })?;
260        Ok(Self::Ecef(position))
261    }
262
263    fn ecef_m(self) -> Result<[f64; 3], TideError> {
264        match self {
265            Self::Ecef(position) => Ok(position.as_array()),
266            Self::Geodetic(position) => Ok(geodetic_to_itrf(position)?.as_array()),
267        }
268    }
269}
270
271/// IERS polar-motion coordinates of the epoch, in arcseconds.
272#[derive(Debug, Clone, Copy, PartialEq)]
273pub struct StationPolarMotion {
274    pub xp_arcsec: f64,
275    pub yp_arcsec: f64,
276}
277
278impl StationPolarMotion {
279    pub const fn from_arcseconds(xp_arcsec: f64, yp_arcsec: f64) -> Self {
280        Self {
281            xp_arcsec,
282            yp_arcsec,
283        }
284    }
285
286    fn polar_motion(self) -> Result<PolarMotion, TideError> {
287        Ok(PolarMotion::from_radians(
288            self.xp_arcsec * ARCSEC_TO_RAD,
289            self.yp_arcsec * ARCSEC_TO_RAD,
290        )?)
291    }
292}
293
294/// UTC epoch for station displacement evaluation.
295#[derive(Debug, Clone, Copy, PartialEq)]
296pub struct StationDisplacementEpoch {
297    pub year: i32,
298    pub month: u8,
299    pub day: u8,
300    pub hour: u8,
301    pub minute: u8,
302    pub second: f64,
303    /// Optional IERS polar motion for pole tide and polar-motion-aware Sun/Moon
304    /// rotation.
305    pub polar_motion: Option<StationPolarMotion>,
306}
307
308impl StationDisplacementEpoch {
309    pub const fn from_utc(
310        year: i32,
311        month: u8,
312        day: u8,
313        hour: u8,
314        minute: u8,
315        second: f64,
316    ) -> Self {
317        Self {
318            year,
319            month,
320            day,
321            hour,
322            minute,
323            second,
324            polar_motion: None,
325        }
326    }
327
328    pub const fn with_polar_motion_arcsec(mut self, xp_arcsec: f64, yp_arcsec: f64) -> Self {
329        self.polar_motion = Some(StationPolarMotion::from_arcseconds(xp_arcsec, yp_arcsec));
330        self
331    }
332
333    fn time_scales(self) -> Result<TimeScales, TideError> {
334        TimeScales::from_utc(
335            self.year,
336            i32::from(self.month),
337            i32::from(self.day),
338            i32::from(self.hour),
339            i32::from(self.minute),
340            self.second,
341        )
342        .map_err(map_time_input)
343    }
344
345    fn validate_utc(self) -> Result<(), TideError> {
346        validate::civil_datetime_with_second_policy(
347            i64::from(self.year),
348            i64::from(self.month),
349            i64::from(self.day),
350            i64::from(self.hour),
351            i64::from(self.minute),
352            self.second,
353            validate::CivilSecondPolicy::Continuous,
354        )
355        .map(|_| ())
356        .map_err(invalid_tide_input)
357    }
358
359    fn fractional_hour(self) -> f64 {
360        f64::from(self.hour) + f64::from(self.minute) / 60.0 + self.second / 3600.0
361    }
362}
363
364/// Switches for the high-level station displacement entry.
365#[derive(Debug, Clone, Copy, PartialEq)]
366pub struct StationDisplacementOptions<'a> {
367    /// Apply the IERS solid Earth tide station displacement.
368    pub solid_earth_tide: bool,
369    /// Apply the IERS pole tide station displacement. Each epoch must carry
370    /// polar motion when this is true.
371    pub pole_tide: bool,
372    /// Optional ocean-loading BLQ coefficients supplied by the caller.
373    pub ocean_loading: Option<&'a OceanLoadingBlq>,
374}
375
376impl Default for StationDisplacementOptions<'_> {
377    fn default() -> Self {
378        Self {
379            solid_earth_tide: true,
380            pole_tide: false,
381            ocean_loading: None,
382        }
383    }
384}
385
386/// Component-resolved station displacement in ITRF/ECEF metres.
387#[derive(Debug, Clone, Copy, PartialEq)]
388pub struct StationDisplacement {
389    /// Sum of all enabled component displacements, in ITRF/ECEF metres.
390    pub ecef_m: [f64; 3],
391    pub solid_earth_tide_ecef_m: Option<[f64; 3]>,
392    pub pole_tide_ecef_m: Option<[f64; 3]>,
393    pub ocean_loading_ecef_m: Option<[f64; 3]>,
394}
395
396impl StationDisplacement {
397    fn zero() -> Self {
398        Self {
399            ecef_m: [0.0; 3],
400            solid_earth_tide_ecef_m: None,
401            pole_tide_ecef_m: None,
402            ocean_loading_ecef_m: None,
403        }
404    }
405
406    fn add_component(total: &mut [f64; 3], component: [f64; 3]) {
407        for i in 0..3 {
408            total[i] += component[i];
409        }
410    }
411}
412
413/// Evaluate the enabled station displacement corrections, returning ITRF/ECEF
414/// metre components.
415///
416/// The solid Earth tide path uses IERS Conventions (2010), Chapter 7 station
417/// displacement with the permanent tide retained. The low-level
418/// [`solid_earth_tide`] routine ships the in-phase degree-2 and degree-3
419/// displacement, the step-1 out-of-phase and latitude-dependence corrections,
420/// and the step-2 diurnal/long-period frequency corrections; it leaves the
421/// optional step-3 permanent-tide removal disabled for ITRF/IGS use. Sun/Moon
422/// positions are generated through the same Earth-fixed analytic ephemeris path
423/// used by the tide-force lane, including caller-supplied polar motion when the
424/// epoch carries it.
425pub fn station_displacement_ecef_m(
426    position: StationDisplacementPosition,
427    epoch: StationDisplacementEpoch,
428    options: StationDisplacementOptions<'_>,
429) -> Result<StationDisplacement, TideError> {
430    let receiver_ecef_m = position.ecef_m()?;
431    epoch.validate_utc()?;
432    let fhr = epoch.fractional_hour();
433    let mut displacement = StationDisplacement::zero();
434
435    if options.solid_earth_tide {
436        let ts = epoch.time_scales()?;
437        let polar_motion = epoch
438            .polar_motion
439            .map(StationPolarMotion::polar_motion)
440            .transpose()?
441            .unwrap_or_default();
442        let sun_moon = sun_moon_ecef_with_polar_motion(&ts, polar_motion)?;
443        let solid = solid_earth_tide(
444            &receiver_ecef_m,
445            epoch.year,
446            i32::from(epoch.month),
447            i32::from(epoch.day),
448            fhr,
449            &sun_moon.sun,
450            &sun_moon.moon,
451        )?;
452        StationDisplacement::add_component(&mut displacement.ecef_m, solid);
453        displacement.solid_earth_tide_ecef_m = Some(solid);
454    }
455
456    if options.pole_tide {
457        let polar = epoch.polar_motion.ok_or(TideError::MissingInput {
458            field: "polar motion",
459        })?;
460        let pole = solid_earth_pole_tide(
461            &receiver_ecef_m,
462            epoch.year,
463            i32::from(epoch.month),
464            i32::from(epoch.day),
465            fhr,
466            polar.xp_arcsec,
467            polar.yp_arcsec,
468        )?;
469        StationDisplacement::add_component(&mut displacement.ecef_m, pole);
470        displacement.pole_tide_ecef_m = Some(pole);
471    }
472
473    if let Some(blq) = options.ocean_loading {
474        let ocean = ocean_tide_loading(
475            &receiver_ecef_m,
476            epoch.year,
477            i32::from(epoch.month),
478            i32::from(epoch.day),
479            fhr,
480            blq,
481        )?;
482        StationDisplacement::add_component(&mut displacement.ecef_m, ocean);
483        displacement.ocean_loading_ecef_m = Some(ocean);
484    }
485
486    Ok(displacement)
487}
488
489/// Evaluate station displacement for many epochs. Each element is equivalent to
490/// a scalar [`station_displacement_ecef_m`] call for the same position, epoch,
491/// and options, so per-epoch failures stay local to their output row.
492pub fn station_displacement_ecef_m_batch(
493    position: StationDisplacementPosition,
494    epochs: &[StationDisplacementEpoch],
495    options: StationDisplacementOptions<'_>,
496) -> Vec<Result<StationDisplacement, TideError>> {
497    epochs
498        .iter()
499        .map(|&epoch| station_displacement_ecef_m(position, epoch, options))
500        .collect()
501}
502
503/// Solid-earth tide displacement of an ITRF station, in metres (ECEF).
504///
505/// Arguments mirror the IERS reference routine:
506/// * `xsta` - geocentric station position (m, ITRF).
507/// * `year`, `month`, `day` - UTC calendar date.
508/// * `fhr` - UTC fractional hour of the day (hour + min/60 + sec/3600).
509/// * `xsun` - geocentric Sun position (m, ECEF).
510/// * `xmon` - geocentric Moon position (m, ECEF).
511///
512/// Returns the displacement vector `dxtide` (m, geocentric ITRF). The permanent
513/// (mean) tide deformation is retained (ITRF/IGS convention).
514///
515/// Returns [`TideError`] when inputs are non-finite or geometrically
516/// degenerate: the station vector must be non-zero and non-polar, and Sun/Moon
517/// vectors must be non-zero.
518pub fn solid_earth_tide(
519    xsta: &[f64; 3],
520    year: i32,
521    month: i32,
522    day: i32,
523    fhr: f64,
524    xsun: &[f64; 3],
525    xmon: &[f64; 3],
526) -> Result<[f64; 3], TideError> {
527    validate_tide_domain(xsta, year, month, day, fhr, xsun, xmon)?;
528    Ok(solid_earth_tide_unchecked(
529        xsta, year, month, day, fhr, xsun, xmon,
530    ))
531}
532
533fn validate_tide_domain(
534    xsta: &[f64; 3],
535    year: i32,
536    month: i32,
537    day: i32,
538    fhr: f64,
539    xsun: &[f64; 3],
540    xmon: &[f64; 3],
541) -> Result<(), TideError> {
542    validate::finite_vec3(*xsta, "station position").map_err(invalid_tide_input)?;
543    validate::civil_datetime_with_second_policy(
544        i64::from(year),
545        i64::from(month),
546        i64::from(day),
547        0,
548        0,
549        0.0,
550        validate::CivilSecondPolicy::Continuous,
551    )
552    .map_err(invalid_tide_input)?;
553    validate::finite_in_range_exclusive_upper(fhr, 0.0, 24.0, "fractional hour")
554        .map_err(invalid_tide_input)?;
555    validate::finite_vec3(*xsun, "sun position").map_err(invalid_tide_input)?;
556    validate::finite_vec3(*xmon, "moon position").map_err(invalid_tide_input)?;
557
558    validate::finite_positive(norm8(xsta), "station radius").map_err(invalid_tide_input)?;
559    let station_horizontal_radius = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt();
560    validate::finite_positive(station_horizontal_radius, "station horizontal radius")
561        .map_err(invalid_tide_input)?;
562    validate::finite_positive(norm8(xsun), "sun radius").map_err(invalid_tide_input)?;
563    validate::finite_positive(norm8(xmon), "moon radius").map_err(invalid_tide_input)?;
564
565    Ok(())
566}
567
568fn solid_earth_tide_unchecked(
569    xsta: &[f64; 3],
570    year: i32,
571    month: i32,
572    day: i32,
573    fhr: f64,
574    xsun: &[f64; 3],
575    xmon: &[f64; 3],
576) -> [f64; 3] {
577    // Nominal second- and third-degree Love and Shida numbers.
578    const H20: f64 = 0.6078;
579    const L20: f64 = 0.0847;
580    const H3: f64 = 0.292;
581    const L3: f64 = 0.015;
582
583    // Scalar product of station vector with Sun/Moon vector (SPROD).
584    let rsta = norm8(xsta);
585    let rsun = norm8(xsun);
586    let rmon = norm8(xmon);
587    let scs = dot(xsta, xsun);
588    let scm = dot(xsta, xmon);
589    let scsun = scs / rsta / rsun;
590    let scmon = scm / rsta / rmon;
591
592    // Latitude-corrected H2 and L2.
593    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
594    let h2 = H20 - 0.0006 * (1.0 - 3.0 / 2.0 * cosphi * cosphi);
595    let l2 = L20 + 0.0002 * (1.0 - 3.0 / 2.0 * cosphi * cosphi);
596
597    // P2 term.
598    let p2sun = 3.0 * (h2 / 2.0 - l2) * scsun * scsun - h2 / 2.0;
599    let p2mon = 3.0 * (h2 / 2.0 - l2) * scmon * scmon - h2 / 2.0;
600
601    // P3 term.
602    let p3sun = 5.0 / 2.0 * (H3 - 3.0 * L3) * scsun.powi(3) + 3.0 / 2.0 * (L3 - H3) * scsun;
603    let p3mon = 5.0 / 2.0 * (H3 - 3.0 * L3) * scmon.powi(3) + 3.0 / 2.0 * (L3 - H3) * scmon;
604
605    // Term in direction of Sun/Moon vector.
606    let x2sun = 3.0 * l2 * scsun;
607    let x2mon = 3.0 * l2 * scmon;
608    let x3sun = 3.0 * L3 / 2.0 * (5.0 * scsun * scsun - 1.0);
609    let x3mon = 3.0 * L3 / 2.0 * (5.0 * scmon * scmon - 1.0);
610
611    // Factors for Sun/Moon (IAU current best estimates).
612    const MASS_RATIO_SUN: f64 = 332946.0482;
613    const MASS_RATIO_MOON: f64 = 0.0123000371;
614    const RE: f64 = SOLID_TIDE_EARTH_RADIUS_M;
615    let fac2sun = MASS_RATIO_SUN * RE * (RE / rsun).powi(3);
616    let fac2mon = MASS_RATIO_MOON * RE * (RE / rmon).powi(3);
617    let fac3sun = fac2sun * (RE / rsun);
618    let fac3mon = fac2mon * (RE / rmon);
619
620    // Total in-phase degree-2/degree-3 displacement.
621    let mut dxtide = [0.0_f64; 3];
622    for i in 0..3 {
623        dxtide[i] = fac2sun * (x2sun * xsun[i] / rsun + p2sun * xsta[i] / rsta)
624            + fac2mon * (x2mon * xmon[i] / rmon + p2mon * xsta[i] / rsta)
625            + fac3sun * (x3sun * xsun[i] / rsun + p3sun * xsta[i] / rsta)
626            + fac3mon * (x3mon * xmon[i] / rmon + p3mon * xsta[i] / rsta);
627    }
628
629    // Out-of-phase corrections (diurnal, semi-diurnal) and latitude dependence.
630    let c = st1idiu(xsta, xsun, xmon, fac2sun, fac2mon);
631    for i in 0..3 {
632        dxtide[i] += c[i];
633    }
634    let c = st1isem(xsta, xsun, xmon, fac2sun, fac2mon);
635    for i in 0..3 {
636        dxtide[i] += c[i];
637    }
638    let c = st1l1(xsta, xsun, xmon, fac2sun, fac2mon);
639    for i in 0..3 {
640        dxtide[i] += c[i];
641    }
642
643    // Step 2 corrections need the date in Julian centuries (TT).
644    let (jjm0, jjm1) = cal2jd(year, month, day);
645    let fhrd = fhr / 24.0;
646    let mut t = ((jjm0 - J2000_JD) + jjm1 + fhrd) / DAYS_PER_JULIAN_CENTURY;
647    let dtt = dat(year, month, day) + TT_MINUS_TAI_S;
648    t += dtt / (SECONDS_PER_DAY * DAYS_PER_JULIAN_CENTURY);
649
650    let c = step2diu(xsta, fhr, t);
651    for i in 0..3 {
652        dxtide[i] += c[i];
653    }
654    let c = step2lon(xsta, t);
655    for i in 0..3 {
656        dxtide[i] += c[i];
657    }
658
659    // Step 3 of the IERS routine, the permanent (zero-frequency) tide removal,
660    // is intentionally not applied, so the permanent (mean) tide deformation is
661    // retained (the ITRF/IGS conform-to-mean-tide convention; see module docs).
662    dxtide
663}
664
665/// Out-of-phase part of the Love numbers, diurnal band (ST1IDIU).
666fn st1idiu(
667    xsta: &[f64; 3],
668    xsun: &[f64; 3],
669    xmon: &[f64; 3],
670    fac2sun: f64,
671    fac2mon: f64,
672) -> [f64; 3] {
673    const DHI: f64 = -0.0025;
674    const DLI: f64 = -0.0007;
675    let rsta = norm8(xsta);
676    let sinphi = xsta[2] / rsta;
677    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
678    let cos2phi = cosphi * cosphi - sinphi * sinphi;
679    let sinla = xsta[1] / cosphi / rsta;
680    let cosla = xsta[0] / cosphi / rsta;
681    let rmon = norm8(xmon);
682    let rsun = norm8(xsun);
683
684    let drsun =
685        -3.0 * DHI * sinphi * cosphi * fac2sun * xsun[2] * (xsun[0] * sinla - xsun[1] * cosla)
686            / (rsun * rsun);
687    let drmon =
688        -3.0 * DHI * sinphi * cosphi * fac2mon * xmon[2] * (xmon[0] * sinla - xmon[1] * cosla)
689            / (rmon * rmon);
690    let dnsun = -3.0 * DLI * cos2phi * fac2sun * xsun[2] * (xsun[0] * sinla - xsun[1] * cosla)
691        / (rsun * rsun);
692    let dnmon = -3.0 * DLI * cos2phi * fac2mon * xmon[2] * (xmon[0] * sinla - xmon[1] * cosla)
693        / (rmon * rmon);
694    let desun = -3.0 * DLI * sinphi * fac2sun * xsun[2] * (xsun[0] * cosla + xsun[1] * sinla)
695        / (rsun * rsun);
696    let demon = -3.0 * DLI * sinphi * fac2mon * xmon[2] * (xmon[0] * cosla + xmon[1] * sinla)
697        / (rmon * rmon);
698
699    let dr = drsun + drmon;
700    let dn = dnsun + dnmon;
701    let de = desun + demon;
702
703    [
704        dr * cosla * cosphi - de * sinla - dn * sinphi * cosla,
705        dr * sinla * cosphi + de * cosla - dn * sinphi * sinla,
706        dr * sinphi + dn * cosphi,
707    ]
708}
709
710/// Out-of-phase part of the Love numbers, semi-diurnal band (ST1ISEM).
711fn st1isem(
712    xsta: &[f64; 3],
713    xsun: &[f64; 3],
714    xmon: &[f64; 3],
715    fac2sun: f64,
716    fac2mon: f64,
717) -> [f64; 3] {
718    const DHI: f64 = -0.0022;
719    const DLI: f64 = -0.0007;
720    let rsta = norm8(xsta);
721    let sinphi = xsta[2] / rsta;
722    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
723    let sinla = xsta[1] / cosphi / rsta;
724    let cosla = xsta[0] / cosphi / rsta;
725    let costwola = cosla * cosla - sinla * sinla;
726    let sintwola = 2.0 * cosla * sinla;
727    let rmon = norm8(xmon);
728    let rsun = norm8(xsun);
729
730    let drsun = -3.0 / 4.0
731        * DHI
732        * cosphi
733        * cosphi
734        * fac2sun
735        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * sintwola - 2.0 * xsun[0] * xsun[1] * costwola)
736        / (rsun * rsun);
737    let drmon = -3.0 / 4.0
738        * DHI
739        * cosphi
740        * cosphi
741        * fac2mon
742        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * sintwola - 2.0 * xmon[0] * xmon[1] * costwola)
743        / (rmon * rmon);
744    let dnsun = 3.0 / 2.0
745        * DLI
746        * sinphi
747        * cosphi
748        * fac2sun
749        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * sintwola - 2.0 * xsun[0] * xsun[1] * costwola)
750        / (rsun * rsun);
751    let dnmon = 3.0 / 2.0
752        * DLI
753        * sinphi
754        * cosphi
755        * fac2mon
756        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * sintwola - 2.0 * xmon[0] * xmon[1] * costwola)
757        / (rmon * rmon);
758    let desun = -3.0 / 2.0
759        * DLI
760        * cosphi
761        * fac2sun
762        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * costwola + 2.0 * xsun[0] * xsun[1] * sintwola)
763        / (rsun * rsun);
764    let demon = -3.0 / 2.0
765        * DLI
766        * cosphi
767        * fac2mon
768        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * costwola + 2.0 * xmon[0] * xmon[1] * sintwola)
769        / (rmon * rmon);
770
771    let dr = drsun + drmon;
772    let dn = dnsun + dnmon;
773    let de = desun + demon;
774
775    [
776        dr * cosla * cosphi - de * sinla - dn * sinphi * cosla,
777        dr * sinla * cosphi + de * cosla - dn * sinphi * sinla,
778        dr * sinphi + dn * cosphi,
779    ]
780}
781
782/// Latitude dependence of the Love numbers, part L^(1) (ST1L1).
783fn st1l1(
784    xsta: &[f64; 3],
785    xsun: &[f64; 3],
786    xmon: &[f64; 3],
787    fac2sun: f64,
788    fac2mon: f64,
789) -> [f64; 3] {
790    const L1D: f64 = 0.0012;
791    const L1SD: f64 = 0.0024;
792    let rsta = norm8(xsta);
793    let sinphi = xsta[2] / rsta;
794    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
795    let sinla = xsta[1] / cosphi / rsta;
796    let cosla = xsta[0] / cosphi / rsta;
797    let rmon = norm8(xmon);
798    let rsun = norm8(xsun);
799
800    // Diurnal band.
801    let mut l1 = L1D;
802    let dnsun = -l1 * sinphi * sinphi * fac2sun * xsun[2] * (xsun[0] * cosla + xsun[1] * sinla)
803        / (rsun * rsun);
804    let dnmon = -l1 * sinphi * sinphi * fac2mon * xmon[2] * (xmon[0] * cosla + xmon[1] * sinla)
805        / (rmon * rmon);
806    let desun = l1
807        * sinphi
808        * (cosphi * cosphi - sinphi * sinphi)
809        * fac2sun
810        * xsun[2]
811        * (xsun[0] * sinla - xsun[1] * cosla)
812        / (rsun * rsun);
813    let demon = l1
814        * sinphi
815        * (cosphi * cosphi - sinphi * sinphi)
816        * fac2mon
817        * xmon[2]
818        * (xmon[0] * sinla - xmon[1] * cosla)
819        / (rmon * rmon);
820
821    let de = 3.0 * (desun + demon);
822    let dn = 3.0 * (dnsun + dnmon);
823
824    let mut xcorsta = [
825        -de * sinla - dn * sinphi * cosla,
826        de * cosla - dn * sinphi * sinla,
827        dn * cosphi,
828    ];
829
830    // Semi-diurnal band.
831    l1 = L1SD;
832    let costwola = cosla * cosla - sinla * sinla;
833    let sintwola = 2.0 * cosla * sinla;
834
835    let dnsun = -l1 / 2.0
836        * sinphi
837        * cosphi
838        * fac2sun
839        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * costwola + 2.0 * xsun[0] * xsun[1] * sintwola)
840        / (rsun * rsun);
841    let dnmon = -l1 / 2.0
842        * sinphi
843        * cosphi
844        * fac2mon
845        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * costwola + 2.0 * xmon[0] * xmon[1] * sintwola)
846        / (rmon * rmon);
847    let desun = -l1 / 2.0
848        * sinphi
849        * sinphi
850        * cosphi
851        * fac2sun
852        * ((xsun[0] * xsun[0] - xsun[1] * xsun[1]) * sintwola - 2.0 * xsun[0] * xsun[1] * costwola)
853        / (rsun * rsun);
854    let demon = -l1 / 2.0
855        * sinphi
856        * sinphi
857        * cosphi
858        * fac2mon
859        * ((xmon[0] * xmon[0] - xmon[1] * xmon[1]) * sintwola - 2.0 * xmon[0] * xmon[1] * costwola)
860        / (rmon * rmon);
861
862    let de = 3.0 * (desun + demon);
863    let dn = 3.0 * (dnsun + dnmon);
864
865    xcorsta[0] += -de * sinla - dn * sinphi * cosla;
866    xcorsta[1] += de * cosla - dn * sinphi * sinla;
867    xcorsta[2] += dn * cosphi;
868    xcorsta
869}
870
871/// In-phase / out-of-phase frequency-dependent corrections, diurnal band
872/// (STEP2DIU). `fhr` is UTC fractional hour, `t` is Julian centuries (TT).
873fn step2diu(xsta: &[f64; 3], fhr: f64, t: f64) -> [f64; 3] {
874    // DATDI(9,31): {l, l', F, D, Omega(Ps), Adr, Adi, Anr, Ani} per wave.
875    #[rustfmt::skip]
876    const DATDI: [[f64; 9]; 31] = [
877        [-3.0, 0.0, 2.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
878        [-3.0, 2.0, 0.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
879        [-2.0, 0.0, 1.0, -1.0, 0.0, -0.02, 0.0, 0.0, 0.0],
880        [-2.0, 0.0, 1.0, 0.0, 0.0, -0.08, 0.0, -0.01, 0.01],
881        [-2.0, 2.0, -1.0, 0.0, 0.0, -0.02, 0.0, 0.0, 0.0],
882        [-1.0, 0.0, 0.0, -1.0, 0.0, -0.10, 0.0, 0.0, 0.0],
883        [-1.0, 0.0, 0.0, 0.0, 0.0, -0.51, 0.0, -0.02, 0.03],
884        [-1.0, 2.0, 0.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0],
885        [0.0, -2.0, 1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0],
886        [0.0, 0.0, -1.0, 0.0, 0.0, 0.02, 0.0, 0.0, 0.0],
887        [0.0, 0.0, 1.0, 0.0, 0.0, 0.06, 0.0, 0.0, 0.0],
888        [0.0, 0.0, 1.0, 1.0, 0.0, 0.01, 0.0, 0.0, 0.0],
889        [0.0, 2.0, -1.0, 0.0, 0.0, 0.01, 0.0, 0.0, 0.0],
890        [1.0, -3.0, 0.0, 0.0, 1.0, -0.06, 0.0, 0.0, 0.0],
891        [1.0, -2.0, 0.0, -1.0, 0.0, 0.01, 0.0, 0.0, 0.0],
892        [1.0, -2.0, 0.0, 0.0, 0.0, -1.23, -0.07, 0.06, 0.01],
893        [1.0, -1.0, 0.0, 0.0, -1.0, 0.02, 0.0, 0.0, 0.0],
894        [1.0, -1.0, 0.0, 0.0, 1.0, 0.04, 0.0, 0.0, 0.0],
895        [1.0, 0.0, 0.0, -1.0, 0.0, -0.22, 0.01, 0.01, 0.0],
896        [1.0, 0.0, 0.0, 0.0, 0.0, 12.00, -0.80, -0.67, -0.03],
897        [1.0, 0.0, 0.0, 1.0, 0.0, 1.73, -0.12, -0.10, 0.0],
898        [1.0, 0.0, 0.0, 2.0, 0.0, -0.04, 0.0, 0.0, 0.0],
899        [1.0, 1.0, 0.0, 0.0, -1.0, -0.50, -0.01, 0.03, 0.0],
900        [1.0, 1.0, 0.0, 0.0, 1.0, 0.01, 0.0, 0.0, 0.0],
901        [0.0, 1.0, 0.0, 1.0, -1.0, -0.01, 0.0, 0.0, 0.0],
902        [1.0, 2.0, -2.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
903        [1.0, 2.0, 0.0, 0.0, 0.0, -0.11, 0.01, 0.01, 0.0],
904        [2.0, -2.0, 1.0, 0.0, 0.0, -0.01, 0.0, 0.0, 0.0],
905        [2.0, 0.0, -1.0, 0.0, 0.0, -0.02, 0.0, 0.0, 0.0],
906        [3.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0],
907        [3.0, 0.0, 0.0, 1.0, 0.0, 0.0, 0.0, 0.0, 0.0],
908    ];
909    let mut s = 218.31664563 + (481267.88194 + (-0.0014663889 + 0.00000185139 * t) * t) * t;
910    let mut tau = fhr * 15.0
911        + 280.4606184
912        + (36000.7700536 + (0.00038793 + -0.0000000258 * t) * t) * t
913        + (-s);
914    let pr = (1.396971278 + (0.000308889 + (0.000000021 + 0.000000007 * t) * t) * t) * t;
915    s += pr;
916    let mut h = 280.46645
917        + (36000.7697489 + (0.00030322222 + (0.000000020 + -0.00000000654 * t) * t) * t) * t;
918    let mut p = 83.35324312
919        + (4069.01363525 + (-0.01032172222 + (-0.0000124991 + 0.00000005263 * t) * t) * t) * t;
920    let mut zns = 234.95544499
921        + (1934.13626197 + (-0.00207561111 + (-0.00000213944 + 0.00000001650 * t) * t) * t) * t;
922    let mut ps = 282.93734098
923        + (1.71945766667 + (0.00045688889 + (-0.00000001778 + -0.00000000334 * t) * t) * t) * t;
924
925    s = s.rem_euclid(360.0);
926    tau = tau.rem_euclid(360.0);
927    h = h.rem_euclid(360.0);
928    p = p.rem_euclid(360.0);
929    zns = zns.rem_euclid(360.0);
930    ps = ps.rem_euclid(360.0);
931
932    let rsta = (xsta[0] * xsta[0] + xsta[1] * xsta[1] + xsta[2] * xsta[2]).sqrt();
933    let sinphi = xsta[2] / rsta;
934    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
935    let cosla = xsta[0] / cosphi / rsta;
936    let sinla = xsta[1] / cosphi / rsta;
937    let zla = xsta[1].atan2(xsta[0]);
938
939    let mut xcorsta = [0.0_f64; 3];
940    for w in &DATDI {
941        let thetaf = (tau + w[0] * s + w[1] * h + w[2] * p + w[3] * zns + w[4] * ps) * DEG_TO_RAD;
942        let dr = w[5] * 2.0 * sinphi * cosphi * (thetaf + zla).sin()
943            + w[6] * 2.0 * sinphi * cosphi * (thetaf + zla).cos();
944        let dn = w[7] * (cosphi * cosphi - sinphi * sinphi) * (thetaf + zla).sin()
945            + w[8] * (cosphi * cosphi - sinphi * sinphi) * (thetaf + zla).cos();
946        let de = w[7] * sinphi * (thetaf + zla).cos() - w[8] * sinphi * (thetaf + zla).sin();
947
948        xcorsta[0] += dr * cosla * cosphi - de * sinla - dn * sinphi * cosla;
949        xcorsta[1] += dr * sinla * cosphi + de * cosla - dn * sinphi * sinla;
950        xcorsta[2] += dr * sinphi + dn * cosphi;
951    }
952    for v in &mut xcorsta {
953        *v /= KM_TO_M;
954    }
955    xcorsta
956}
957
958/// In-phase / out-of-phase frequency-dependent corrections, long-period band
959/// (STEP2LON). `t` is Julian centuries (TT).
960fn step2lon(xsta: &[f64; 3], t: f64) -> [f64; 3] {
961    #[rustfmt::skip]
962    const DATDI: [[f64; 9]; 5] = [
963        [0.0, 0.0, 0.0, 1.0, 0.0, 0.47, 0.23, 0.16, 0.07],
964        [0.0, 2.0, 0.0, 0.0, 0.0, -0.20, -0.12, -0.11, -0.05],
965        [1.0, 0.0, -1.0, 0.0, 0.0, -0.11, -0.08, -0.09, -0.04],
966        [2.0, 0.0, 0.0, 0.0, 0.0, -0.13, -0.11, -0.15, -0.07],
967        [2.0, 0.0, 0.0, 1.0, 0.0, -0.05, -0.05, -0.06, -0.03],
968    ];
969    let mut s = 218.31664563 + (481267.88194 + (-0.0014663889 + 0.00000185139 * t) * t) * t;
970    let pr = (1.396971278 + (0.000308889 + (0.000000021 + 0.000000007 * t) * t) * t) * t;
971    s += pr;
972    let mut h = 280.46645
973        + (36000.7697489 + (0.00030322222 + (0.000000020 + -0.00000000654 * t) * t) * t) * t;
974    let mut p = 83.35324312
975        + (4069.01363525 + (-0.01032172222 + (-0.0000124991 + 0.00000005263 * t) * t) * t) * t;
976    let mut zns = 234.95544499
977        + (1934.13626197 + (-0.00207561111 + (-0.00000213944 + 0.00000001650 * t) * t) * t) * t;
978    let mut ps = 282.93734098
979        + (1.71945766667 + (0.00045688889 + (-0.00000001778 + -0.00000000334 * t) * t) * t) * t;
980
981    let rsta = (xsta[0] * xsta[0] + xsta[1] * xsta[1] + xsta[2] * xsta[2]).sqrt();
982    let sinphi = xsta[2] / rsta;
983    let cosphi = (xsta[0] * xsta[0] + xsta[1] * xsta[1]).sqrt() / rsta;
984    let cosla = xsta[0] / cosphi / rsta;
985    let sinla = xsta[1] / cosphi / rsta;
986
987    s = s.rem_euclid(360.0);
988    h = h.rem_euclid(360.0);
989    p = p.rem_euclid(360.0);
990    zns = zns.rem_euclid(360.0);
991    ps = ps.rem_euclid(360.0);
992
993    let mut xcorsta = [0.0_f64; 3];
994    for w in &DATDI {
995        let thetaf = (w[0] * s + w[1] * h + w[2] * p + w[3] * zns + w[4] * ps) * DEG_TO_RAD;
996        let dr = w[5] * (3.0 * sinphi * sinphi - 1.0) / 2.0 * thetaf.cos()
997            + w[7] * (3.0 * sinphi * sinphi - 1.0) / 2.0 * thetaf.sin();
998        let dn = w[6] * (cosphi * sinphi * 2.0) * thetaf.cos()
999            + w[8] * (cosphi * sinphi * 2.0) * thetaf.sin();
1000        let de = 0.0;
1001
1002        xcorsta[0] += dr * cosla * cosphi - de * sinla - dn * sinphi * cosla;
1003        xcorsta[1] += dr * sinla * cosphi + de * cosla - dn * sinphi * sinla;
1004        xcorsta[2] += dr * sinphi + dn * cosphi;
1005    }
1006    for v in &mut xcorsta {
1007        *v /= KM_TO_M;
1008    }
1009    xcorsta
1010}
1011
1012/// Gregorian calendar date -> (MJD epoch 2400000.5, MJD) (SOFA CAL2JD).
1013///
1014/// This is a SOFA parity adapter, deliberately NOT routed through
1015/// [`crate::astro::time::civil`]: the solid-Earth/ocean/pole tide models are
1016/// validated bit-for-bit against the IERS/SOFA reference (the
1017/// `ocean_loading_oracle` test), so the calendar-to-MJD step must reproduce
1018/// SOFA's `iauCal2jd` exactly. It is kept local under this tides-specific name
1019/// so it is not mistaken for a duplicate of the canonical civil conversions and
1020/// is never consolidated into them.
1021fn cal2jd(iy: i32, im: i32, id: i32) -> (f64, f64) {
1022    let my = (im - 14) / 12;
1023    let iypmy = iy + my;
1024    let djm0 = 2400000.5;
1025    let djm = ((1461 * (iypmy + 4800)) / 4 + (367 * (im - 2 - 12 * my)) / 12
1026        - (3 * ((iypmy + 4900) / 100)) / 4
1027        + id
1028        - 2432076) as f64;
1029    (djm0, djm)
1030}
1031
1032/// TAI-UTC (Delta(AT)) in seconds for the given date (SOFA DAT, leap-second
1033/// table only; the four golden dates are all post-1972 so the pre-1972 drift
1034/// branch is not exercised, but it is retained for completeness).
1035fn dat(iy: i32, im: i32, _id: i32) -> f64 {
1036    // Post-1972 leap-second table: (year, month, Delta(AT) seconds).
1037    const IDAT: [(i32, i32, f64); 28] = [
1038        (1972, 1, 10.0),
1039        (1972, 7, 11.0),
1040        (1973, 1, 12.0),
1041        (1974, 1, 13.0),
1042        (1975, 1, 14.0),
1043        (1976, 1, 15.0),
1044        (1977, 1, 16.0),
1045        (1978, 1, 17.0),
1046        (1979, 1, 18.0),
1047        (1980, 1, 19.0),
1048        (1981, 7, 20.0),
1049        (1982, 7, 21.0),
1050        (1983, 7, 22.0),
1051        (1985, 7, 23.0),
1052        (1988, 1, 24.0),
1053        (1990, 1, 25.0),
1054        (1991, 1, 26.0),
1055        (1992, 7, 27.0),
1056        (1993, 7, 28.0),
1057        (1994, 7, 29.0),
1058        (1996, 1, 30.0),
1059        (1997, 7, 31.0),
1060        (1999, 1, 32.0),
1061        (2006, 1, 33.0),
1062        (2009, 1, 34.0),
1063        (2012, 7, 35.0),
1064        (2015, 7, 36.0),
1065        (2017, 1, 37.0),
1066    ];
1067    let m = 12 * iy + im;
1068    let mut da = IDAT[0].2;
1069    for &(y, mo, d) in &IDAT {
1070        if m >= 12 * y + mo {
1071            da = d;
1072        }
1073    }
1074    da
1075}
sidereon_core/tides/mod.rs

sidereon_core/tides/
mod.rs
