1 /*-
2 * SPDX-License-Identifier: BSD-3-Clause
3 *
4 * Copyright (c) 2019-2020 The DragonFly Project. All rights reserved.
5 *
6 * This code is derived from software contributed to The DragonFly Project
7 * by Aaron LI <aly@aaronly.me>
8 *
9 * Redistribution and use in source and binary forms, with or without
10 * modification, are permitted provided that the following conditions
11 * are met:
12 *
13 * 1. Redistributions of source code must retain the above copyright
14 * notice, this list of conditions and the following disclaimer.
15 * 2. Redistributions in binary form must reproduce the above copyright
16 * notice, this list of conditions and the following disclaimer in
17 * the documentation and/or other materials provided with the
18 * distribution.
19 * 3. Neither the name of The DragonFly Project nor the names of its
20 * contributors may be used to endorse or promote products derived
21 * from this software without specific, prior written permission.
22 *
23 * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS
24 * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT
25 * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS
26 * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE
27 * COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT,
28 * INCIDENTAL, SPECIAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING,
29 * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES;
30 * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED
31 * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY,
32 * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT
33 * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF
34 * SUCH DAMAGE.
35 *
36 * Reference:
37 * Calendrical Calculations, The Ultimate Edition (4th Edition)
38 * by Edward M. Reingold and Nachum Dershowitz
39 * 2018, Cambridge University Press
40 */
41
42 #include <err.h>
43 #include <math.h>
44 #include <stdio.h>
45
46 #include "basics.h"
47 #include "gregorian.h"
48 #include "utils.h"
49
50 /*
51 * Determine the day of week of the fixed date $rd.
52 * Ref: Sec.(1.12), Eq.(1.60)
53 */
54 int
dayofweek_from_fixed(int rd)55 dayofweek_from_fixed(int rd)
56 {
57 /* NOTE: R.D. 1 is Monday */
58 return mod(rd, 7);
59 }
60
61 /*
62 * Calculate the fixed date of the day-of-week $dow on or before
63 * the fixed date $rd.
64 * Ref: Sec.(1.12), Eq.(1.62)
65 */
66 int
kday_onbefore(int dow,int rd)67 kday_onbefore(int dow, int rd)
68 {
69 return rd - dayofweek_from_fixed(rd - dow);
70 }
71
72 /*
73 * Calculate the fixed date of the day-of-week $dow after
74 * the fixed date $rd.
75 * Ref: Sec.(1.12), Eq.(1.68)
76 */
77 int
kday_after(int dow,int rd)78 kday_after(int dow, int rd)
79 {
80 return kday_onbefore(dow, rd + 7);
81 }
82
83 /*
84 * Calculate the fixed date of the day-of-week $dow nearest to
85 * the fixed date $rd.
86 * Ref: Sec.(1.12), Eq.(1.66)
87 */
88 int
kday_nearest(int dow,int rd)89 kday_nearest(int dow, int rd)
90 {
91 return kday_onbefore(dow, rd + 3);
92 }
93
94 /*
95 * Calculate the $n-th occurrence of a given day-of-week $dow counting
96 * from either after (if $n > 0) or before (if $n < 0) the given $date.
97 * Ref: Sec.(2.5), Eq.(2.33)
98 */
99 int
nth_kday(int n,int dow,struct date * date)100 nth_kday(int n, int dow, struct date *date)
101 {
102 if (n == 0)
103 errx(1, "%s: invalid n = 0!", __func__);
104
105 int rd = fixed_from_gregorian(date);
106 int kday;
107 if (n > 0)
108 kday = kday_onbefore(dow, rd - 1); /* Sec.(1.12), Eq.(1.67) */
109 else
110 kday = kday_onbefore(dow, rd + 7); /* Sec.(1.12), Eq.(1.68) */
111 return 7 * n + kday;
112 }
113
114 /*
115 * Calculate the ephemeris correction (fraction of day) required for
116 * converting between Universal Time and Dynamical Time at moment $t.
117 * Ref: Sec.(14.2), Eq.(14.15)
118 */
119 double
ephemeris_correction(double t)120 ephemeris_correction(double t)
121 {
122 int rd = (int)floor(t);
123 int year = gregorian_year_from_fixed(rd);
124 int y2000 = year - 2000;
125 int y1700 = year - 1700;
126 int y1600 = year - 1600;
127 double y1820 = (year - 1820) / 100.0;
128 double y1000 = (year - 1000) / 100.0;
129 double y0 = year / 100.0;
130
131 double coef2006[] = { 62.92, 0.32217, 0.005589 };
132 double coef1987[] = { 63.86, 0.3345, -0.060374, 0.0017275,
133 0.000651814, 0.00002373599 };
134 double coef1900[] = { -0.00002, 0.000297, 0.025184, -0.181133,
135 0.553040, -0.861938, 0.677066, -0.212591 };
136 double coef1800[] = { -0.000009, 0.003844, 0.083563, 0.865736,
137 4.867575, 15.845535, 31.332267, 38.291999,
138 28.316289, 11.636204, 2.043794 };
139 double coef1700[] = { 8.118780842, -0.005092142,
140 0.003336121, -0.0000266484 };
141 double coef1600[] = { 120.0, -0.9808, -0.01532, 0.000140272128 };
142 double coef500[] = { 1574.2, -556.01, 71.23472, 0.319781,
143 -0.8503463, -0.005050998, 0.0083572073 };
144 double coef0[] = { 10583.6, -1014.41, 33.78311, -5.952053,
145 -0.1798452, 0.022174192, 0.0090316521 };
146
147 double c_other = (-20.0 + 32.0 * y1820 * y1820) / 86400.0;
148 struct date date1 = { 1900, 1, 1 };
149 struct date date2 = { year, 7, 1 };
150 double c = gregorian_date_difference(&date1, &date2) / 36525.0;
151
152 if (year > 2150) {
153 return c_other;
154 } else if (year >= 2051) {
155 return c_other + 0.5628 * (2150 - year) / 86400.0;
156 } else if (year >= 2006) {
157 return poly(y2000, coef2006, nitems(coef2006)) / 86400.0;
158 } else if (year >= 1987) {
159 return poly(y2000, coef1987, nitems(coef1987)) / 86400.0;
160 } else if (year >= 1900) {
161 return poly(c, coef1900, nitems(coef1900));
162 } else if (year >= 1800) {
163 return poly(c, coef1800, nitems(coef1800));
164 } else if (year >= 1700) {
165 return poly(y1700, coef1700, nitems(coef1700)) / 86400.0;
166 } else if (year >= 1600) {
167 return poly(y1600, coef1600, nitems(coef1600)) / 86400.0;
168 } else if (year >= 500) {
169 return poly(y1000, coef500, nitems(coef500)) / 86400.0;
170 } else if (year > -500) {
171 return poly(y0, coef0, nitems(coef0)) / 86400.0;
172 } else {
173 return c_other;
174 }
175 }
176
177 /*
178 * Convert from Universal Time (UT) to Dynamical Time (DT).
179 * Ref: Sec.(14.2), Eq.(14.16)
180 */
181 double
dynamical_from_universal(double t)182 dynamical_from_universal(double t)
183 {
184 return t + ephemeris_correction(t);
185 }
186
187 /*
188 * Convert from dynamical time to universal time.
189 * Ref: Sec.(14.2), Eq.(14.17)
190 */
191 double
universal_from_dynamical(double t)192 universal_from_dynamical(double t)
193 {
194 return t - ephemeris_correction(t);
195 }
196
197 /*
198 * Calculate the number (and fraction) of uniform-length centuries
199 * (36525 days) since noon on January 1, 2000 (Gregorian) at moment $t.
200 * Ref: Sec.(14.2), Eq.(14.18)
201 */
202 double
julian_centuries(double t)203 julian_centuries(double t)
204 {
205 double dt = dynamical_from_universal(t);
206 double j2000 = 0.5 + gregorian_new_year(2000);
207 return (dt - j2000) / 36525.0;
208 }
209
210 /*
211 * Calculate the mean sidereal time of day expressed as hour angle
212 * at moment $t.
213 * Ref: Sec.(14.3), Eq.(14.27)
214 */
215 double
sidereal_from_moment(double t)216 sidereal_from_moment(double t)
217 {
218 double j2000 = 0.5 + gregorian_new_year(2000);
219 int century_days = 36525;
220 double c = (t - j2000) / century_days;
221 double coef[] = { 280.46061837, 360.98564736629 * century_days,
222 0.000387933, -1.0 / 38710000.0 };
223 return mod_f(poly(c, coef, nitems(coef)), 360);
224 }
225
226 /*
227 * Ref: Sec.(14.3), Eq.(14.20)
228 */
229 double
equation_of_time(double t)230 equation_of_time(double t)
231 {
232 double c = julian_centuries(t);
233 double epsilon = obliquity(t);
234 double y = pow(tan_deg(epsilon/2), 2);
235
236 double lambda = 280.46645 + 36000.76983 * c + 0.0003032 * c*c;
237 double anomaly = (357.52910 + 35999.05030 * c -
238 0.0001559 * c*c - 0.00000048 * c*c*c);
239 double eccentricity = (0.016708617 - 0.000042037 * c -
240 0.0000001236 * c*c);
241
242 double equation = (y * sin_deg(2*lambda) -
243 2 * eccentricity * sin_deg(anomaly) +
244 4 * eccentricity * y * sin_deg(anomaly) * cos_deg(2*lambda) -
245 1.25 * eccentricity*eccentricity * sin_deg(2*anomaly) -
246 0.5 * y*y * sin_deg(4*lambda)) / (2*M_PI);
247
248 double vmax = 0.5; /* i.e., 12 hours */
249 if (fabs(equation) < vmax)
250 return equation;
251 else
252 return (equation > 0) ? vmax : -vmax;
253 }
254
255 /*
256 * Calculate the sundial time from local time $t at location of
257 * longitude $longitude.
258 * Ref: Sec.(14.3), Eq.(14.21)
259 */
260 double
apparent_from_local(double t,double longitude)261 apparent_from_local(double t, double longitude)
262 {
263 double ut = t - longitude / 360.0; /* local time -> universal time */
264 return t + equation_of_time(ut);
265 }
266
267 /*
268 * Calculate the local time from sundial time $t at location of
269 * longitude $longitude.
270 * Ref: Sec.(14.3), Eq.(14.22)
271 */
272 double
local_from_apparent(double t,double longitude)273 local_from_apparent(double t, double longitude)
274 {
275 double ut = t - longitude / 360.0; /* local time -> universal time */
276 return t - equation_of_time(ut);
277 }
278
279 /*
280 * Calculate the obliquity of ecliptic at moment $t
281 * Ref: Sec.(14.4), Eq.(14.28)
282 */
283 double
obliquity(double t)284 obliquity(double t)
285 {
286 double c = julian_centuries(t);
287 double coef[] = {
288 0.0,
289 angle2deg(0, 0, -46.8150),
290 angle2deg(0, 0, -0.00059),
291 angle2deg(0, 0, 0.001813),
292 };
293 double correction = poly(c, coef, nitems(coef));
294 return angle2deg(23, 26, 21.448) + correction;
295 }
296
297 /*
298 * Calculate the declination at moment $t of an object at latitude $beta
299 * and longitude $lambda.
300 * Ref: Sec.(14.4), Eq.(14.29)
301 */
302 double
declination(double t,double beta,double lambda)303 declination(double t, double beta, double lambda)
304 {
305 double epsilon = obliquity(t);
306 return arcsin_deg(sin_deg(beta) * cos_deg(epsilon) +
307 cos_deg(beta) * sin_deg(epsilon) * sin_deg(lambda));
308 }
309
310 /*
311 * Calculate the right ascension at moment $t of an object at latitude $beta
312 * and longitude $lambda.
313 * Ref: Sec.(14.4), Eq.(14.30)
314 */
315 double
right_ascension(double t,double beta,double lambda)316 right_ascension(double t, double beta, double lambda)
317 {
318 double epsilon = obliquity(t);
319 double x = cos_deg(lambda);
320 double y = (sin_deg(lambda) * cos_deg(epsilon) -
321 tan_deg(beta) * sin_deg(epsilon));
322 return arctan_deg(y, x);
323 }
324
325 /*
326 * Calculate the refraction angle (in degrees) at a location of elevation
327 * $elevation.
328 * Ref: Sec.(14.7), Eq.(14.75)
329 */
330 double
refraction(double elevation)331 refraction(double elevation)
332 {
333 double h = (elevation > 0) ? elevation : 0.0;
334 double radius = 6.372e6; /* Earth radius in meters */
335 /* depression of visible horizon */
336 double dip = arccos_deg(radius / (radius + h));
337 /* depression contributed by an elevation of $h */
338 double dip2 = (19.0/3600.0) * sqrt(h);
339 /* average effect of refraction */
340 double avg = 34.0 / 60.0;
341
342 return avg + dip + dip2;
343 }
344
345 /*
346 * Determine the day of year of the fixed date $rd.
347 */
348 int
dayofyear_from_fixed(int rd)349 dayofyear_from_fixed(int rd)
350 {
351 int year = gregorian_year_from_fixed(rd);
352 struct date date = { year - 1, 12, 31 };
353
354 return rd - fixed_from_gregorian(&date);
355 }
356
357 /*
358 * Format the given time $t to 'HH:MM:SS' style.
359 */
360 int
format_time(char * buf,size_t size,double t)361 format_time(char *buf, size_t size, double t)
362 {
363 int hh, mm, ss, i;
364
365 t -= floor(t);
366 i = (int)round(t * 24*60*60);
367
368 hh = i / (60*60);
369 i %= 60*60;
370 mm = i / 60;
371 ss = i % 60;
372
373 return snprintf(buf, size, "%02d:%02d:%02d", hh, mm, ss);
374 }
375
376 /*
377 * Format the given timezone (in fraction of days) $zone to '[+-]HH:MM' style.
378 */
379 int
format_zone(char * buf,size_t size,double zone)380 format_zone(char *buf, size_t size, double zone)
381 {
382 bool positive;
383 int hh, mm, i;
384
385 positive = (zone >= 0);
386 i = (int)round(fabs(zone) * 24*60);
387 hh = i / 60;
388 mm = i % 60;
389
390 return snprintf(buf, size, "%c%02d:%02d",
391 (positive ? '+' : '-'), hh, mm);
392 }
393
394 /*
395 * Format the location to style: 'dd°mm'ss" [NS], dd°mm'ss" [EW], mm.m m'
396 */
397 int
format_location(char * buf,size_t size,const struct location * loc)398 format_location(char *buf, size_t size, const struct location *loc)
399 {
400 int d1, d2, m1, m2, s1, s2, i;
401 bool north, east;
402
403 north = (loc->latitude >= 0);
404 i = (int)round(fabs(loc->latitude) * 60*60);
405 d1 = i / (60*60);
406 i %= 60*60;
407 m1 = i / 60;
408 s1 = i % 60;
409
410 east = (loc->longitude >= 0);
411 i = (int)round(fabs(loc->longitude) * 60*60);
412 d2 = i / (60*60);
413 i %= 60*60;
414 m2 = i / 60;
415 s2 = i % 60;
416
417 return snprintf(buf, size, "%d°%d'%d\" %c, %d°%d'%d\" %c, %.1lf m",
418 d1, m1, s1, (north ? 'N' : 'S'),
419 d2, m2, s2, (east ? 'E' : 'W'),
420 loc->elevation);
421 }
422