KHolidays

sunriseset.cpp
1/*
2 This file is part of the kholidays library.
3
4 Adapted from the Javascript found at https://www.esrl.noaa.gov/gmd/grad/solcalc
5 which is in the public domain.
6
7 SPDX-FileCopyrightText: 2012 Allen Winter <winter@kde.org>
8
9 SPDX-License-Identifier: LGPL-2.0-or-later
10*/
11
12#include "sunriseset.h"
13#include "util.h"
14#include <cmath>
15static const double MaxLat = 89.99;
16static const double MaxLong = 179.99;
17
18using namespace KHolidays;
19using namespace SunRiseSet;
20
21static double degToRad(double degree)
22{
23 return (degree * PI) / 180.00;
24}
25
26static double radToDeg(double rad)
27{
28 return (rad * 180.0) / PI;
29}
30
31// convert Julian Day to centuries since J2000.0.
32static double calcTimeJulianCent(int jday)
33{
34 return (double(jday) - 2451545.0) / 36525.0;
35}
36
37// calculate the mean obliquity of the ecliptic (in degrees)
38static double calcMeanObliquityOfEcliptic(double t)
39{
40 double seconds = 21.448 - t * (46.8150 + t * (0.00059 - t * 0.001813));
41 double e0 = 23.0 + (26.0 + (seconds / 60.0)) / 60.0;
42 return e0; // in degrees
43}
44
45// calculate the corrected obliquity of the ecliptic (in degrees)
46static double calcObliquityCorrection(double t)
47{
48 double e0 = calcMeanObliquityOfEcliptic(t);
49 double omega = 125.04 - 1934.136 * t;
50 double e = e0 + 0.00256 * cos(degToRad(omega));
51 return e; // in degrees
52}
53
54// calculate the Geometric Mean Longitude of the Sun (in degrees)
55static double calcGeomMeanLongSun(double t)
56{
57 double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
58 while (L0 > 360.0) {
59 L0 -= 360.0;
60 }
61 while (L0 < 0.0) {
62 L0 += 360.0;
63 }
64 return L0; // in degrees
65}
66
67// calculate the Geometric Mean Anomaly of the Sun (in degrees)
68static double calcGeomMeanAnomalySun(double t)
69{
70 double M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
71 return M; // in degrees
72}
73
74// calculate the equation of center for the sun (in degrees)
75static double calcSunEqOfCenter(double t)
76{
77 double m = calcGeomMeanAnomalySun(t);
78 double mrad = degToRad(m);
79 double sinm = sin(mrad);
80 double sin2m = sin(mrad + mrad);
81 double sin3m = sin(mrad + mrad + mrad);
82 double C = (sinm * (1.914602 - t * (0.004817 + 0.000014 * t)) //
83 + sin2m * (0.019993 - 0.000101 * t) //
84 + sin3m * 0.000289);
85 return C; // in degrees
86}
87
88// calculate the true longitude of the sun (in degrees)
89static double calcSunTrueLong(double t)
90{
91 double l0 = calcGeomMeanLongSun(t);
92 double c = calcSunEqOfCenter(t);
93 double O = l0 + c;
94 return O; // in degrees
95}
96
97// calculate the apparent longitude of the sun (in degrees)
98static double calcSunApparentLong(double t)
99{
100 double o = calcSunTrueLong(t);
101 double omega = 125.04 - 1934.136 * t;
102 double lambda = o - 0.00569 - 0.00478 * sin(degToRad(omega));
103 return lambda; // in degrees
104}
105
106// calculate the declination of the sun (in degrees)
107static double calcSunDeclination(double t)
108{
109 double e = calcObliquityCorrection(t);
110 double lambda = calcSunApparentLong(t);
111
112 double sint = sin(degToRad(e)) * sin(degToRad(lambda));
113 double theta = radToDeg(asin(sint));
114 return theta; // in degrees
115}
116
117// calculate the eccentricity of earth's orbit (unitless)
118static double calcEccentricityEarthOrbit(double t)
119{
120 double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
121 return e; // unitless
122}
123
124// calculate the difference between true solar time and mean solar time
125//(output: equation of time in minutes of time)
126static double calcEquationOfTime(double t)
127{
128 double epsilon = calcObliquityCorrection(t);
129 double l0 = calcGeomMeanLongSun(t);
130 double e = calcEccentricityEarthOrbit(t);
131 double m = calcGeomMeanAnomalySun(t);
132
133 double y = tan(degToRad(epsilon) / 2.0);
134 y *= y;
135
136 double sin2l0 = sin(2.0 * degToRad(l0));
137 double sinm = sin(degToRad(m));
138 double cos2l0 = cos(2.0 * degToRad(l0));
139 double sin4l0 = sin(4.0 * degToRad(l0));
140 double sin2m = sin(2.0 * degToRad(m));
141
142 /* clang-format off */
143 double Etime = (y * sin2l0
144 - 2.0 * e * sinm
145 + 4.0 * e * y * sinm * cos2l0
146 - 0.5 * y * y * sin4l0
147 - 1.25 * e * e * sin2m);
148 /* clang-format on */
149 return radToDeg(Etime) * 4.0; // in minutes of time
150}
151
152// sun positions (in degree relative to the horizon) for certain events
153static constexpr const double Sunrise = -0.833; // see https://en.wikipedia.org/wiki/Sunrise
154static constexpr const double CivilTwilight = -6.0; // see https://en.wikipedia.org/wiki/Twilight
155
156static constexpr const double Up = 1.0;
157static constexpr const double Down = -1.0;
158
159// calculate the hour angle of the sun at angle @p sunHeight for the latitude (in radians)
160static double calcHourAngleSunrise(double latitude, double solarDec, double sunHeight)
161{
162 double latRad = degToRad(latitude);
163 double sdRad = degToRad(solarDec);
164 double HAarg = (cos(degToRad(90.0 - sunHeight)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad));
165 double HA = acos(HAarg);
166 return HA; // in radians (for sunset, use -HA)
167}
168
169static QTime calcSunEvent(const QDate &date, double latitude, double longitude, double sunHeight, double direction)
170{
171 latitude = qMax(qMin(MaxLat, latitude), -MaxLat);
172 longitude = qMax(qMin(MaxLong, longitude), -MaxLong);
173
174 double t = calcTimeJulianCent(date.toJulianDay());
175 double eqTime = calcEquationOfTime(t);
176 double solarDec = calcSunDeclination(t);
177 double hourAngle = direction * calcHourAngleSunrise(latitude, solarDec, sunHeight);
178 double delta = longitude + radToDeg(hourAngle);
179 if (std::isnan(delta)) {
180 return {};
181 }
182 QTime timeUTC(0, 0);
183 timeUTC = timeUTC.addSecs((720 - (4.0 * delta) - eqTime) * 60);
184
185 // round to nearest minute
186 if (timeUTC.second() > 29) {
187 return QTime(timeUTC.hour(), timeUTC.minute()).addSecs(60);
188 }
189 return QTime(timeUTC.hour(), timeUTC.minute());
190}
191
192QTime KHolidays::SunRiseSet::utcSunrise(const QDate &date, double latitude, double longitude)
193{
194 return calcSunEvent(date, latitude, longitude, Sunrise, Up);
195}
196
197QTime KHolidays::SunRiseSet::utcSunset(const QDate &date, double latitude, double longitude)
198{
199 return calcSunEvent(date, latitude, longitude, Sunrise, Down);
200}
201
202QTime SunRiseSet::utcDawn(const QDate &date, double latitude, double longitude)
203{
204 return calcSunEvent(date, latitude, longitude, CivilTwilight, Up);
205}
206
207QTime SunRiseSet::utcDusk(const QDate &date, double latitude, double longitude)
208{
209 return calcSunEvent(date, latitude, longitude, CivilTwilight, Down);
210}
211
212// see https://en.wikipedia.org/wiki/Solar_zenith_angle
213bool SunRiseSet::isPolarDay(const QDate &date, double latitude)
214{
215 const double t = calcTimeJulianCent(date.toJulianDay());
216 const double solarDec = calcSunDeclination(t);
217 const double maxSolarZenithAngle = 180.0 - std::abs(latitude + solarDec);
218
219 return maxSolarZenithAngle <= 90.0 - Sunrise;
220}
221
222bool SunRiseSet::isPolarTwilight(const QDate &date, double latitude)
223{
224 const double t = calcTimeJulianCent(date.toJulianDay());
225 const double solarDec = calcSunDeclination(t);
226 const double minSolarZenithAngle = std::abs(latitude - solarDec);
227
228 return minSolarZenithAngle > 90.0 - Sunrise && minSolarZenithAngle <= 90.0 - CivilTwilight;
229}
230
231bool SunRiseSet::isPolarNight(const QDate &date, double latitude)
232{
233 const double t = calcTimeJulianCent(date.toJulianDay());
234 const double solarDec = calcSunDeclination(t);
235 const double minSolarZenithAngle = std::abs(latitude - solarDec);
236
237 return minSolarZenithAngle > 90.0 - CivilTwilight;
238}
KHOLIDAYS_EXPORT bool isPolarTwilight(const QDate &date, double latitude)
Checks whether it is polar twilight on day date at latitude.
KHOLIDAYS_EXPORT QTime utcSunset(const QDate &date, double latitude, double longitude)
Compute the sunset time (UTC) for a date and Earth location.
KHOLIDAYS_EXPORT QTime utcSunrise(const QDate &date, double latitude, double longitude)
Compute the sunrise time (UTC) for a date and Earth location.
KHOLIDAYS_EXPORT bool isPolarDay(const QDate &date, double latitude)
Checks whether it is polar day on day date at latitude.
KHOLIDAYS_EXPORT bool isPolarNight(const QDate &date, double latitude)
Checks whether it is polar night on day date at latitude.
KHOLIDAYS_EXPORT QTime utcDawn(const QDate &date, double latitude, double longitude)
Compute the civil dawn time (UTC) for a date and Earth location.
KHOLIDAYS_EXPORT QTime utcDusk(const QDate &date, double latitude, double longitude)
Compute the civil dawn time (UTC) for a date and Earth location.
constexpr double radToDeg(double rad)
constexpr double degToRad(double deg)
qint64 toJulianDay() const const
QTime addSecs(int s) const const
This file is part of the KDE documentation.
Documentation copyright © 1996-2024 The KDE developers.
Generated on Tue Mar 26 2024 11:15:37 by doxygen 1.10.0 written by Dimitri van Heesch, © 1997-2006

KDE's Doxygen guidelines are available online.