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