Differential D5931 Diff 23750 libcolorcorrect/suncalc.cpp

Changeset View

Standalone View

View Options

libcolorcorrect/suncalc.cpp

  • This file was added.
1​/********************************************************************
2​Copyright 2017 Roman Gilg <subdiff@gmail.com>
3​
4​This program is free software; you can redistribute it and/or modify
5​it under the terms of the GNU General Public License as published by
6​the Free Software Foundation; either version 2 of the License, or
7​(at your option) any later version.
8​
9​This program is distributed in the hope that it will be useful,
10​but WITHOUT ANY WARRANTY; without even the implied warranty of
11​MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
12​GNU General Public License for more details.
13​
14​You should have received a copy of the GNU General Public License
15​along with this program. If not, see <http://www.gnu.org/licenses/>.
16​*********************************************************************/
17​#include "suncalc.h"
18​#include "colorcorrectconstants.h"
19​
20​#include <QDateTime>
21​#include <qmath.h>
22​
23​namespace ColorCorrect {
24​
25​static const double TWILIGHT_NAUT = -12.0;
26​static const double TWILIGHT_CIVIL = -6.0;
27​static const double SUN_RISE_SET = -0.833;
28​static const double SUN_HIGH = 2.0;
29​
30​QVariantMap calculateSunTimings(double latitude, double longitude, bool morning)
31​{
32​ // calculations based on http://aa.quae.nl/en/reken/zonpositie.html
33​ // accuracy: +/- 5min
34​
35​ // positioning
36​ const double rad = M_PI / 180.;
37​ const double earthObliquity = 23.4397; // epsilon
38​
39​ const double lat = latitude; // phi
40​ const double lng = -longitude; // lw
41​
42​ // times
43​ QDate prompt = QDate::currentDate();
44​ const double juPrompt = prompt.toJulianDay(); // J
45​ const double ju2000 = 2451545.; // J2000
46​
47​ // geometry
48​ auto mod360 = [](double number) -> double {
49​ return std::fmod(number, 360.);
50​ };
51​
52​ auto sin = [&rad](double angle) -> double {
53​ return std::sin(angle * rad);
54​ };
55​ auto cos = [&rad](double angle) -> double {
56​ return std::cos(angle * rad);
57​ };
58​ auto asin = [&rad](double val) -> double {
59​ return std::asin(val) / rad;
60​ };
61​ auto acos = [&rad](double val) -> double {
62​ return std::acos(val) / rad;
63​ };
64​
65​ auto anomaly = [&](const double date) -> double { // M
66​ return mod360(357.5291 + 0.98560028 * (date - ju2000));
67​ };
68​
69​ auto center = [&sin](double anomaly) -> double { // C
70​ return 1.9148 * sin(anomaly) + 0.02 * sin(2 * anomaly) + 0.0003 * sin(3 * anomaly);
71​ };
72​
73​ auto ecliptLngMean = [](double anom) -> double { // Mean ecliptical longitude L_sun = Mean Anomaly + Perihelion + 180°
74​ return anom + 282.9372; // anom + 102.9372 + 180°
75​ };
76​
77​ auto ecliptLng = [&](double anom) -> double { // lambda = L_sun + C
78​ return ecliptLngMean(anom) + center(anom);
79​ };
80​
81​ auto declination = [&](const double date) -> double { // delta
82​ const double anom = anomaly(date);
83​ const double eclLng = ecliptLng(anom);
84​
85​ return mod360(asin(sin(earthObliquity) * sin(eclLng)));
86​ };
87​
88​ // sun hour angle at specific angle
89​ auto hourAngle = [&](const double date, double angle) -> double { // H_t
90​ const double decl = declination(date);
91​ const double ret0 = (sin(angle) - sin(lat) * sin(decl)) / (cos(lat) * cos(decl));
92​
93​ double ret = mod360(acos( ret0 ));
94​ if (180. < ret) {
95​ ret = ret - 360.;
96​ }
97​ return ret;
98​ };
99​
100​ /*
101​ * Sun positions
102​ */
103​
104​ // transit is at noon
105​ auto getTransit = [&](const double date) -> double { // Jtransit
106​ const double juMeanSolTime = juPrompt - ju2000 - 0.0009 - lng / 360.; // n_x = J - J_2000 - J_0 - l_w / 360°
107​ const double juTrEstimate = date + qRound64(juMeanSolTime) - juMeanSolTime; // J_x = J + n - n_x
108​ const double anom = anomaly(juTrEstimate); // M
109​ const double eclLngM = ecliptLngMean(anom); // L_sun
110​
111​ return juTrEstimate + 0.0053 * sin(anom) - 0.0068 * sin(2 * eclLngM);
112​ };
113​
114​ auto getSunMorning = [&hourAngle](const double angle, const double transit) -> double {
115​ return transit - hourAngle(transit, angle) / 360.;
116​ };
117​
118​ auto getSunEvening = [&hourAngle](const double angle, const double transit) -> double {
119​ return transit + hourAngle(transit, angle) / 360.;
120​ };
121​
122​ /*
123​ * Begin calculations
124​ */
125​
126​ // noon - sun at the highest point
127​ const double juNoon = getTransit(juPrompt);
128​
129​ double begin, end;
130​ if (morning) {
131​ begin = getSunMorning(TWILIGHT_CIVIL, juNoon);
132​ end = getSunMorning(SUN_HIGH, juNoon);
133​ } else {
134​ begin = getSunEvening(SUN_HIGH, juNoon);
135​ end = getSunEvening(TWILIGHT_CIVIL, juNoon);
136​ }
137​ // transform to QDateTime
138​ begin += 0.5;
139​ end += 0.5;
140​
141​ QTime timeBegin, timeEnd;
142​
143​ if (std::isnan(begin)) {
144​ timeBegin = QTime();
145​ } else {
146​ double timePart = begin - (int)begin;
147​ timeBegin = QTime::fromMSecsSinceStartOfDay((int)( timePart * MSC_DAY ));
148​ }
149​ if (std::isnan(end)) {
150​ timeEnd = QTime();
151​ } else {
152​ double timePart = end - (int)end;
153​ timeEnd = QTime::fromMSecsSinceStartOfDay((int)( timePart * MSC_DAY ));
154​ }
155​
156​ QDateTime dateTimeBegin(prompt, timeBegin, Qt::UTC);
157​ QDateTime dateTimeEnd(prompt, timeEnd, Qt::UTC);
158​
159​ QVariantMap map;
160​ map.insert("begin", dateTimeBegin.toLocalTime());
161​ map.insert("end", dateTimeEnd.toLocalTime());
162​ return map;
163​}
164​
165​QVariantMap SunCalc::getMorningTimings(double latitude, double longitude)
166​{
167​ return calculateSunTimings(latitude, longitude, true);
168​}
169​
170​QVariantMap SunCalc::getEveningTimings(double latitude, double longitude)
171​{
172​ return calculateSunTimings(latitude, longitude, false);
173​}
174​
175​}