Differential D5928 Diff 23698 colorcorrection/suncalc.cpp

Changeset 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		}