/*- * SPDX-License-Identifier: BSD-3-Clause * * Copyright (c) 2019-2020 The DragonFly Project. All rights reserved. * * This code is derived from software contributed to The DragonFly Project * by Aaron LI * * Redistribution and use in source and binary forms, with or without * modification, are permitted provided that the following conditions * are met: * * 1. Redistributions of source code must retain the above copyright * notice, this list of conditions and the following disclaimer. * 2. Redistributions in binary form must reproduce the above copyright * notice, this list of conditions and the following disclaimer in * the documentation and/or other materials provided with the * distribution. * 3. Neither the name of The DragonFly Project nor the names of its * contributors may be used to endorse or promote products derived * from this software without specific, prior written permission. * * THIS SOFTWARE IS PROVIDED BY THE COPYRIGHT HOLDERS AND CONTRIBUTORS * ``AS IS'' AND ANY EXPRESS OR IMPLIED WARRANTIES, INCLUDING, BUT NOT * LIMITED TO, THE IMPLIED WARRANTIES OF MERCHANTABILITY AND FITNESS * FOR A PARTICULAR PURPOSE ARE DISCLAIMED. IN NO EVENT SHALL THE * COPYRIGHT HOLDERS OR CONTRIBUTORS BE LIABLE FOR ANY DIRECT, INDIRECT, * INCIDENTAL, SPECIAL, EXEMPLARY OR CONSEQUENTIAL DAMAGES (INCLUDING, * BUT NOT LIMITED TO, PROCUREMENT OF SUBSTITUTE GOODS OR SERVICES; * LOSS OF USE, DATA, OR PROFITS; OR BUSINESS INTERRUPTION) HOWEVER CAUSED * AND ON ANY THEORY OF LIABILITY, WHETHER IN CONTRACT, STRICT LIABILITY, * OR TORT (INCLUDING NEGLIGENCE OR OTHERWISE) ARISING IN ANY WAY OUT * OF THE USE OF THIS SOFTWARE, EVEN IF ADVISED OF THE POSSIBILITY OF * SUCH DAMAGE. * * Reference: * Calendrical Calculations, The Ultimate Edition (4th Edition) * by Edward M. Reingold and Nachum Dershowitz * 2018, Cambridge University Press */ #include #include #include #include #include "basics.h" #include "gregorian.h" #include "sun.h" #include "utils.h" /* * Time for the "mean sun" traval from one mean vernal equinox to the next. * Ref: Sec.(14.4), Eq.(14.31) */ const double mean_tropical_year = 365.242189; /* * Calculate the longitudinal nutation (in degrees) at moment $t. * Ref: Sec.(14.4), Eq.(14.34) */ double nutation(double t) { double c = julian_centuries(t); double coefsA[] = { 124.90, -1934.134, 0.002063 }; double coefsB[] = { 201.11, 72001.5377, 0.00057 }; double A = poly(c, coefsA, nitems(coefsA)); double B = poly(c, coefsB, nitems(coefsB)); return -0.004778 * sin_deg(A) - 0.0003667 * sin_deg(B); } /* * Calculate the aberration (in degrees) at moment $t. * Ref: Sec.(14.4), Eq.(14.35) */ double aberration(double t) { double c = julian_centuries(t); double A = 177.63 + 35999.01848 * c; return 0.0000974 * cos_deg(A) - 0.005575; } /* * Argument data used by 'solar_longitude()' for calculating the solar * longitude. * Ref: Sec.(14.4), Table(14.1) */ static const struct solar_longitude_arg { int x; double y; double z; } solar_longitude_data[] = { { 403406, 270.54861, 0.9287892 }, { 195207, 340.19128, 35999.1376958 }, { 119433, 63.91854, 35999.4089666 }, { 112392, 331.26220, 35998.7287385 }, { 3891, 317.843 , 71998.20261 }, { 2819, 86.631 , 71998.4403 }, { 1721, 240.052 , 36000.35726 }, { 660, 310.26 , 71997.4812 }, { 350, 247.23 , 32964.4678 }, { 334, 260.87 , -19.4410 }, { 314, 297.82 , 445267.1117 }, { 268, 343.14 , 45036.8840 }, { 242, 166.79 , 3.1008 }, { 234, 81.53 , 22518.4434 }, { 158, 3.50 , -19.9739 }, { 132, 132.75 , 65928.9345 }, { 129, 182.95 , 9038.0293 }, { 114, 162.03 , 3034.7684 }, { 99, 29.8 , 33718.148 }, { 93, 266.4 , 3034.448 }, { 86, 249.2 , -2280.773 }, { 78, 157.6 , 29929.992 }, { 72, 257.8 , 31556.493 }, { 68, 185.1 , 149.588 }, { 64, 69.9 , 9037.750 }, { 46, 8.0 , 107997.405 }, { 38, 197.1 , -4444.176 }, { 37, 250.4 , 151.771 }, { 32, 65.3 , 67555.316 }, { 29, 162.7 , 31556.080 }, { 28, 341.5 , -4561.540 }, { 27, 291.6 , 107996.706 }, { 27, 98.5 , 1221.655 }, { 25, 146.7 , 62894.167 }, { 24, 110.0 , 31437.369 }, { 21, 5.2 , 14578.298 }, { 21, 342.6 , -31931.757 }, { 20, 230.9 , 34777.243 }, { 18, 256.1 , 1221.999 }, { 17, 45.3 , 62894.511 }, { 14, 242.9 , -4442.039 }, { 13, 115.2 , 107997.909 }, { 13, 151.8 , 119.066 }, { 13, 285.3 , 16859.071 }, { 12, 53.3 , -4.578 }, { 10, 126.6 , 26895.292 }, { 10, 205.7 , -39.127 }, { 10, 85.9 , 12297.536 }, { 10, 146.1 , 90073.778 }, }; /* * Calculate the longitude (in degrees) of Sun at moment $t. * Ref: Sec.(14.4), Eq.(14.33) */ double solar_longitude(double t) { double c = julian_centuries(t); double sum = 0.0; const struct solar_longitude_arg *arg; for (size_t i = 0; i < nitems(solar_longitude_data); i++) { arg = &solar_longitude_data[i]; sum += arg->x * sin_deg(arg->y + arg->z * c); } double lambda = (282.7771834 + 36000.76953744 * c + 0.000005729577951308232 * sum); double ab = aberration(t); double nu = nutation(t); return mod_f(lambda + ab + nu, 360); } /* * Calculate the moment (in universal time) of the first time at or after * the given moment $t when the solar longitude will be $lambda degree. * Ref: Sec.(14.5), Eq.(14.36) */ double solar_longitude_atafter(double lambda, double t) { double rate = mean_tropical_year / 360.0; double lon = solar_longitude(t); double tau = t + rate * mod_f(lambda - lon, 360); /* estimate range (within 5 days) */ double a = (t > tau - 5) ? t : tau - 5; double b = tau + 5; return invert_angular(solar_longitude, lambda, a, b); } /* * Calculate the approximate moment at or before the given moment $t when * the solar longitude just exceeded the given degree $lambda. * Ref: Sec.(14.5), Eq.(14.42) */ double estimate_prior_solar_longitude(double lambda, double t) { double rate = mean_tropical_year / 360.0; /* first approximation */ double lon = solar_longitude(t); double tau = t - rate * mod_f(lon - lambda, 360); /* refine the estimate to within a day */ lon = solar_longitude(tau); double delta = mod3_f(lon - lambda, -180, 180); double t2 = tau - rate * delta; return (t < t2) ? t : t2; } /* * Calculate the geocentric altitude of Sun at moment $t and at location * ($latitude, $longitude), ignoring parallax and refraction. * Ref: Sec.(14.4), Eq.(14.41) */ double solar_altitude(double t, double latitude, double longitude) { double lambda = solar_longitude(t); double alpha = right_ascension(t, 0, lambda); double delta = declination(t, 0, lambda); double theta = sidereal_from_moment(t); double H = mod_f(theta + longitude - alpha, 360); double v = (sin_deg(latitude) * sin_deg(delta) + cos_deg(latitude) * cos_deg(delta) * cos_deg(H)); return mod3_f(arcsin_deg(v), -180, 180); } /* * Calculate the sine of angle between positions of Sun at local time $t * and when its depression angle is $alpha degrees at location ($latitude, * $longitude). * Ref: Sec.(14.7), Eq.(14.69) */ static double sine_offset(double t, double latitude, double longitude, double alpha) { double ut = t - longitude / 360.0; /* local -> universal time */ double lambda = solar_longitude(ut); double delta = declination(ut, 0, lambda); return (tan_deg(latitude) * tan_deg(delta) + sin_deg(alpha) / cos_deg(delta) / cos_deg(latitude)); } /* * Approximate the moment in local time near the given moment $t when * the depression angle of Sun is $alpha (negative if above horizon) at * location ($latitude, $longitude). If $morning is true, then searching * for the morning event; otherwise for the evening event. * NOTE: Return an NaN if the depression angle cannot be reached. * Ref: Sec.(14.7), Eq.(14.68) */ static double approx_depression_moment(double t, double latitude, double longitude, double alpha, bool morning) { double t2 = floor(t); /* midnight */ if (alpha < 0) t2 += 0.5; /* midday */ else if (morning == false) t2 += 1.0; /* next day */ double try = sine_offset(t, latitude, longitude, alpha); double value = ((fabs(try) > 1) ? sine_offset(t2, latitude, longitude, alpha) : try); if (fabs(value) > 1) { return NAN; } else { double offset = mod3_f(arcsin_deg(value) / 360.0, -0.5, 0.5); double t3 = floor(t); if (morning) t3 += 6.0/24.0 - offset; else t3 += 18.0/24.0 + offset; return local_from_apparent(t3, longitude); } } /* * Calculate the moment in local time near the given moment $tapprox when * the depression angle of Sun is $alpha (negative if above horizon) at * location ($latitude, $longitude). If $morning is true, then searching * for the morning event; otherwise for the evening event. * NOTE: Return an NaN if the depression angle cannot be reached. * Ref: Sec.(14.7), Eq.(14.70) */ static double depression_moment(double tapprox, double latitude, double longitude, double alpha, bool morning) { const double eps = 30.0 / 3600 / 24; /* accuracy of 30 seconds */ double t = approx_depression_moment(tapprox, latitude, longitude, alpha, morning); if (isnan(t)) return NAN; else if (fabs(t - tapprox) < eps) return t; else return depression_moment(t, latitude, longitude, alpha, morning); } /* * Calculate the moment of sunrise in standard time on fixed date $rd at * location $loc. * NOTE: Return an NaN if no sunrise. * Ref: Sec.(14.7), Eq.(14.72,14.76) */ double sunrise(int rd, const struct location *loc) { double sun_radius = 16.0 / 60.0; /* 16 arcminutes */ double alpha = refraction(loc->elevation) + sun_radius; double tapprox = (double)rd + (6.0/24.0); double lt = depression_moment(tapprox, loc->latitude, loc->longitude, alpha, true); if (isnan(lt)) return NAN; else return lt - loc->longitude / 360.0 + loc->zone; } /* * Calculate the moment of sunset in standard time on fixed date $rd at * location $loc. * NOTE: Return an NaN if no sunset. * Ref: Sec.(14.7), Eq.(14.74,14.77) */ double sunset(int rd, const struct location *loc) { double sun_radius = 16.0 / 60.0; /* 16 arcminutes */ double alpha = refraction(loc->elevation) + sun_radius; double tapprox = (double)rd + (18.0/24.0); double lt = depression_moment(tapprox, loc->latitude, loc->longitude, alpha, false); if (isnan(lt)) return NAN; else return lt - loc->longitude / 360.0 + loc->zone; } /**************************************************************************/ /* Equinoxes and solstices */ static const struct solar_event { const char *name; int longitude; /* longitude of Sun */ int month; /* month of the event */ } SOLAR_EVENTS[] = { { "March Equinox", 0, 3 }, { "June Solstice", 90, 6 }, { "September Equinox", 180, 9 }, { "December Solstice", 270, 12 }, }; /* * Print Sun information at the given moment $t (in standard time) * and events in the year. */ void show_sun_info(double t, const struct location *loc) { char buf[64]; int rd = (int)floor(t); /* * Sun position */ double t_u = t - loc->zone; /* universal time */ double lon = solar_longitude(t_u); double alt = solar_altitude(t_u, loc->latitude, loc->longitude); printf("Sun position: %.4lf° (longitude), %.4lf° (altitude)\n", lon, alt); /* * Sun rise and set */ double moments[2] = { sunrise(rd, loc), sunset(rd, loc) }; const char *names[2] = { "Sunrise", "Sunset" }; for (size_t i = 0; i < nitems(moments); i++) { if (isnan(moments[i])) snprintf(buf, sizeof(buf), "(null)"); else format_time(buf, sizeof(buf), moments[i]); printf("%-7s: %s\n", names[i], buf); } /* * Equinoxes and solstices */ const struct solar_event *event; int lambda, day_approx; int year = gregorian_year_from_fixed(rd); struct date date = { year, 1, 1 }; printf("\nSolar events in year %d:\n", year); for (size_t i = 0; i < nitems(SOLAR_EVENTS); i++) { event = &SOLAR_EVENTS[i]; lambda = event->longitude; date.month = event->month; date.day = 1; day_approx = fixed_from_gregorian(&date); t = solar_longitude_atafter(lambda, day_approx) + loc->zone; gregorian_from_fixed((int)floor(t), &date); format_time(buf, sizeof(buf), t); printf("%-17s: %3d°, %d-%02d-%02d %s\n", event->name, lambda, date.year, date.month, date.day, buf); } }