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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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 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