Differential D5928 Diff 14725 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
		38		const qint64 ju2000 = 2451545; // J2000
		39
		40		// positioning
		41		const double rad = M_PI / 180;
		42		const double earthObliquity = 23.4397 * rad; // epsilon
		43		const double perihelionLng = 102.9372 * rad; // BigPi
		44
		45		const double lat = rad * latitude; // phi
		46		double lng = - rad * longitude; // lw
		47
		48		// times
		49		const double j0 = 0.0009; // J0
		50
		51		// substract 0.5, since the first full day ends at noon
		52		const double juPrompt = prompt.toJulianDay() /- 0.5/; // J
		53		double juCycle = juPrompt - ju2000 - j0 - lng / (2 * M_PI); // nX
		54
		55		// geometry
		56		auto modulo = [](double number, double base) {
		57		if (number < 0) {
		58		base = -base;
		59		}
		60		return number - ((int)(number / base)) * base;
		61		};
		62		auto observer = [lng, ju2000, rad, modulo](double date) { // theta
		63		double ret0 = 280.147 + 360.9856235 * (date - ju2000) - lng / rad;
		64		double ret = std::fmod(ret0, 360) * rad;
		65		return ret;
		66		};
		67
		68		auto anomaly = [rad,modulo](double date) { // M
		69		double ret = rad * (357.5291 + 0.98560028 * (date - ju2000));
		70		return modulo(ret, 2 * M_PI);
		71		};
		72
		73		auto center = [rad, modulo](double anomaly) { // C
		74		double ret = modulo(rad * (1.9148 * qSin(anomaly) + 0.02 * qSin(2 * anomaly) + 0.0003 * qSin(3 * anomaly)), 2 * M_PI);
		75		return ret;
		76		};
		77
		78		auto ecliptLngMean = [perihelionLng, modulo](double anom) { // Lsun = M + BigPi + 180°
		79		double ret = modulo(anom + perihelionLng + M_PI, 2 * M_PI);
		80		return ret;
		81		};
		82
		83		auto ecliptLng = [ecliptLngMean, center](double anom) { // lambda = Lsun + C
		84		double ret = ecliptLngMean(anom) + center(anom);
		85		return ret;
		86		};
		87
		88		auto declination = [earthObliquity, anomaly, ecliptLng, modulo](double date) { // delta
		89		double anom = anomaly(date);
		90		double eclLng = ecliptLng(anom);
		91		double ret = qAsin(qSin(earthObliquity) * qSin(eclLng));
		92		return ret;
		93		};
		94
		95		// TODO: This is a bit ineffcient, because we calculate anom and eclLng in declination as well
		96		auto ascension = [earthObliquity, anomaly, ecliptLng, declination,rad,modulo](double date) { // delta
		97		double anom = anomaly(date);
		98		double eclLng = ecliptLng(anom);
		99		double ret = modulo(qAcos(qCos(eclLng) / qCos(declination(date))), 2 * M_PI);
		100		return ret;
		101		};
		102
		103		// sun hour angle at noon (angle is 0)
		104		auto hourAngle = [observer, rad, ascension, modulo](double date) { // H
		105		double ret0 = observer(date) - ascension(date);
		106		double ret = modulo(ret0, 2 * M_PI);
		107		if (M_PI < ret) {
		108		ret = ret - 2 * M_PI;
		109		}
		110		return ret; // theta - alpha
		111		};
		112		// sun hour angle at specific angle
		113		auto hourAngleT = [lat, declination, modulo](double date, double angle) { // M + BigPi + 180°
		114		double decl = declination(date);
		115		double ret0 = (qSin(angle) - qSin(lat) * qSin(decl)) / (qCos(lat) * qCos(decl));
		116		double ret = modulo(qAcos( ret0 ), 2 * M_PI);
		117		if (M_PI < ret) {
		118		ret = ret - 2 * M_PI;
		119		}
		120		return ret;
		121		};
		122
		123		// sun position getters
		124		auto getTransit = [juCycle, anomaly, ecliptLngMean, hourAngle,hourAngleT](double date) { // Jtransit
		125		double estimate = date + qRound64(juCycle) - juCycle; // Jx
		126		double anom = anomaly(estimate); // M
		127
		128		double jTr = estimate + 0.0053 * qSin(anom) - 0.0068 * qSin(2 * ecliptLngMean(anom));
		129		// improving
		130		for (int i = 0; i < 5; i++) {
		131		double hA = hourAngle(jTr);
		132		jTr = jTr - (hA / (2 * M_PI));
		133		if (hA < 0.0001) {
		134		break;
		135		}
		136		}
		137		return jTr;
		138		};
		139
		140		auto getSunMorning = [hourAngle, hourAngleT](double angle, double transit) {
		141		double jRise = transit - hourAngleT(transit, angle) / (2 * M_PI);
		142
		143		// improving
		144		for (int i = 0; i < 5; i++) {
		145		double hA = hourAngle(jRise);
		146		double hAT = hourAngleT(jRise, angle);
		147
		148		jRise = jRise - (hA + hAT) /(2 * M_PI);
		149
		150		if (hA + hAT < 0.0001) {
		151		break;
		152		}
		153		}
		154		return jRise;
		155
		156		};
		157
		158		auto getSunEvening = [hourAngle, hourAngleT](double angle, double transit) {
		159		double jSet = transit + hourAngleT(transit, angle) / (2 * M_PI);
		160
		161		// improving
		162		for (int i = 0; i < 5; i++) {
		163		double hA = hourAngle(jSet);
		164		double hAT = hourAngleT(jSet, angle);
		165
		166		jSet = jSet - (hA - hAT) /(2 * M_PI);
		167		if (hA - hAT < 0.0001) {
		168		break;
		169		}
		170		}
		171		return jSet;
		172		};
		173
		174		/*
		175		* Begin calculations
		176		*/
		177
		178		// noon - sun at the highest point
		179		double juNoon = getTransit(juPrompt);
		180
		181		double begin, end;
		182		if (morning) {
		183		begin = getSunMorning(TWILIGHT_CIVIL * rad, juNoon);
		184		end = getSunMorning(SUN_HIGH * rad, juNoon);
		185		} else {
		186		begin = getSunEvening(SUN_HIGH * rad, juNoon);
		187		end = getSunEvening(TWILIGHT_CIVIL * rad, juNoon);
		188		}
		189		// transform to QDateTime
		190		begin += 0.5;
		191		end += 0.5;
		192
		193		QTime timeBegin, timeEnd;
		194
		195		if (std::isnan(begin)) {
		196		timeBegin = QTime();
		197		} else {
		198		double timePart = begin - (int)begin;
		199		timeBegin = QTime::fromMSecsSinceStartOfDay((int)( timePart * MSC_DAY ));
		200		}
		201		if (std::isnan(end)) {
		202		timeEnd = QTime();
		203		} else {
		204		double timePart = end - (int)end;
		205		timeEnd = QTime::fromMSecsSinceStartOfDay((int)( timePart * MSC_DAY ));
		206		}
		207
		208		return QPair<QTime,QTime> (timeBegin, timeEnd);
		209		}
		210
		211		}
		212		}