12#include "sunriseset.h"
15static const double MaxLat = 89.99;
16static const double MaxLong = 179.99;
18using namespace KHolidays;
19using namespace SunRiseSet;
23 return (degree * PI) / 180.00;
28 return (rad * 180.0) / PI;
32static double calcTimeJulianCent(
int jday)
34 return (
double(jday) - 2451545.0) / 36525.0;
38static double calcMeanObliquityOfEcliptic(
double t)
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;
46static double calcObliquityCorrection(
double t)
48 double e0 = calcMeanObliquityOfEcliptic(t);
49 double omega = 125.04 - 1934.136 * t;
50 double e = e0 + 0.00256 * cos(
degToRad(omega));
55static double calcGeomMeanLongSun(
double t)
57 double L0 = 280.46646 + t * (36000.76983 + 0.0003032 * t);
68static double calcGeomMeanAnomalySun(
double t)
70 double M = 357.52911 + t * (35999.05029 - 0.0001537 * t);
75static double calcSunEqOfCenter(
double t)
77 double m = calcGeomMeanAnomalySun(t);
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)
89static double calcSunTrueLong(
double t)
91 double l0 = calcGeomMeanLongSun(t);
92 double c = calcSunEqOfCenter(t);
98static double calcSunApparentLong(
double t)
100 double o = calcSunTrueLong(t);
101 double omega = 125.04 - 1934.136 * t;
102 double lambda = o - 0.00569 - 0.00478 * sin(
degToRad(omega));
107static double calcSunDeclination(
double t)
109 double e = calcObliquityCorrection(t);
110 double lambda = calcSunApparentLong(t);
113 double theta =
radToDeg(asin(sint));
118static double calcEccentricityEarthOrbit(
double t)
120 double e = 0.016708634 - t * (0.000042037 + 0.0000001267 * t);
126static double calcEquationOfTime(
double t)
128 double epsilon = calcObliquityCorrection(t);
129 double l0 = calcGeomMeanLongSun(t);
130 double e = calcEccentricityEarthOrbit(t);
131 double m = calcGeomMeanAnomalySun(t);
133 double y = tan(
degToRad(epsilon) / 2.0);
136 double sin2l0 = sin(2.0 *
degToRad(l0));
138 double cos2l0 = cos(2.0 *
degToRad(l0));
139 double sin4l0 = sin(4.0 *
degToRad(l0));
140 double sin2m = sin(2.0 *
degToRad(m));
143 double Etime = (y * sin2l0
145 + 4.0 * e * y * sinm * cos2l0
146 - 0.5 * y * y * sin4l0
147 - 1.25 * e * e * sin2m);
153static constexpr const double Sunrise = -0.833;
154static constexpr const double CivilTwilight = -6.0;
156static constexpr const double Up = 1.0;
157static constexpr const double Down = -1.0;
160static double calcHourAngleSunrise(
double latitude,
double solarDec,
double sunHeight)
164 double HAarg = (cos(
degToRad(90.0 - sunHeight)) / (cos(latRad) * cos(sdRad)) - tan(latRad) * tan(sdRad));
165 double HA = acos(HAarg);
169static QTime calcSunEvent(
const QDate &date,
double latitude,
double longitude,
double sunHeight,
double direction)
171 latitude = qMax(qMin(MaxLat, latitude), -MaxLat);
172 longitude = qMax(qMin(MaxLong, longitude), -MaxLong);
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)) {
183 timeUTC = timeUTC.addSecs((720 - (4.0 * delta) - eqTime) * 60);
186 if (timeUTC.second() > 29) {
187 return QTime(timeUTC.hour(), timeUTC.minute()).
addSecs(60);
189 return QTime(timeUTC.hour(), timeUTC.minute());
194 return calcSunEvent(date, latitude, longitude, Sunrise, Up);
199 return calcSunEvent(date, latitude, longitude, Sunrise, Down);
204 return calcSunEvent(date, latitude, longitude, CivilTwilight, Up);
209 return calcSunEvent(date, latitude, longitude, CivilTwilight, Down);
215 const double t = calcTimeJulianCent(date.
toJulianDay());
216 const double solarDec = calcSunDeclination(t);
217 const double maxSolarZenithAngle = 180.0 - std::abs(latitude + solarDec);
219 return maxSolarZenithAngle <= 90.0 - Sunrise;
224 const double t = calcTimeJulianCent(date.
toJulianDay());
225 const double solarDec = calcSunDeclination(t);
226 const double minSolarZenithAngle = std::abs(latitude - solarDec);
228 return minSolarZenithAngle > 90.0 - Sunrise && minSolarZenithAngle <= 90.0 - CivilTwilight;
233 const double t = calcTimeJulianCent(date.
toJulianDay());
234 const double solarDec = calcSunDeclination(t);
235 const double minSolarZenithAngle = std::abs(latitude - solarDec);
237 return minSolarZenithAngle > 90.0 - CivilTwilight;
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