1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
175
176
177
178
179
use crate::astronomy::eval_polynomial;
/// Returns the difference (ΔT = TD - UT) between dynamic time (TD) and
/// universal time (UT) in seconds at a specific Julian year.
pub fn delta_t(y: f64) -> f64 {
// polynomials are from http://eclipse.gsfc.nasa.gov/SEcat5/deltatpoly.html
// and http://www.staff.science.uu.nl/~gent0113/deltat/deltat_old.htm
// using the ΔT values derived from historical record and from direct observations,
// a series of polynomial expressions have been created to simplify the evaluation of ΔT
// this gives `y` for the middle of the month, which is accurate enough given the precision
// in the known values of ΔT
// values derived from [Table 1] and [Table 2] from https://eclipse.gsfc.nasa.gov/SEcat5/deltat.html
if y < -500.0 {
// ΔT = -20 + 32 * u^2
// where: u = (y - 1820)/100
eval_polynomial((y - 1820.0) * 0.01, &[-20.0, 0.0, 32.0])
} else if y < 500.0 {
// Between years -500 and +500, we use the data from [Table 1], except that for the
// year -500 we changed value 17190 to 17203.7 to avoid a discontinuity with the
// previous formula at that epoch. The value for ΔT is given by a polynomial of
// the 6th degree, which reproduces the values in [Table 1] with an error not
// larger than 4 seconds:
// ΔT = 10583.6 - 1014.41 * u + 33.78311 * u^2 - 5.952053 * u^3 - 0.1798452 * u^4 + 0.022174192 * u^5 + 0.0090316521 * u^6
// where: u = y/100
eval_polynomial(y * 0.01, &[
10583.6,
-1014.41,
33.78311,
-5.952053,
-0.1798452,
0.022174192,
0.0090316521,
])
} else if y < 1600.0 {
// ΔT = 1574.2 - 556.01 * u + 71.23472 * u^2 + 0.319781 * u^3 - 0.8503463 * u^4 - 0.005050998 * u^5 + 0.0083572073 * u^6
// where: u = (y - 1000)/100
eval_polynomial((y - 1000.0) * 0.01, &[
1574.2,
-556.01,
71.23472,
0.319781,
-0.8503463,
-0.005050998,
0.0083572073,
])
} else if y < 1700.0 {
// ΔT = 120 - 0.9808 * u - 0.01532 * u^2 + u^3 / 7129
// where: u = y - 1600
eval_polynomial(y - 1600.0, &[120.0, -0.9808, -0.01532, 1.0 / 7129.0])
} else if y < 1800.0 {
// ΔT = 8.83 + 0.1603 * u - 0.0059285 * u^2 + 0.00013336 * u^3 - u^4 / 1174000
// where: u = y - 1700
eval_polynomial(y - 1700.0, &[
8.83,
0.1603,
-0.0059285,
0.00013336,
-1.0 / 1174000.0,
])
} else if y < 1860.0 {
// ΔT = 13.72 - 0.332447 * u + 0.0068612 * u^2 + 0.0041116 * u^3 - 0.00037436 * u^4 + 0.0000121272 * u^5 - 0.0000001699 * u^6 + 0.000000000875 * u^7
// where: u = y - 1800
eval_polynomial(y - 1800.0, &[
13.72,
-0.332447,
0.0068612,
0.0041116,
-0.00037436,
0.0000121272,
-0.0000001699,
0.000000000875,
])
} else if y < 1900.0 {
// ΔT = 7.62 + 0.5737 * u - 0.251754 * u^2 + 0.01680668 * u^3 - 0.0004473624 * u^4 + u^5 / 233174
// where: u = y - 1860
eval_polynomial(y - 1860.0, &[
7.62,
0.5737,
-0.251754,
0.01680668,
-0.0004473624,
1.0 / 233174.0,
])
} else if y < 1920.0 {
// ΔT = -2.79 + 1.494119 * u - 0.0598939 * u^2 + 0.0061966 * u^3 - 0.000197 * u^4
// where: u = y - 1900
eval_polynomial(y - 1900.0, &[
-2.79, 1.494119, -0.0598939, 0.0061966, -0.000197,
])
} else if y < 1941.0 {
// ΔT = 21.20 + 0.84493*u - 0.076100 * u^2 + 0.0020936 * u^3
// where: u = y - 1920
eval_polynomial(y - 1920.0, &[21.20, 0.84493, 0.076100, 0.0020936])
} else if y < 1961.0 {
// ΔT = 29.07 + 0.407*u - u^2/233 + u^3 / 2547
// where: u = y - 1950
eval_polynomial(y - 1950.0, &[29.07, 0.407, 1.0 / 233.0, 1.0 / 2547.0])
} else if y < 1986.0 {
// ΔT = 45.45 + 1.067*u - u^2/260 - u^3 / 718
// where: u = y - 1975
eval_polynomial(y - 1975.0, &[45.45, 1.067, 1.0 / 260.0, 1.0 / 718.0])
} else if y < 2005.0 {
// ΔT = 63.86 + 0.3345 * u - 0.060374 * u^2 + 0.0017275 * u^3 + 0.000651814 * u^4 + 0.00002373599 * u^5
// where: u = y - 2000
eval_polynomial(y - 2000.0, &[
63.86,
0.3345,
-0.060374,
0.0017275,
0.000651814,
0.00002373599,
])
} else if y < 2050.0 {
// ΔT = 62.92 + 0.32217 * u + 0.005589 * u^2
// where: u = y - 2000
eval_polynomial(y - 2000.0, &[62.92, 0.32217, 0.005589])
} else if y < 2150.0 {
// ΔT = -20 + 32 * ((y-1820)/100)^2 - 0.5628 * (2150 - y)
0.5628f64.mul_add(
-(2150.0 - y),
32.0f64.mul_add(((y - 1820.0) * 0.01).powi(2), -20.0),
)
} else {
// ΔT = -20 + 32 * u^2
// where: u = (y-1820)/100
eval_polynomial((y - 1820.0) * 0.01, &[-20.0, 0.0, 32.0])
}
}
#[cfg(test)]
mod tests {
use float_cmp::assert_approx_eq;
use test_case::test_case;
use super::delta_t;
#[test_case(-501, 11, 17205.536050000002)]
#[test_case(-480, 3, 16844.22790443188)]
#[test_case(-240, 5, 13284.287479640398)]
#[test_case(0, 7, 10578.106269426406)]
#[test_case(371, 3, 6984.200618855546)]
#[test_case(501, 6, 5695.587627439579)]
#[test_case(985, 12, 1653.6763788571984)]
#[test_case(1000, 12, 1568.8781133240834)]
#[test_case(1120, 2, 1007.880950144438)]
#[test_case(1400, 1, 321.6919250528026)]
#[test_case(1750, 1, 13.375979008768494)]
#[test_case(1820, 3, 11.806505586706487)]
#[test_case(1905, 4, 4.2029658668397705)]
#[test_case(1921, 7, 22.691141267419063)]
#[test_case(1943, 8, 26.54807704872117)]
#[test_case(1987, 2, 55.37182043227025)]
#[test_case(2007, 8, 65.70149420312501)]
#[test_case(2052, 1, 97.16772222222194)]
#[test_case(2162, 3, 354.7409388888892)]
fn expect_delta_t_to_match(year: i32, month: u8, expected: f64) {
let y = f64::from(year) + ((f64::from(month) - 0.5) / 12.0);
assert_approx_eq!(f64, delta_t(y), expected, ulps = 1);
}
}