usr.bin/calendar/sun.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 <math.h>
  43 #include <stdbool.h>
  44 #include <stddef.h>
  45 #include <stdio.h>
  46
  47 #include "basics.h"
  48 #include "gregorian.h"
  49 #include "sun.h"
  50 #include "utils.h"
  51
  52 /*
  53  * Time for the "mean sun" traval from one mean vernal equinox to the next.
  54  * Ref: Sec.(14.4), Eq.(14.31)
  55  */
  56 const double mean_tropical_year = 365.242189;
  57
  58 /*
  59  * Calculate the longitudinal nutation (in degrees) at moment $t.
  60  * Ref: Sec.(14.4), Eq.(14.34)
  61  */
  62 double
  63 nutation(double t)
  64 {
  65         double c = julian_centuries(t);
  66         double coefsA[] = { 124.90, -1934.134, 0.002063 };
  67         double coefsB[] = { 201.11, 72001.5377, 0.00057 };
  68         double A = poly(c, coefsA, nitems(coefsA));
  69         double B = poly(c, coefsB, nitems(coefsB));
  70         return -0.004778 * sin_deg(A) - 0.0003667 * sin_deg(B);
  71 }
  72
  73 /*
  74  * Calculate the aberration (in degrees) at moment $t.
  75  * Ref: Sec.(14.4), Eq.(14.35)
  76  */
  77 double
  78 aberration(double t)
  79 {
  80         double c = julian_centuries(t);
  81         double A = 177.63 + 35999.01848 * c;
  82         return 0.0000974 * cos_deg(A) - 0.005575;
  83 }
  84
  85 /*
  86  * Argument data used by 'solar_longitude()' for calculating the solar
  87  * longitude.
  88  * Ref: Sec.(14.4), Table(14.1)
  89  */
  90 static const struct solar_longitude_arg {
  91         int     x;
  92         double  y;
  93         double  z;
  94 } solar_longitude_data[] = {
  95         { 403406, 270.54861,      0.9287892 },
  96         { 195207, 340.19128,  35999.1376958 },
  97         { 119433,  63.91854,  35999.4089666 },
  98         { 112392, 331.26220,  35998.7287385 },
  99         {   3891, 317.843  ,  71998.20261   },
 100         {   2819,  86.631  ,  71998.4403    },
 101         {   1721, 240.052  ,  36000.35726   },
 102         {    660, 310.26   ,  71997.4812    },
 103         {    350, 247.23   ,  32964.4678    },
 104         {    334, 260.87   ,    -19.4410    },
 105         {    314, 297.82   , 445267.1117    },
 106         {    268, 343.14   ,  45036.8840    },
 107         {    242, 166.79   ,      3.1008    },
 108         {    234,  81.53   ,  22518.4434    },
 109         {    158,   3.50   ,    -19.9739    },
 110         {    132, 132.75   ,  65928.9345    },
 111         {    129, 182.95   ,   9038.0293    },
 112         {    114, 162.03   ,   3034.7684    },
 113         {     99,  29.8    ,  33718.148     },
 114         {     93, 266.4    ,   3034.448     },
 115         {     86, 249.2    ,  -2280.773     },
 116         {     78, 157.6    ,  29929.992     },
 117         {     72, 257.8    ,  31556.493     },
 118         {     68, 185.1    ,    149.588     },
 119         {     64,  69.9    ,   9037.750     },
 120         {     46,   8.0    , 107997.405     },
 121         {     38, 197.1    ,  -4444.176     },
 122         {     37, 250.4    ,    151.771     },
 123         {     32,  65.3    ,  67555.316     },
 124         {     29, 162.7    ,  31556.080     },
 125         {     28, 341.5    ,  -4561.540     },
 126         {     27, 291.6    , 107996.706     },
 127         {     27,  98.5    ,   1221.655     },
 128         {     25, 146.7    ,  62894.167     },
 129         {     24, 110.0    ,  31437.369     },
 130         {     21,   5.2    ,  14578.298     },
 131         {     21, 342.6    , -31931.757     },
 132         {     20, 230.9    ,  34777.243     },
 133         {     18, 256.1    ,   1221.999     },
 134         {     17,  45.3    ,  62894.511     },
 135         {     14, 242.9    ,  -4442.039     },
 136         {     13, 115.2    , 107997.909     },
 137         {     13, 151.8    ,    119.066     },
 138         {     13, 285.3    ,  16859.071     },
 139         {     12,  53.3    ,     -4.578     },
 140         {     10, 126.6    ,  26895.292     },
 141         {     10, 205.7    ,    -39.127     },
 142         {     10,  85.9    ,  12297.536     },
 143         {     10, 146.1    ,  90073.778     },
 144 };
 145
 146 /*
 147 * Calculate the longitude (in degrees) of Sun at moment $t.
 148  * Ref: Sec.(14.4), Eq.(14.33)
 149  */
 150 double
 151 solar_longitude(double t)
 152 {
 153         double c = julian_centuries(t);
 154
 155         double sum = 0.0;
 156         const struct solar_longitude_arg *arg;
 157         for (size_t i = 0; i < nitems(solar_longitude_data); i++) {
 158                 arg = &solar_longitude_data[i];
 159                 sum += arg->x * sin_deg(arg->y + arg->z * c);
 160         }
 161         double lambda = (282.7771834 + 36000.76953744 * c +
 162                          0.000005729577951308232 * sum);
 163
 164         double ab = aberration(t);
 165         double nu = nutation(t);
 166
 167         return mod_f(lambda + ab + nu, 360);
 168 }
 169
 170 /*
 171  * Calculate the moment (in universal time) of the first time at or after
 172  * the given moment $t when the solar longitude will be $lambda degree.
 173  * Ref: Sec.(14.5), Eq.(14.36)
 174  */
 175 double
 176 solar_longitude_atafter(double lambda, double t)
 177 {
 178         double rate = mean_tropical_year / 360.0;
 179         double lon = solar_longitude(t);
 180         double tau = t + rate * mod_f(lambda - lon, 360);
 181
 182         /* estimate range (within 5 days) */
 183         double a = (t > tau - 5) ? t : tau - 5;
 184         double b = tau + 5;
 185
 186         return invert_angular(solar_longitude, lambda, a, b);
 187 }
 188
 189 /*
 190  * Calculate the approximate moment at or before the given moment $t when
 191  * the solar longitude just exceeded the given degree $lambda.
 192  * Ref: Sec.(14.5), Eq.(14.42)
 193  */
 194 double
 195 estimate_prior_solar_longitude(double lambda, double t)
 196 {
 197         double rate = mean_tropical_year / 360.0;
 198
 199         /* first approximation */
 200         double lon = solar_longitude(t);
 201         double tau = t - rate * mod_f(lon - lambda, 360);
 202
 203         /* refine the estimate to within a day */
 204         lon = solar_longitude(tau);
 205         double delta = mod3_f(lon - lambda, -180, 180);
 206         double t2 = tau - rate * delta;
 207
 208         return (t < t2) ? t : t2;
 209 }
 210
 211 /*
 212  * Calculate the geocentric altitude of Sun at moment $t and at location
 213  * ($latitude, $longitude), ignoring parallax and refraction.
 214  * Ref: Sec.(14.4), Eq.(14.41)
 215  */
 216 double
 217 solar_altitude(double t, double latitude, double longitude)
 218 {
 219         double lambda = solar_longitude(t);
 220         double alpha = right_ascension(t, 0, lambda);
 221         double delta = declination(t, 0, lambda);
 222         double theta = sidereal_from_moment(t);
 223         double H = mod_f(theta + longitude - alpha, 360);
 224
 225         double v = (sin_deg(latitude) * sin_deg(delta) +
 226                     cos_deg(latitude) * cos_deg(delta) * cos_deg(H));
 227         return mod3_f(arcsin_deg(v), -180, 180);
 228 }
 229
 230 /*
 231  * Calculate the sine of angle between positions of Sun at local time $t
 232  * and when its depression angle is $alpha degrees at location ($latitude,
 233  * $longitude).
 234  * Ref: Sec.(14.7), Eq.(14.69)
 235  */
 236 static double
 237 sine_offset(double t, double latitude, double longitude, double alpha)
 238 {
 239         double ut = t - longitude / 360.0;  /* local -> universal time */
 240         double lambda = solar_longitude(ut);
 241         double delta = declination(ut, 0, lambda);
 242
 243         return (tan_deg(latitude) * tan_deg(delta) +
 244                 sin_deg(alpha) / cos_deg(delta) / cos_deg(latitude));
 245 }
 246
 247 /*
 248  * Approximate the moment in local time near the given moment $t when
 249  * the depression angle of Sun is $alpha (negative if above horizon) at
 250  * location ($latitude, $longitude).  If $morning is true, then searching
 251  * for the morning event; otherwise for the evening event.
 252  * NOTE: Return an NaN if the depression angle cannot be reached.
 253  * Ref: Sec.(14.7), Eq.(14.68)
 254  */
 255 static double
 256 approx_depression_moment(double t, double latitude, double longitude,
 257                          double alpha, bool morning)
 258 {
 259         double t2 = floor(t);  /* midnight */
 260         if (alpha < 0)
 261                 t2 += 0.5;  /* midday */
 262         else if (morning == false)
 263                 t2 += 1.0;  /* next day */
 264
 265         double try = sine_offset(t, latitude, longitude, alpha);
 266         double value = ((fabs(try) > 1) ?
 267                         sine_offset(t2, latitude, longitude, alpha) :
 268                         try);
 269
 270         if (fabs(value) > 1) {
 271                 return NAN;
 272         } else {
 273                 double offset = mod3_f(arcsin_deg(value) / 360.0, -0.5, 0.5);
 274                 double t3 = floor(t);
 275                 if (morning)
 276                         t3 += 6.0/24.0 - offset;
 277                 else
 278                         t3 += 18.0/24.0 + offset;
 279                 return local_from_apparent(t3, longitude);
 280         }
 281 }
 282
 283 /*
 284  * Calculate the moment in local time near the given moment $tapprox when
 285  * the depression angle of Sun is $alpha (negative if above horizon) at
 286  * location ($latitude, $longitude).  If $morning is true, then searching
 287  * for the morning event; otherwise for the evening event.
 288  * NOTE: Return an NaN if the depression angle cannot be reached.
 289  * Ref: Sec.(14.7), Eq.(14.70)
 290  */
 291 static double
 292 depression_moment(double tapprox, double latitude, double longitude,
 293                   double alpha, bool morning)
 294 {
 295         const double eps = 30.0 / 3600 / 24;  /* accuracy of 30 seconds */
 296         double t = approx_depression_moment(tapprox, latitude, longitude,
 297                                             alpha, morning);
 298         if (isnan(t))
 299                 return NAN;
 300         else if (fabs(t - tapprox) < eps)
 301                 return t;
 302         else
 303                 return depression_moment(t, latitude, longitude,
 304                                          alpha, morning);
 305 }
 306
 307 /*
 308  * Calculate the moment of sunrise in standard time on fixed date $rd at
 309  * location $loc.
 310  * NOTE: Return an NaN if no sunrise.
 311  * Ref: Sec.(14.7), Eq.(14.72,14.76)
 312  */
 313 double
 314 sunrise(int rd, const struct location *loc)
 315 {
 316         double sun_radius = 16.0 / 60.0;  /* 16 arcminutes */
 317         double alpha = refraction(loc->elevation) + sun_radius;
 318         double tapprox = (double)rd + (6.0/24.0);
 319         double lt = depression_moment(tapprox, loc->latitude, loc->longitude,
 320                                       alpha, true);
 321         if (isnan(lt))
 322                 return NAN;
 323         else
 324                 return lt - loc->longitude / 360.0 + loc->zone;
 325 }
 326
 327 /*
 328  * Calculate the moment of sunset in standard time on fixed date $rd at
 329  * location $loc.
 330  * NOTE: Return an NaN if no sunset.
 331  * Ref: Sec.(14.7), Eq.(14.74,14.77)
 332  */
 333 double
 334 sunset(int rd, const struct location *loc)
 335 {
 336         double sun_radius = 16.0 / 60.0;  /* 16 arcminutes */
 337         double alpha = refraction(loc->elevation) + sun_radius;
 338         double tapprox = (double)rd + (18.0/24.0);
 339         double lt = depression_moment(tapprox, loc->latitude, loc->longitude,
 340                                       alpha, false);
 341         if (isnan(lt))
 342                 return NAN;
 343         else
 344                 return lt - loc->longitude / 360.0 + loc->zone;
 345 }
 346
 347 /**************************************************************************/
 348
 349 /* Equinoxes and solstices */
 350 static const struct solar_event {
 351         const char      *name;
 352         int             longitude;  /* longitude of Sun */
 353         int             month;  /* month of the event */
 354 } SOLAR_EVENTS[] = {
 355         { "March Equinox",       0,  3 },
 356         { "June Solstice",      90,  6 },
 357         { "September Equinox", 180,  9 },
 358         { "December Solstice", 270, 12 },
 359 };
 360
 361 /*
 362  * Print Sun information at the given moment $t (in standard time)
 363  * and events in the year.
 364  */
 365 void
 366 show_sun_info(double t, const struct location *loc)
 367 {
 368         char buf[64];
 369         int rd = (int)floor(t);
 370
 371         /*
 372          * Sun position
 373          */
 374         double t_u = t - loc->zone;  /* universal time */
 375         double lon = solar_longitude(t_u);
 376         double alt = solar_altitude(t_u, loc->latitude, loc->longitude);
 377         printf("Sun position: %.4lf° (longitude), %.4lf° (altitude)\n",
 378                lon, alt);
 379
 380         /*
 381          * Sun rise and set
 382          */
 383         double moments[2] = { sunrise(rd, loc), sunset(rd, loc) };
 384         const char *names[2] = { "Sunrise", "Sunset" };
 385         for (size_t i = 0; i < nitems(moments); i++) {
 386                 if (isnan(moments[i]))
 387                         snprintf(buf, sizeof(buf), "(null)");
 388                 else
 389                         format_time(buf, sizeof(buf), moments[i]);
 390                 printf("%-7s: %s\n", names[i], buf);
 391         }
 392
 393         /*
 394          * Equinoxes and solstices
 395          */
 396         const struct solar_event *event;
 397         int lambda, day_approx;
 398         int year = gregorian_year_from_fixed(rd);
 399         struct date date = { year, 1, 1 };
 400
 401         printf("\nSolar events in year %d:\n", year);
 402         for (size_t i = 0; i < nitems(SOLAR_EVENTS); i++) {
 403                 event = &SOLAR_EVENTS[i];
 404                 lambda = event->longitude;
 405                 date.month = event->month;
 406                 date.day = 1;
 407                 day_approx = fixed_from_gregorian(&date);
 408                 t = solar_longitude_atafter(lambda, day_approx) + loc->zone;
 409                 gregorian_from_fixed((int)floor(t), &date);
 410                 format_time(buf, sizeof(buf), t);
 411                 printf("%-17s: %3d°, %d-%02d-%02d %s\n",
 412                        event->name, lambda,
 413                        date.year, date.month, date.day, buf);
 414         }
 415 }