sofars/ts/dat.rs
1/// For a given UTC date, calculate Delta(AT) = TAI-UTC.
2/// ```
3/// :------------------------------------------:
4/// : :
5/// : IMPORTANT :
6/// : :
7/// : A new version of this function must be :
8/// : produced whenever a new leap second is :
9/// : announced. There are four items to :
10/// : change on each such occasion: :
11/// : :
12/// : 1) A new line must be added to the set :
13/// : of statements that initialize the :
14/// : array "changes". :
15/// : :
16/// : 2) The constant IYV must be set to the :
17/// : current year. :
18/// : :
19/// : 3) The "Latest leap second" comment :
20/// : below must be set to the new leap :
21/// : second date. :
22/// : :
23/// : 4) The "This revision" comment, later, :
24/// : must be set to the current date. :
25/// : :
26/// : Change (2) must also be carried out :
27/// : whenever the function is re-issued, :
28/// : even if no leap seconds have been :
29/// : added. :
30/// : :
31/// : Latest leap second: 2016 December 31 :
32/// : :
33/// :__________________________________________:
34/// ```
35/// This function is part of the International Astronomical Union's
36/// SOFA (Standards of Fundamental Astronomy) software collection.
37///
38/// Status: user-replaceable support function.
39///
40/// Given:
41/// ```
42/// iy int UTC: year (Notes 1 and 2)
43/// im int month (Note 2)
44/// id int day (Notes 2 and 3)
45/// fd double fraction of day (Note 4)
46/// ```
47/// Returned:
48/// ```
49/// deltat double TAI minus UTC, seconds
50/// ```
51/// Returned (function value):
52/// ```
53/// int status (Note 5):
54/// 1 = dubious year (Note 1)
55/// 0 = OK
56/// -1 = bad year
57/// -2 = bad month
58/// -3 = bad day (Note 3)
59/// -4 = bad fraction (Note 4)
60/// -5 = internal error (Note 5)
61/// ```
62/// Notes:
63///
64/// 1) UTC began at 1960 January 1.0 (JD 2436934.5) and it is improper
65/// to call the function with an earlier date. If this is attempted,
66/// zero is returned together with a warning status.
67///
68/// Because leap seconds cannot, in principle, be predicted in
69/// advance, a reliable check for dates beyond the valid range is
70/// impossible. To guard against gross errors, a year five or more
71/// after the release year of the present function (see the constant
72/// IYV) is considered dubious. In this case a warning status is
73/// returned but the result is computed in the normal way.
74///
75/// For both too-early and too-late years, the warning status is +1.
76/// This is distinct from the error status -1, which signifies a year
77/// so early that JD could not be computed.
78///
79/// 2) If the specified date is for a day which ends with a leap second,
80/// the TAI-UTC value returned is for the period leading up to the
81/// leap second. If the date is for a day which begins as a leap
82/// second ends, the TAI-UTC returned is for the period following the
83/// leap second.
84///
85/// 3) The day number must be in the normal calendar range, for example
86/// 1 through 30 for April. The "almanac" convention of allowing
87/// such dates as January 0 and December 32 is not supported in this
88/// function, in order to avoid confusion near leap seconds.
89///
90/// 4) The fraction of day is used only for dates before the
91/// introduction of leap seconds, the first of which occurred at the
92/// end of 1971. It is tested for validity (0 to 1 is the valid
93/// range) even if not used; if invalid, zero is used and status -4
94/// is returned. For many applications, setting fd to zero is
95/// acceptable; the resulting error is always less than 3 ms (and
96/// occurs only pre-1972).
97///
98/// 5) The status value returned in the case where there are multiple
99/// errors refers to the first error detected. For example, if the
100/// month and day are 13 and 32 respectively, status -2 (bad month)
101/// will be returned. The "internal error" status refers to a
102/// case that is impossible but causes some compilers to issue a
103/// warning.
104///
105/// 6) In cases where a valid result is not available, zero is returned.
106///
107/// References:
108///
109/// 1) For dates from 1961 January 1 onwards, the expressions from the
110/// file ftp://maia.usno.navy.mil/ser7/tai-utc.dat are used.
111///
112/// 2) The 5ms timestep at 1961 January 1 is taken from 2.58.1 (p87) of
113/// the 1992 Explanatory Supplement.
114///
115/// Called:
116/// iauCal2jd Gregorian calendar to JD
117pub fn dat(iy: i32, im: i32, id: i32, fd: f64) -> Result<f64, i32> {
118 // Release year for this version of iauDat
119 const IYV: i32 = 2023;
120
121 // Reference dates (MJD) and drift rates (s/day), pre leap seconds
122 const DRIFT: &[(f64, f64)] = &[
123 (37300.0, 0.0012960),
124 (37300.0, 0.0012960),
125 (37300.0, 0.0012960),
126 (37665.0, 0.0011232),
127 (37665.0, 0.0011232),
128 (38761.0, 0.0012960),
129 (38761.0, 0.0012960),
130 (38761.0, 0.0012960),
131 (38761.0, 0.0012960),
132 (38761.0, 0.0012960),
133 (38761.0, 0.0012960),
134 (38761.0, 0.0012960),
135 (39126.0, 0.0025920),
136 (39126.0, 0.0025920),
137 ];
138
139 // Number of Delta(AT) expressions before leap seconds were introduced
140 const NERA1: usize = DRIFT.len();
141
142 // Dates and Delta(AT)s
143 const CHANGES: &[(i32, i32, f64)] = &[
144 (1960, 1, 1.4178180),
145 (1961, 1, 1.4228180),
146 (1961, 8, 1.3728180),
147 (1962, 1, 1.8458580),
148 (1963, 11, 1.9458580),
149 (1964, 1, 3.2401300),
150 (1964, 4, 3.3401300),
151 (1964, 9, 3.4401300),
152 (1965, 1, 3.5401300),
153 (1965, 3, 3.6401300),
154 (1965, 7, 3.7401300),
155 (1965, 9, 3.8401300),
156 (1966, 1, 4.3131700),
157 (1968, 2, 4.2131700),
158 (1972, 1, 10.0),
159 (1972, 7, 11.0),
160 (1973, 1, 12.0),
161 (1974, 1, 13.0),
162 (1975, 1, 14.0),
163 (1976, 1, 15.0),
164 (1977, 1, 16.0),
165 (1978, 1, 17.0),
166 (1979, 1, 18.0),
167 (1980, 1, 19.0),
168 (1981, 7, 20.0),
169 (1982, 7, 21.0),
170 (1983, 7, 22.0),
171 (1985, 7, 23.0),
172 (1988, 1, 24.0),
173 (1990, 1, 25.0),
174 (1991, 1, 26.0),
175 (1992, 7, 27.0),
176 (1993, 7, 28.0),
177 (1994, 7, 29.0),
178 (1996, 1, 30.0),
179 (1997, 7, 31.0),
180 (1999, 1, 32.0),
181 (2006, 1, 33.0),
182 (2009, 1, 34.0),
183 (2012, 7, 35.0),
184 (2015, 7, 36.0),
185 (2017, 1, 37.0),
186 ];
187
188 // Number of Delta(AT) changes
189 const NDAT: usize = CHANGES.len();
190
191 // Miscellaneous local variables
192 let mut i: usize;
193 let m: i32;
194 let mut da: f64;
195
196 // Initialize the result to zero.
197 let mut deltat: f64 = 0.0;
198
199 // If invalid fraction of a day, set error status and give up.
200 if fd < 0.0 || fd > 1.0 {
201 return Err(-4);
202 }
203
204 use crate::cal::cal2jd;
205 // Convert the date into an MJD.
206 let j = cal2jd(iy, im, id);
207
208 // If invalid year, month, or day, give up.
209 if j.is_err() {
210 return Err(j.err().unwrap());
211 }
212
213 let (_, djm) = j.unwrap();
214
215 // If pre-UTC year, set warning status and give up.
216 if iy < CHANGES[0].0 {
217 return Err(1);
218 }
219
220 let j: i32;
221 // If suspiciously late year, set warning status but proceed.
222 if iy > IYV + 5 {
223 j = 1;
224 }
225
226 // Combine year and month to form a date-ordered integer...
227 m = 12 * iy + im;
228
229 // ...and use it to find the preceding table entry.
230 i = NDAT - 1;
231 while i > 0 && m < (12 * CHANGES[i].0 + CHANGES[i].1) {
232 i -= 1;
233 }
234
235 // Prevent underflow warnings.
236 if i == 0 && m < (12 * CHANGES[i].0 + CHANGES[i].1) {
237 return Err(-5);
238 }
239
240 // Get the Delta(AT).
241 da = CHANGES[i].2;
242
243 // If pre-1972, adjust for drift.
244 if i < NERA1 {
245 da += (djm + fd - DRIFT[i].0) * DRIFT[i].1;
246 }
247
248 // Return the Delta(AT) value.
249 deltat = da;
250
251 // Return the status.
252 Ok(deltat)
253}