File indexing completed on 2024-04-28 16:54:30
0001 /* 0002 SPDX-FileCopyrightText: 2017 Roman Gilg <subdiff@gmail.com> 0003 0004 SPDX-License-Identifier: GPL-2.0-or-later 0005 */ 0006 0007 #include "suncalc.h" 0008 #include "colorcorrectconstants.h" 0009 0010 #include <QDateTime> 0011 #include <qmath.h> 0012 0013 namespace ColorCorrect 0014 { 0015 static const double TWILIGHT_NAUT = -12.0; 0016 static const double TWILIGHT_CIVIL = -6.0; 0017 static const double SUN_RISE_SET = -0.833; 0018 static const double SUN_HIGH = 2.0; 0019 0020 QVariantMap calculateSunTimings(double latitude, double longitude, bool morning) 0021 { 0022 // calculations based on https://aa.quae.nl/en/reken/zonpositie.html 0023 // accuracy: +/- 5min 0024 0025 // positioning 0026 const double rad = M_PI / 180.; 0027 const double earthObliquity = 23.4397; // epsilon 0028 0029 const double lat = latitude; // phi 0030 const double lng = -longitude; // lw 0031 0032 // times 0033 QDate prompt = QDate::currentDate(); 0034 const double juPrompt = prompt.toJulianDay(); // J 0035 const double ju2000 = 2451545.; // J2000 0036 0037 // geometry 0038 auto mod360 = [](double number) -> double { 0039 return std::fmod(number, 360.); 0040 }; 0041 0042 auto sin = [&rad](double angle) -> double { 0043 return std::sin(angle * rad); 0044 }; 0045 auto cos = [&rad](double angle) -> double { 0046 return std::cos(angle * rad); 0047 }; 0048 auto asin = [&rad](double val) -> double { 0049 return std::asin(val) / rad; 0050 }; 0051 auto acos = [&rad](double val) -> double { 0052 return std::acos(val) / rad; 0053 }; 0054 0055 auto anomaly = [&](const double date) -> double { // M 0056 return mod360(357.5291 + 0.98560028 * (date - ju2000)); 0057 }; 0058 0059 auto center = [&sin](double anomaly) -> double { // C 0060 return 1.9148 * sin(anomaly) + 0.02 * sin(2 * anomaly) + 0.0003 * sin(3 * anomaly); 0061 }; 0062 0063 auto ecliptLngMean = [](double anom) -> double { // Mean ecliptical longitude L_sun = Mean Anomaly + Perihelion + 180° 0064 return anom + 282.9372; // anom + 102.9372 + 180° 0065 }; 0066 0067 auto ecliptLng = [&](double anom) -> double { // lambda = L_sun + C 0068 return ecliptLngMean(anom) + center(anom); 0069 }; 0070 0071 auto declination = [&](const double date) -> double { // delta 0072 const double anom = anomaly(date); 0073 const double eclLng = ecliptLng(anom); 0074 0075 return mod360(asin(sin(earthObliquity) * sin(eclLng))); 0076 }; 0077 0078 // sun hour angle at specific angle 0079 auto hourAngle = [&](const double date, double angle) -> double { // H_t 0080 const double decl = declination(date); 0081 const double ret0 = (sin(angle) - sin(lat) * sin(decl)) / (cos(lat) * cos(decl)); 0082 0083 double ret = mod360(acos(ret0)); 0084 if (180. < ret) { 0085 ret = ret - 360.; 0086 } 0087 return ret; 0088 }; 0089 0090 /* 0091 * Sun positions 0092 */ 0093 0094 // transit is at noon 0095 auto getTransit = [&](const double date) -> double { // Jtransit 0096 const double juMeanSolTime = juPrompt - ju2000 - 0.0009 - lng / 360.; // n_x = J - J_2000 - J_0 - l_w / 360° 0097 const double juTrEstimate = date + qRound64(juMeanSolTime) - juMeanSolTime; // J_x = J + n - n_x 0098 const double anom = anomaly(juTrEstimate); // M 0099 const double eclLngM = ecliptLngMean(anom); // L_sun 0100 0101 return juTrEstimate + 0.0053 * sin(anom) - 0.0068 * sin(2 * eclLngM); 0102 }; 0103 0104 auto getSunMorning = [&hourAngle](const double angle, const double transit) -> double { 0105 return transit - hourAngle(transit, angle) / 360.; 0106 }; 0107 0108 auto getSunEvening = [&hourAngle](const double angle, const double transit) -> double { 0109 return transit + hourAngle(transit, angle) / 360.; 0110 }; 0111 0112 /* 0113 * Begin calculations 0114 */ 0115 0116 // noon - sun at the highest point 0117 const double juNoon = getTransit(juPrompt); 0118 0119 double begin, end; 0120 if (morning) { 0121 begin = getSunMorning(TWILIGHT_CIVIL, juNoon); 0122 end = getSunMorning(SUN_HIGH, juNoon); 0123 } else { 0124 begin = getSunEvening(SUN_HIGH, juNoon); 0125 end = getSunEvening(TWILIGHT_CIVIL, juNoon); 0126 } 0127 // transform to QDateTime 0128 begin += 0.5; 0129 end += 0.5; 0130 0131 QTime timeBegin, timeEnd; 0132 0133 if (std::isnan(begin)) { 0134 timeBegin = QTime(); 0135 } else { 0136 double timePart = begin - (int)begin; 0137 timeBegin = QTime::fromMSecsSinceStartOfDay((int)(timePart * MSC_DAY)); 0138 } 0139 if (std::isnan(end)) { 0140 timeEnd = QTime(); 0141 } else { 0142 double timePart = end - (int)end; 0143 timeEnd = QTime::fromMSecsSinceStartOfDay((int)(timePart * MSC_DAY)); 0144 } 0145 0146 QDateTime dateTimeBegin(prompt, timeBegin, Qt::UTC); 0147 QDateTime dateTimeEnd(prompt, timeEnd, Qt::UTC); 0148 0149 QVariantMap map; 0150 map.insert("begin", dateTimeBegin.toLocalTime()); 0151 map.insert("end", dateTimeEnd.toLocalTime()); 0152 return map; 0153 } 0154 0155 QVariantMap SunCalc::getMorningTimings(double latitude, double longitude) 0156 { 0157 return calculateSunTimings(latitude, longitude, true); 0158 } 0159 0160 QVariantMap SunCalc::getEveningTimings(double latitude, double longitude) 0161 { 0162 return calculateSunTimings(latitude, longitude, false); 0163 } 0164 0165 }