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