553 lines
21 KiB
Rust
553 lines
21 KiB
Rust
//! Tests for the IAU 2000/2006 Earth orientation models.
|
||
//!
|
||
//! The expected values are those of the IAU SOFA library and its ERFA
|
||
//! derivative for the same instants. They are recorded here so that a change
|
||
//! to a coefficient or a rotation convention shows up as a failing
|
||
//! assertion.
|
||
//!
|
||
//! IAU 2000/2006 地球定向模型的测试。
|
||
//!
|
||
//! 期望值是 IAU SOFA 库及其 ERFA 衍生物对相同瞬时的值。它们被记录在此，使对某个
|
||
//! 系数或转动约定的改动表现为断言失败。
|
||
|
||
use apsis_frames::iau2000::{
|
||
PolarMotion, bias_precession_nutation_matrix, celestial_to_terrestrial_matrix,
|
||
equation_of_the_equinoxes, equation_of_the_equinoxes_complementary_terms, general_precession,
|
||
greenwich_apparent_sidereal_time, mean_obliquity, nutation_2000a,
|
||
nutation_2000a_with_2006_adjustment, precession_angles,
|
||
};
|
||
use core::f64::consts::PI;
|
||
|
||
/// Converts arcseconds to radians.
|
||
///
|
||
/// 把角秒换算为弧度。
|
||
const ARCSECOND: f64 = PI / (180.0 * 3_600.0);
|
||
|
||
/// Returns the Terrestrial Time Julian date for a Modified Julian Date on
|
||
/// the UTC scale, using the offsets in force from 2017.
|
||
///
|
||
/// 使用自 2017 年起生效的偏移，返回 UTC 时标上某约化儒略日对应的地球时儒略日。
|
||
fn tt_julian_date(utc_mjd: f64) -> f64 {
|
||
utc_mjd + 2_400_000.5 + 69.184 / 86_400.0
|
||
}
|
||
|
||
#[test]
|
||
fn the_mean_obliquity_at_j2000_matches_the_reference() {
|
||
// The IAU 2006 obliquity at J2000.0 is 84381.406 arcseconds.
|
||
// J2000.0 处的 IAU 2006 黄赤交角为 84381.406 角秒。
|
||
let obliquity = mean_obliquity(0.0);
|
||
assert!(
|
||
(obliquity / ARCSECOND - 84_381.406).abs() < 1e-9,
|
||
"{} arcseconds",
|
||
obliquity / ARCSECOND
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_nutation_at_j2000_matches_the_reference() {
|
||
// The IAU 2000A nutation at J2000.0, from ERFA.
|
||
// J2000.0 处的 IAU 2000A 章动，取自 ERFA。
|
||
let t = 0.0;
|
||
let nutation = nutation_2000a(t);
|
||
assert!(
|
||
(nutation.longitude - -6.754422426417298e-05).abs() < 1e-16,
|
||
"longitude {}",
|
||
nutation.longitude
|
||
);
|
||
assert!(
|
||
(nutation.obliquity - -2.7970831192374137e-05).abs() < 1e-16,
|
||
"obliquity {}",
|
||
nutation.obliquity
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_nutation_2006_adjustment_is_a_scaling() {
|
||
// The IAU 2006 adjustment scales the IAU 2000A nutation rather than
|
||
// adding an independent term, so the two must stay in proportion.
|
||
// IAU 2006 调整是对 IAU 2000A 章动缩放而非加上独立项，故二者必须保持比例。
|
||
for t in [0.0, 0.5, 1.0, 2.0] {
|
||
let base = nutation_2000a(t);
|
||
let adjusted = nutation_2000a_with_2006_adjustment(t);
|
||
let expected_longitude = base.longitude * (1.0 + 0.469_7e-6 - 2.777_4e-6 * t);
|
||
let expected_obliquity = base.obliquity * (1.0 - 2.777_4e-6 * t);
|
||
assert!(
|
||
(adjusted.longitude - expected_longitude).abs() < 1e-20,
|
||
"t {t}: {} against {expected_longitude}",
|
||
adjusted.longitude
|
||
);
|
||
assert!(
|
||
(adjusted.obliquity - expected_obliquity).abs() < 1e-20,
|
||
"t {t}: {} against {expected_obliquity}",
|
||
adjusted.obliquity
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_nutation_matches_the_reference_over_a_century() {
|
||
// Recorded from ERFA at each date. The values are in radians.
|
||
// 在每个日期从 ERFA 记录。各值以弧度计。
|
||
let cases = [
|
||
(40_587.0, 2.1896610891495316e-05, 4.070477585626008e-05),
|
||
(51_544.5, -6.754419567422628e-05, -2.797091091919891e-05),
|
||
(60_000.0, -4.496357815269203e-05, 3.753534415176348e-05),
|
||
];
|
||
for (mjd, longitude, obliquity) in cases {
|
||
let tt = tt_julian_date(mjd);
|
||
let t = (tt - 2_451_545.0) / 36_525.0;
|
||
let nutation = nutation_2000a(t);
|
||
assert!(
|
||
(nutation.longitude - longitude).abs() < 1e-14,
|
||
"MJD {mjd}: longitude {} against {longitude}",
|
||
nutation.longitude
|
||
);
|
||
assert!(
|
||
(nutation.obliquity - obliquity).abs() < 1e-14,
|
||
"MJD {mjd}: obliquity {} against {obliquity}",
|
||
nutation.obliquity
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_precession_angles_at_j2000_match_the_reference() {
|
||
// The P03 precession angles at J2000.0.
|
||
// J2000.0 处的 P03 岁差角。
|
||
let angles = precession_angles(0.0);
|
||
assert!(
|
||
(angles.gamma_bar - -0.052_928 * ARCSECOND).abs() < 1e-18,
|
||
"gamma bar {}",
|
||
angles.gamma_bar
|
||
);
|
||
assert!(
|
||
(angles.phi_bar - 84_381.412_819 * ARCSECOND).abs() < 1e-15,
|
||
"phi bar {}",
|
||
angles.phi_bar
|
||
);
|
||
assert!(
|
||
(angles.psi_bar - -0.041_775 * ARCSECOND).abs() < 1e-18,
|
||
"psi bar {}",
|
||
angles.psi_bar
|
||
);
|
||
assert!(
|
||
(angles.epsilon - 84_381.406 * ARCSECOND).abs() < 1e-15,
|
||
"epsilon {}",
|
||
angles.epsilon
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_bias_precession_nutation_matrix_matches_the_reference() {
|
||
// The GCRS-to-true-of-date matrix at J2000.0 on the Terrestrial Time
|
||
// scale, from ERFA's `pnm06a`. The values are the first row.
|
||
// 地球时时标上 J2000.0 处的 GCRS 到该日真坐标矩阵，取自 ERFA 的 `pnm06a`。
|
||
// 各值是第一行。
|
||
let tt = 51_544.5 + 2_400_000.5 + 69.184 / 86_400.0;
|
||
let t = (tt - 2_451_545.0) / 36_525.0;
|
||
let matrix = bias_precession_nutation_matrix(t);
|
||
|
||
let expected = [
|
||
[
|
||
0.9999999977211409,
|
||
6.189934767298694e-05,
|
||
2.69478892058386e-05,
|
||
],
|
||
[
|
||
-6.1900102296513e-05,
|
||
0.999999997692101,
|
||
2.8003132867160385e-05,
|
||
],
|
||
[
|
||
-2.6946155767988307e-05,
|
||
-2.8004800880443437e-05,
|
||
0.9999999992448179,
|
||
],
|
||
];
|
||
for (row, expected_row) in expected.iter().enumerate() {
|
||
for (column, expected_value) in expected_row.iter().enumerate() {
|
||
let value = matrix.get(row, column).expect("the indices are in range");
|
||
assert!(
|
||
(value - expected_value).abs() < 1e-14,
|
||
"at ({row}, {column}): {value} against {expected_value}"
|
||
);
|
||
}
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_bias_precession_nutation_matrix_is_a_rotation() {
|
||
for mjd in [40_587.0, 51_544.5, 60_000.0] {
|
||
let tt = tt_julian_date(mjd);
|
||
let t = (tt - 2_451_545.0) / 36_525.0;
|
||
let matrix = bias_precession_nutation_matrix(t);
|
||
assert!(matrix.is_rotation(1e-12), "MJD {mjd}");
|
||
assert!(
|
||
(matrix.determinant() - 1.0).abs() < 1e-12,
|
||
"MJD {mjd}: determinant {}",
|
||
matrix.determinant()
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_celestial_to_terrestrial_matrix_is_a_rotation() {
|
||
for mjd in [40_587.0, 51_544.5, 60_000.0] {
|
||
let tt = tt_julian_date(mjd);
|
||
let matrix = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
|
||
assert!(matrix.is_rotation(1e-9), "MJD {mjd}");
|
||
assert!(
|
||
(matrix.determinant() - 1.0).abs() < 1e-9,
|
||
"MJD {mjd}: determinant {}",
|
||
matrix.determinant()
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_celestial_to_terrestrial_matrix_matches_the_reference() {
|
||
// The matrix at J2000.0 on the Terrestrial Time scale, from ERFA with
|
||
// zero polar motion. Each entry is checked to a nanoradian.
|
||
// 地球时时标上 J2000.0 处的矩阵，取自极移为零的 ERFA。每个元素检查到纳弧度。
|
||
let tt = 51_544.5 + 2_400_000.5 + 69.184 / 86_400.0;
|
||
let matrix = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
|
||
|
||
// The first row, from ERFA's `c2t06a` composed with the equinox based
|
||
// sidereal time.
|
||
// 第一行，取自 ERFA 的 `c2t06a` 与基于春分点的恒星时之合成。
|
||
let expected_first_row = [
|
||
0.18651846184811946,
|
||
-0.982451455739228,
|
||
-2.248740188309595e-05,
|
||
];
|
||
for (column, expected_value) in expected_first_row.iter().enumerate() {
|
||
let value = matrix.get(0, column).expect("the index is in range");
|
||
assert!(
|
||
(value - expected_value).abs() < 1e-9,
|
||
"at (0, {column}): {value} against {expected_value}"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_equation_of_the_equinoxes_includes_the_complementary_terms() {
|
||
// The complementary terms are a few thousandths of an arcsecond and must
|
||
// not be dropped: their omission would show as a millisecond of
|
||
// sidereal time.
|
||
// 补充项是几千分之一角秒，不得省略：略去它们会表现为一毫秒的恒星时。
|
||
let t = 0.0;
|
||
let complementary = equation_of_the_equinoxes_complementary_terms(t);
|
||
assert!(
|
||
complementary.abs() > 1e-9,
|
||
"the complementary terms are {complementary}"
|
||
);
|
||
// At J2000.0 the leading term is 2640.96 microarcseconds times the sine
|
||
// of the Moon's ascending node.
|
||
// 在 J2000.0 处首项是 2640.96 微角秒乘以月球升交点角的正弦。
|
||
let expected = 1.021330096302465e-08;
|
||
assert!(
|
||
(complementary - expected).abs() < 1e-15,
|
||
"{complementary} against {expected}"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_general_precession_is_quadratic_in_time() {
|
||
// The argument is (0.024381750 + 0.00000538691 t) t, so at small times
|
||
// it grows linearly with the leading coefficient.
|
||
// 该辐角为 (0.024381750 + 0.00000538691 t) t，故在小时间处它随首系数线性增长。
|
||
for t in [0.1, 0.5, 1.0] {
|
||
let value = general_precession(t);
|
||
let expected = (0.024_381_750 + 0.000_005_386_91 * t) * t;
|
||
assert!((value - expected).abs() < 1e-18, "t {t}");
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn sidereal_time_advances_by_about_a_turn_per_day() {
|
||
// Greenwich apparent sidereal time must increase by a full turn plus a
|
||
// small amount each day. The small amount is about 3.94 minutes of
|
||
// time, which is 0.9856 degrees.
|
||
// 格林尼治视恒星时每天必须增加一整圈再加上一小量。该小量约为 3.94 分钟时间，
|
||
// 即 0.9856 度。
|
||
let first = greenwich_apparent_sidereal_time(60_000.0, 60_000.0);
|
||
let second = greenwich_apparent_sidereal_time(60_001.0, 60_001.0);
|
||
let advance = (second - first).rem_euclid(2.0 * PI);
|
||
let degrees = advance.to_degrees();
|
||
assert!(
|
||
(degrees - 0.985_6).abs() < 0.01,
|
||
"advance {degrees} degrees"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn zero_polar_motion_is_the_identity_rotation() {
|
||
let matrix = PolarMotion::NONE.matrix();
|
||
for row in 0..3 {
|
||
for column in 0..3 {
|
||
let value = matrix.get(row, column).expect("the indices are in range");
|
||
let expected = if row == column { 1.0 } else { 0.0 };
|
||
assert!((value - expected).abs() < 1e-15, "at ({row}, {column})");
|
||
}
|
||
}
|
||
assert!(PolarMotion::NONE.is_valid());
|
||
}
|
||
|
||
#[test]
|
||
fn a_polar_motion_offsets_the_rotation() {
|
||
// A polar motion of one arcsecond must tilt the matrix by about that
|
||
// much, which is 4.85 microradians.
|
||
// 一角秒的极移必须使矩阵倾斜约那么多，即 4.85 微弧度。
|
||
let motion = PolarMotion::new(ARCSECOND, ARCSECOND, 0.0);
|
||
let matrix = motion.matrix();
|
||
assert!(matrix.is_rotation(1e-12));
|
||
// The off-diagonal term is the sine of the displacement.
|
||
// 非对角项是该位移的正弦。
|
||
let value = matrix.get(0, 2).expect("the index is in range");
|
||
assert!(
|
||
(value - ARCSECOND.sin()).abs() < 1e-15,
|
||
"at (0, 2): {value}"
|
||
);
|
||
assert!(motion.is_valid());
|
||
assert!(!PolarMotion::new(f64::NAN, 0.0, 0.0).is_valid());
|
||
}
|
||
|
||
#[test]
|
||
fn the_polar_motion_shifts_the_earth_fixed_frame() {
|
||
// Applying a polar motion must change the terrestrial position of a
|
||
// fixed inertial direction.
|
||
// 施加极移必须改变某固定惯性方向的地球位置。
|
||
let tt = 60_000.0;
|
||
let without = celestial_to_terrestrial_matrix(tt, tt, PolarMotion::NONE);
|
||
let with = celestial_to_terrestrial_matrix(
|
||
tt,
|
||
tt,
|
||
PolarMotion::new(0.1 * ARCSECOND, -0.2 * ARCSECOND, 0.0),
|
||
);
|
||
let difference = (without.transpose() * with).transpose() * with;
|
||
// The two must differ, and the difference must be of the order of the
|
||
// applied angles.
|
||
// 二者必须不同，且差值必须与所施加的角度同量级。
|
||
let mut largest = 0.0_f64;
|
||
for row in 0..3 {
|
||
for column in 0..3 {
|
||
let a = without.get(row, column).expect("in range");
|
||
let b = difference.get(row, column).expect("in range");
|
||
largest = largest.max((a - b).abs());
|
||
}
|
||
}
|
||
assert!(largest > 1e-8, "the largest difference is {largest}");
|
||
assert!(largest < 1e-5, "the largest difference is {largest}");
|
||
}
|
||
|
||
#[test]
|
||
fn the_nutation_reaches_the_expected_magnitude() {
|
||
// The nutation in longitude reaches about seventeen arcseconds and the
|
||
// nutation in obliquity about nine, so over a range of dates the
|
||
// largest value found must be in that range.
|
||
// 黄经章动达约十七角秒，黄赤交角章动约九角秒，故在一段日期上所能找到的最大值
|
||
// 必须落在该范围内。
|
||
let mut largest_longitude = 0.0_f64;
|
||
let mut largest_obliquity = 0.0_f64;
|
||
for step in 0..400 {
|
||
// A span of about fifty years in steps of a fortnight.
|
||
// 约五十年、步长为两周的跨度。
|
||
let t = f64::from(step) * 0.125;
|
||
let nutation = nutation_2000a(t);
|
||
largest_longitude = largest_longitude.max(nutation.longitude.abs());
|
||
largest_obliquity = largest_obliquity.max(nutation.obliquity.abs());
|
||
}
|
||
let longitude_arcseconds = largest_longitude / ARCSECOND;
|
||
let obliquity_arcseconds = largest_obliquity / ARCSECOND;
|
||
assert!(
|
||
(19.0..20.0).contains(&longitude_arcseconds),
|
||
"longitude {longitude_arcseconds} arcseconds"
|
||
);
|
||
assert!(
|
||
(8.0..10.0).contains(&obliquity_arcseconds),
|
||
"obliquity {obliquity_arcseconds} arcseconds"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_nutation_is_periodic() {
|
||
// The dominant terms have periods of about eighteen years and half a
|
||
// year, so the value must return close to its initial value after
|
||
// eighteen years and must oscillate in between.
|
||
// 主导项的周期约为十八年与半年，故该值必须在十八年后回到接近其初值，并在其间
|
||
// 振荡。
|
||
let initial = nutation_2000a(0.0).longitude;
|
||
let after_eighteen_years = nutation_2000a(0.18).longitude;
|
||
assert!(
|
||
(after_eighteen_years - initial).abs() < 1.5e-05,
|
||
"{after_eighteen_years} against {initial}"
|
||
);
|
||
|
||
// The sign must change at some point within that span.
|
||
// 符号必须在该跨度内某处改变。
|
||
let mut seen_positive = false;
|
||
let mut seen_negative = false;
|
||
for step in 0..200 {
|
||
let t = f64::from(step) * 0.001;
|
||
let value = nutation_2000a(t).longitude;
|
||
seen_positive |= value > 0.0;
|
||
seen_negative |= value < 0.0;
|
||
}
|
||
assert!(seen_positive && seen_negative);
|
||
}
|
||
|
||
#[test]
|
||
fn the_equation_of_the_equinoxes_follows_the_nutation() {
|
||
// The equation of the equinoxes is dominated by the nutation in
|
||
// longitude projected onto the true equator, so it shares its sign.
|
||
// 分点方程式主要由投影到真赤道上的黄经章动主导，故它与后者同号。
|
||
for t in [0.02, 0.05, 0.1, 0.2] {
|
||
let nutation = nutation_2000a_with_2006_adjustment(t);
|
||
let equation = equation_of_the_equinoxes(t);
|
||
assert!(
|
||
nutation.longitude.signum() == equation.signum(),
|
||
"t {t}: {} against {}",
|
||
nutation.longitude,
|
||
equation
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn a_matrix_at_two_epochs_separated_by_a_century_differs_by_the_precession() {
|
||
// The precession over a century is about 1.4 degrees, so the two
|
||
// matrices must differ by an angle of that order.
|
||
// 一世纪的岁差约为 1.4 度，故两个矩阵必须相差该量级的角。
|
||
let first = bias_precession_nutation_matrix(0.0);
|
||
let second = bias_precession_nutation_matrix(1.0);
|
||
let difference = first.transpose() * second;
|
||
// The trace of a rotation matrix is 1 + 2 cos(angle), so the angle
|
||
// follows from it.
|
||
// 转动矩阵的迹为 1 + 2 cos(angle)，故该角由它得出。
|
||
let cosine = ((difference.trace() - 1.0) / 2.0).clamp(-1.0, 1.0);
|
||
let degrees = cosine.acos().to_degrees();
|
||
assert!(
|
||
(1.0..2.0).contains(°rees),
|
||
"the rotation angle is {degrees} degrees"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_1982_sidereal_time_matches_the_reference() {
|
||
// Greenwich mean sidereal time from the 1982 expression, which is the
|
||
// one the SGP4 theory was fitted against. The values are from ERFA.
|
||
// 1982 表达式给出的格林尼治平恒星时，即 SGP4 理论所拟合者。各值取自 ERFA。
|
||
let cases = [
|
||
(51_544.5, 4.894_961_212_823_059),
|
||
(60_476.999_999_810_5, 4.620_858_198_832_686),
|
||
(40_587.0, 1.749_337_177_344_97),
|
||
];
|
||
for (mjd, expected) in cases {
|
||
let actual = apsis_frames::iau2000::greenwich_mean_sidereal_time_1982(2_400_000.5, mjd);
|
||
assert!(
|
||
(actual - expected).abs() < 1e-15,
|
||
"MJD {mjd}: {actual} against {expected}"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_1982_sidereal_time_is_a_full_turn_per_day() {
|
||
// Sidereal time advances by one turn plus the excess of the sidereal
|
||
// day over the solar day, which is 0.9856 degrees per day.
|
||
// 恒星时推进一整圈再加上恒星日超出太阳日的部分，即每天 0.9856 度。
|
||
let first = apsis_frames::iau2000::greenwich_mean_sidereal_time_1982(2_400_000.5, 60_000.0);
|
||
let second = apsis_frames::iau2000::greenwich_mean_sidereal_time_1982(2_400_000.5, 60_001.0);
|
||
let advance = (second - first).rem_euclid(2.0 * PI);
|
||
assert!(
|
||
(advance.to_degrees() - 0.985_6).abs() < 0.001,
|
||
"advance {} degrees",
|
||
advance.to_degrees()
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_teme_transformation_to_earth_fixed_matches_the_reference() {
|
||
use apsis_frames::iau2000::{PolarMotion, teme_to_terrestrial_matrix};
|
||
use apsis_math::Vec3;
|
||
|
||
// The reference values are ERFA's `gmst82` applied as a third axis
|
||
// rotation with the identity polar motion.
|
||
// 参考值是把 ERFA 的 `gmst82` 作为第三轴转动、配以单位极移施加所得。
|
||
let cases = [
|
||
(
|
||
60_481.5,
|
||
[7_000.0, 1_000.0, 2_000.0],
|
||
[1_098.716_242_792, -6_985.185_940_104, 2_000.0],
|
||
),
|
||
(
|
||
51_544.5,
|
||
[7_000.0, -3_000.0, 1_000.0],
|
||
[4_221.057_369_867, 6_338.980_571_061, 1_000.0],
|
||
),
|
||
];
|
||
for (mjd, position, expected) in cases {
|
||
let matrix = teme_to_terrestrial_matrix(2_400_000.5, mjd, PolarMotion::NONE);
|
||
let result = matrix * Vec3::new(position[0], position[1], position[2]);
|
||
let expected = Vec3::new(expected[0], expected[1], expected[2]);
|
||
assert!(
|
||
(result - expected).norm() < 1e-6,
|
||
"MJD {mjd}: {result} against {expected}"
|
||
);
|
||
}
|
||
}
|
||
|
||
#[test]
|
||
fn the_teme_transformation_preserves_length() {
|
||
use apsis_frames::iau2000::{PolarMotion, teme_to_terrestrial_matrix};
|
||
use apsis_math::Vec3;
|
||
|
||
let matrix = teme_to_terrestrial_matrix(2_400_000.5, 60_000.0, PolarMotion::NONE);
|
||
assert!(matrix.is_rotation(1e-12));
|
||
let vector = Vec3::new(7_000.0, -1_000.0, 500.0);
|
||
assert!(
|
||
((matrix * vector).norm() - vector.norm()).abs() < 1e-9,
|
||
"the length changed"
|
||
);
|
||
}
|
||
|
||
#[test]
|
||
fn the_teme_frame_uses_the_mean_equinox_not_the_true_one() {
|
||
// The distinction matters. TEME's first axis is the mean equinox, so it
|
||
// pairs with a mean sidereal time. Substituting the apparent value, or
|
||
// the modern mean value, rotates the result by a fraction of a degree,
|
||
// which is many kilometres of ground position.
|
||
// 该区别要紧。TEME 的第一轴是平春分点，故它与平恒星时配对。换成视值或现代平值
|
||
// 会把结果转动零点几度，即许多千米的地面位置。
|
||
use apsis_frames::iau2000::{
|
||
greenwich_apparent_sidereal_time, greenwich_mean_sidereal_time,
|
||
greenwich_mean_sidereal_time_1982,
|
||
};
|
||
|
||
let mjd = 60_481.5;
|
||
let tt = mjd + 2_400_000.5 + 69.184 / 86_400.0;
|
||
let sidereal_1982 = greenwich_mean_sidereal_time_1982(2_400_000.5, mjd);
|
||
let apparent = greenwich_apparent_sidereal_time(tt, tt);
|
||
let modern_mean = greenwich_mean_sidereal_time(tt, tt);
|
||
|
||
// The apparent value includes the equation of the equinoxes, which is a
|
||
// few thousandths of a degree.
|
||
// 视值包含分点方程式，即千分之几度。
|
||
let equation = (apparent - modern_mean).abs().to_degrees();
|
||
assert!(
|
||
equation > 1e-4 && equation < 1e-2,
|
||
"the equation of the equinoxes is {equation} degrees"
|
||
);
|
||
|
||
// The 1982 expression, which TEME is fitted against, differs from the
|
||
// modern mean value by a fraction of a degree, because the precession
|
||
// model changed between them.
|
||
// TEME 所拟合的 1982 表达式与现代平值相差零点几度，因为二者之间的岁差模型
|
||
// 改变了。
|
||
let secular = (sidereal_1982 - modern_mean).abs().to_degrees();
|
||
assert!(
|
||
secular > 0.1 && secular < 1.0,
|
||
"the difference is {secular} degrees"
|
||
);
|
||
}
|