File indexing completed on 2022-12-06 19:13:35

0001 /*
0002     This file is part of the kholidays library.
0003 
0004     Adapted from the Javascript found at https://www.esrl.noaa.gov/gmd/grad/solcalc
0005     which is in the public domain.
0006 
0007     SPDX-FileCopyrightText: 2012 Allen Winter <winter@kde.org>
0008 
0009     SPDX-License-Identifier: LGPL-2.0-or-later
0010 */
0011 
0012 #include "sunriseset.h"
0013 #include "util.h"
0014 #include <cmath>
0015 static const double MaxLat = 89.99;
0016 static const double MaxLong = 179.99;
0017 
0018 using namespace KHolidays;
0019 using namespace SunRiseSet;
0020 
0021 static double degToRad(double degree)
0022 {
0023     return (degree * PI) / 180.00;
0024 }
0025 
0026 static double radToDeg(double rad)
0027 {
0028     return (rad * 180.0) / PI;
0029 }
0030 
0031 // convert Julian Day to centuries since J2000.0.
0032 static double calcTimeJulianCent(int jday)
0033 {
0034     return (double(jday) - 2451545.0) / 36525.0;
0035 }
0036 
0037 // calculate the mean obliquity of the ecliptic (in degrees)
0038 static double calcMeanObliquityOfEcliptic(double t)
0039 {
0040     double seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * 0.001813));
0041     double e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
0042     return e0; // in degrees
0043 }
0044 
0045 // calculate the corrected obliquity of the ecliptic (in degrees)
0046 static double calcObliquityCorrection(double t)
0047 {
0048     double e0 = calcMeanObliquityOfEcliptic(t);
0049     double omega = 125.04 - 1934.136 * t;
0050     double e = e0 + 0.00256 * cos(degToRad(omega));
0051     return e; // in degrees
0052 }
0053 
0054 // calculate the Geometric Mean Longitude of the Sun (in degrees)
0055 static double calcGeomMeanLongSun(double t)
0056 {
0057     double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
0058     while (L0 > 360.0) {
0059         L0 -= 360.0;
0060     }
0061     while (L0 < 0.0) {
0062         L0 += 360.0;
0063     }
0064     return L0; // in degrees
0065 }
0066 
0067 // calculate the Geometric Mean Anomaly of the Sun (in degrees)
0068 static double calcGeomMeanAnomalySun(double t)
0069 {
0070     double M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
0071     return M; // in degrees
0072 }
0073 
0074 // calculate the equation of center for the sun (in degrees)
0075 static double calcSunEqOfCenter(double t)
0076 {
0077     double m = calcGeomMeanAnomalySun(t);
0078     double mrad = degToRad(m);
0079     double sinm = sin(mrad);
0080     double sin2m = sin(mrad + mrad);
0081     double sin3m = sin(mrad + mrad + mrad);
0082     double C = (sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) //
0083                 + sin2m * (0.019993 - 0.000101 * t) //
0084                 + sin3m * 0.000289);
0085     return C; // in degrees
0086 }
0087 
0088 // calculate the true longitude of the sun (in degrees)
0089 static double calcSunTrueLong(double t)
0090 {
0091     double l0 = calcGeomMeanLongSun(t);
0092     double c = calcSunEqOfCenter(t);
0093     double O = l0 + c;
0094     return O; // in degrees
0095 }
0096 
0097 // calculate the apparent longitude of the sun (in degrees)
0098 static double calcSunApparentLong(double t)
0099 {
0100     double o = calcSunTrueLong(t);
0101     double omega = 125.04 - 1934.136 * t;
0102     double lambda = o - 0.00569 - 0.00478 * sin(degToRad(omega));
0103     return lambda; // in degrees
0104 }
0105 
0106 // calculate the declination of the sun (in degrees)
0107 static double calcSunDeclination(double t)
0108 {
0109     double e = calcObliquityCorrection(t);
0110     double lambda = calcSunApparentLong(t);
0111 
0112     double sint = sin(degToRad(e)) * sin(degToRad(lambda));
0113     double theta = radToDeg(asin(sint));
0114     return theta; // in degrees
0115 }
0116 
0117 // calculate the eccentricity of earth's orbit (unitless)
0118 static double calcEccentricityEarthOrbit(double t)
0119 {
0120     double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
0121     return e; // unitless
0122 }
0123 
0124 // calculate the difference between true solar time and mean solar time
0125 //(output: equation of time in minutes of time)
0126 static double calcEquationOfTime(double t)
0127 {
0128     double epsilon = calcObliquityCorrection(t);
0129     double l0 = calcGeomMeanLongSun(t);
0130     double e = calcEccentricityEarthOrbit(t);
0131     double m = calcGeomMeanAnomalySun(t);
0132 
0133     double y = tan(degToRad(epsilon) / 2.0);
0134     y *= y;
0135 
0136     double sin2l0 = sin(2.0 * degToRad(l0));
0137     double sinm = sin(degToRad(m));
0138     double cos2l0 = cos(2.0 * degToRad(l0));
0139     double sin4l0 = sin(4.0 * degToRad(l0));
0140     double sin2m = sin(2.0 * degToRad(m));
0141 
0142     /* clang-format off */
0143     double Etime = (y * sin2l0
0144                     - 2.0 * e * sinm
0145                     + 4.0 * e * y * sinm * cos2l0
0146                     - 0.5 * y * y * sin4l0
0147                     - 1.25 * e * e * sin2m);
0148     /* clang-format on */
0149     return radToDeg(Etime) * 4.0; // in minutes of time
0150 }
0151 
0152 // sun positions (in degree relative to the horizon) for certain events
0153 static constexpr const double Sunrise = -0.833; // see https://en.wikipedia.org/wiki/Sunrise
0154 static constexpr const double CivilTwilight = -6.0; // see https://en.wikipedia.org/wiki/Twilight
0155 
0156 static constexpr const double Up = 1.0;
0157 static constexpr const double Down = -1.0;
0158 
0159 // calculate the hour angle of the sun at angle @p sunHeight for the latitude (in radians)
0160 static double calcHourAngleSunrise(double latitude, double solarDec, double sunHeight)
0161 {
0162     double latRad = degToRad(latitude);
0163     double sdRad = degToRad(solarDec);
0164     double HAarg = (cos(degToRad(90.0 - sunHeight)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad));
0165     double HA = acos(HAarg);
0166     return HA; // in radians (for sunset, use -HA)
0167 }
0168 
0169 static QTime calcSunEvent(const QDate &date, double latitude, double longitude, double sunHeight, double direction)
0170 {
0171     latitude = qMax(qMin(MaxLat, latitude), -MaxLat);
0172     longitude = qMax(qMin(MaxLong, longitude), -MaxLong);
0173 
0174     double t = calcTimeJulianCent(date.toJulianDay());
0175     double eqTime = calcEquationOfTime(t);
0176     double solarDec = calcSunDeclination(t);
0177     double hourAngle = direction * calcHourAngleSunrise(latitude, solarDec, sunHeight);
0178     double delta = longitude + radToDeg(hourAngle);
0179     if (std::isnan(delta)) {
0180         return {};
0181     }
0182     QTime timeUTC(0, 0);
0183     timeUTC = timeUTC.addSecs((720 - (4.0 * delta) - eqTime) * 60);
0184 
0185     // round to nearest minute
0186     if (timeUTC.second() > 29) {
0187         return QTime(timeUTC.hour(), timeUTC.minute()).addSecs(60);
0188     }
0189     return QTime(timeUTC.hour(), timeUTC.minute());
0190 }
0191 
0192 QTime KHolidays::SunRiseSet::utcSunrise(const QDate &date, double latitude, double longitude)
0193 {
0194     return calcSunEvent(date, latitude, longitude, Sunrise, Up);
0195 }
0196 
0197 QTime KHolidays::SunRiseSet::utcSunset(const QDate &date, double latitude, double longitude)
0198 {
0199     return calcSunEvent(date, latitude, longitude, Sunrise, Down);
0200 }
0201 
0202 QTime SunRiseSet::utcDawn(const QDate &date, double latitude, double longitude)
0203 {
0204     return calcSunEvent(date, latitude, longitude, CivilTwilight, Up);
0205 }
0206 
0207 QTime SunRiseSet::utcDusk(const QDate &date, double latitude, double longitude)
0208 {
0209     return calcSunEvent(date, latitude, longitude, CivilTwilight, Down);
0210 }
0211 
0212 // see https://en.wikipedia.org/wiki/Solar_zenith_angle
0213 bool SunRiseSet::isPolarDay(const QDate &date, double latitude)
0214 {
0215     const double t = calcTimeJulianCent(date.toJulianDay());
0216     const double solarDec = calcSunDeclination(t);
0217     const double maxSolarZenithAngle = 180.0 - std::abs(latitude + solarDec);
0218 
0219     return maxSolarZenithAngle <= 90.0 - Sunrise;
0220 }
0221 
0222 bool SunRiseSet::isPolarTwilight(const QDate &date, double latitude)
0223 {
0224     const double t = calcTimeJulianCent(date.toJulianDay());
0225     const double solarDec = calcSunDeclination(t);
0226     const double minSolarZenithAngle = std::abs(latitude - solarDec);
0227 
0228     return minSolarZenithAngle > 90.0 - Sunrise && minSolarZenithAngle <= 90.0 - CivilTwilight;
0229 }
0230 
0231 bool SunRiseSet::isPolarNight(const QDate &date, double latitude)
0232 {
0233     const double t = calcTimeJulianCent(date.toJulianDay());
0234     const double solarDec = calcSunDeclination(t);
0235     const double minSolarZenithAngle = std::abs(latitude - solarDec);
0236 
0237     return minSolarZenithAngle > 90.0 - CivilTwilight;
0238 }