Differential D5931 Diff 20970 libcolorcorrect/suncalc.cpp

Changeset View

Standalone View

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