usr.bin/calendar/basics.c

   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 }