kwin/colorcorrection/suncalc.cpp

/*
    KWin - the KDE window manager
    This file is part of the KDE project.

    SPDX-FileCopyrightText: 2017 Roman Gilg <subdiff@gmail.com>

    SPDX-License-Identifier: GPL-2.0-or-later
*/
#include "suncalc.h"
#include "constants.h"

#include <QDateTime>
#include <QTimeZone>
#include <QtMath>

namespace KWin {
namespace ColorCorrect {

#define TWILIGHT_NAUT       -12.0
#define TWILIGHT_CIVIL       -6.0
#define SUN_RISE_SET         -0.833
#define SUN_HIGH              2.0

static QTime convertToLocalTime(const QDateTime &when, const QTime &utcTime)
{
    const QTimeZone timeZone = QTimeZone::systemTimeZone();
    const int utcOffset = timeZone.offsetFromUtc(when);
    return utcTime.addSecs(utcOffset);
}

QPair<QDateTime, QDateTime> calculateSunTimings(const QDateTime &dateTime, double latitude, double longitude, bool morning)
{
    // calculations based on https://aa.quae.nl/en/reken/zonpositie.html
    // accuracy: +/- 5min

    // positioning
    const double rad = M_PI / 180.;
    const double earthObliquity = 23.4397; // epsilon

    const double lat = latitude;    // phi
    const double lng = -longitude;  // lw

    // times
    const QDateTime utcDateTime = dateTime.toUTC();
    const double juPrompt = utcDateTime.date().toJulianDay();   // J
    const double ju2000 = 2451545.;                 // J2000

    // geometry
    auto mod360 = [](double number) -> double {
        return std::fmod(number, 360.);
    };

    auto sin = [&rad](double angle) -> double {
        return std::sin(angle * rad);
    };
    auto cos = [&rad](double angle) -> double {
        return std::cos(angle * rad);
    };
    auto asin = [&rad](double val) -> double {
        return std::asin(val) / rad;
    };
    auto acos = [&rad](double val) -> double {
        return std::acos(val) / rad;
    };

    auto anomaly = [&](const double date) -> double { // M
        return mod360(357.5291 + 0.98560028 * (date - ju2000));
    };

    auto center = [&sin](double anomaly) -> double { // C
        return 1.9148 * sin(anomaly) + 0.02 * sin(2 * anomaly) + 0.0003 * sin(3 * anomaly);
    };

    auto ecliptLngMean = [](double anom) -> double { // Mean ecliptical longitude L_sun = Mean Anomaly + Perihelion + 180°
        return anom + 282.9372; // anom + 102.9372 + 180°
    };

    auto ecliptLng = [&](double anom) -> double { // lambda = L_sun + C
        return ecliptLngMean(anom) + center(anom);
    };

    auto declination = [&](const double date) -> double { // delta
        const double anom = anomaly(date);
        const double eclLng = ecliptLng(anom);

        return mod360(asin(sin(earthObliquity) * sin(eclLng)));
    };

    // sun hour angle at specific angle
    auto hourAngle = [&](const double date, double angle) -> double { // H_t
        const double decl = declination(date);
        const double ret0 = (sin(angle) - sin(lat) * sin(decl)) / (cos(lat) * cos(decl));

        double ret = mod360(acos( ret0 ));
        if (180. < ret) {
            ret = ret - 360.;
        }
        return ret;
    };

    /*
     * Sun positions
     */

    // transit is at noon
    auto getTransit = [&](const double date) -> double {  // Jtransit
        const double juMeanSolTime = juPrompt - ju2000 - 0.0009 - lng / 360.;       // n_x = J - J_2000 - J_0 - l_w / 360°
        const double juTrEstimate = date + qRound64(juMeanSolTime) - juMeanSolTime; // J_x = J + n - n_x
        const double anom = anomaly(juTrEstimate);                                  // M
        const double eclLngM = ecliptLngMean(anom);                                 // L_sun

        return juTrEstimate + 0.0053 * sin(anom) - 0.0068 * sin(2 * eclLngM);
    };

    auto getSunMorning = [&hourAngle](const double angle, const double transit) -> double {
        return transit - hourAngle(transit, angle) / 360.;
    };

    auto getSunEvening = [&hourAngle](const double angle, const double transit) -> double {
        return transit + hourAngle(transit, angle) / 360.;
    };

    /*
     * Begin calculations
     */

    // noon - sun at the highest point
    const double juNoon = getTransit(juPrompt);

    double begin, end;
    if (morning) {
        begin = getSunMorning(TWILIGHT_CIVIL, juNoon);
        end = getSunMorning(SUN_HIGH, juNoon);
    } else {
        begin = getSunEvening(SUN_HIGH, juNoon);
        end = getSunEvening(TWILIGHT_CIVIL, juNoon);
    }
    // transform to QDateTime
    begin += 0.5;
    end += 0.5;

    QDateTime dateTimeBegin;
    QDateTime dateTimeEnd;

    if (!std::isnan(begin)) {
        const double dayFraction = begin - int(begin);
        const QTime utcTime = QTime::fromMSecsSinceStartOfDay(dayFraction * MSC_DAY);
        const QTime localTime = convertToLocalTime(dateTime, utcTime);
        dateTimeBegin = QDateTime(dateTime.date(), localTime);
    }

    if (!std::isnan(end)) {
        const double dayFraction = end - int(end);
        const QTime utcTime = QTime::fromMSecsSinceStartOfDay(dayFraction * MSC_DAY);
        const QTime localTime = convertToLocalTime(dateTime, utcTime);
        dateTimeEnd = QDateTime(dateTime.date(), localTime);
    }

    return { dateTimeBegin, dateTimeEnd };
}

}
}