plasma-workspace/dataengines/time/solarsystem.cpp

/*
    SPDX-FileCopyrightText: 2009 Petri Damsten <damu@iki.fi>

    SPDX-License-Identifier: GPL-2.0-or-later
*/

#include "solarsystem.h"
#include <QList>
#include <math.h>

/*
 *   Mathematics, ideas, public domain code used for these classes from:
 *   https://www.stjarnhimlen.se/comp/tutorial.html
 *   https://www.stjarnhimlen.se/comp/riset.html
 *   https://www.esrl.noaa.gov/gmd/grad/solcalc/azel.html
 *   https://www.esrl.noaa.gov/gmd/grad/solcalc/sunrise.html
 *   http://web.archive.org/web/20080309162302/http://bodmas.org/astronomy/riset.html
 *   moontool.c by John Walker
 *   Wikipedia
 */

Sun::Sun()
    : SolarSystemObject()
{
}

void Sun::calcForDateTime(const QDateTime &local, int offset)
{
    SolarSystemObject::calcForDateTime(local, offset);

    N = 0.0;
    i = 0.0;
    w = rev(282.9404 + 4.70935E-5 * m_day);
    a = 1.0;
    e = rev(0.016709 - 1.151E-9 * m_day);
    M = rev(356.0470 + 0.9856002585 * m_day);

    calc();
}

void Sun::rotate(double *y, double *z)
{
    *y *= cosd(m_obliquity);
    *z *= sind(m_obliquity);
}

Moon::Moon(Sun *sun)
    : m_sun(sun)
{
}

void Moon::calcForDateTime(const QDateTime &local, int offset)
{
    if (m_sun->dateTime() != local) {
        m_sun->calcForDateTime(local, offset);
    }

    SolarSystemObject::calcForDateTime(local, offset);

    N = rev(125.1228 - 0.0529538083 * m_day);
    i = 5.1454;
    w = rev(318.0634 + 0.1643573223 * m_day);
    a = 60.2666;
    e = 0.054900;
    M = rev(115.3654 + 13.0649929509 * m_day);

    calc();
}

bool Moon::calcPerturbations(double *lo, double *la, double *r)
{
    double Ms = m_sun->meanAnomaly();
    double D = L - m_sun->meanLongitude();
    double F = L - N;
    // clang-format off
    *lo +=  -1.274 * sind(M - 2 * D)
            +0.658 * sind(2 * D)
            -0.186 * sind(Ms)
            -0.059 * sind(2 * M - 2 * D)
            -0.057 * sind(M - 2 * D + Ms)
            +0.053 * sind(M + 2 * D)
            +0.046 * sind(2 * D - Ms)
            +0.041 * sind(M - Ms)
            -0.035 * sind(D)
            -0.031 * sind(M + Ms)
            -0.015 * sind(2 * F - 2 * D)
            +0.011 * sind(M - 4 * D);
    *la +=  -0.173 * sind(F - 2 * D)
            -0.055 * sind(M - F - 2 * D)
            -0.046 * sind(M + F - 2 * D)
            +0.033 * sind(F + 2 * D)
            +0.017 * sind(2 * M + F);
    *r +=   -0.58 * cosd(M - 2 * D)
            -0.46 * cosd(2 * D);
    // clang-format on
    return true;
}

void Moon::topocentricCorrection(double *RA, double *dec)
{
    double HA = rev(siderealTime() - *RA);
    double gclat = m_latitude - 0.1924 * sind(2 * m_latitude);
    double rho = 0.99833 + 0.00167 * cosd(2 * m_latitude);
    double mpar = asind(1 / rad);
    double g = atand(tand(gclat) / cosd(HA));

    *RA -= mpar * rho * cosd(gclat) * sind(HA) / cosd(*dec);
    *dec -= mpar * rho * sind(gclat) * sind(g - *dec) / sind(g);
}

double Moon::phase()
{
    return rev(m_eclipticLongitude - m_sun->lambda());
}

void Moon::rotate(double *y, double *z)
{
    double t = *y;
    *y = t * cosd(m_obliquity) - *z * sind(m_obliquity);
    *z = t * sind(m_obliquity) + *z * cosd(m_obliquity);
}

void SolarSystemObject::calc()
{
    double x, y, z;
    double la, r;

    L = rev(N + w + M);
    double E0 = 720.0;
    double E = M + (180.0 / M_PI) * e * sind(M) * (1.0 + e * cosd(M));
    for (int j = 0; fabs(E0 - E) > 0.005 && j < 10; ++j) {
        E0 = E;
        E = E0 - (E0 - (180.0 / M_PI) * e * sind(E0) - M) / (1 - e * cosd(E0));
    }
    x = a * (cosd(E) - e);
    y = a * sind(E) * sqrt(1.0 - e * e);
    r = sqrt(x * x + y * y);
    double v = rev(atan2d(y, x));
    m_lambda = rev(v + w);
    x = r * (cosd(N) * cosd(m_lambda) - sind(N) * sind(m_lambda) * cosd(i));
    y = r * (sind(N) * cosd(m_lambda) + cosd(N) * sind(m_lambda) * cosd(i));
    z = r * sind(m_lambda);
    if (!qFuzzyCompare(i, 0.0)) {
        z *= sind(i);
    }
    toSpherical(x, y, z, &m_eclipticLongitude, &la, &r);
    if (calcPerturbations(&m_eclipticLongitude, &la, &r)) {
        toRectangular(m_eclipticLongitude, la, r, &x, &y, &z);
    }
    rotate(&y, &z);
    toSpherical(x, y, z, &RA, &dec, &rad);
    topocentricCorrection(&RA, &dec);

    HA = rev(siderealTime() - RA);
    x = cosd(HA) * cosd(dec) * sind(m_latitude) - sind(dec) * cosd(m_latitude);
    y = sind(HA) * cosd(dec);
    z = cosd(HA) * cosd(dec) * cosd(m_latitude) + sind(dec) * sind(m_latitude);
    m_azimuth = atan2d(y, x) + 180.0;
    m_altitude = asind(z);
}

double SolarSystemObject::siderealTime()
{
    double UT = m_utc.time().hour() + m_utc.time().minute() / 60.0 + m_utc.time().second() / 3600.0;
    double GMST0 = rev(282.9404 + 4.70935E-5 * m_day + 356.0470 + 0.9856002585 * m_day + 180.0);
    return GMST0 + UT * 15.0 + m_longitude;
}

void SolarSystemObject::calcForDateTime(const QDateTime &local, int offset)
{
    m_local = local;
    m_utc = local.addSecs(-offset);
    m_day = 367 * m_utc.date().year() - (7 * (m_utc.date().year() + ((m_utc.date().month() + 9) / 12))) / 4 + (275 * m_utc.date().month()) / 9
        + m_utc.date().day() - 730530;
    m_day += m_utc.time().hour() / 24.0 + m_utc.time().minute() / (24.0 * 60.0) + m_utc.time().second() / (24.0 * 60.0 * 60.0);
    m_obliquity = 23.4393 - 3.563E-7 * m_day;
}

SolarSystemObject::SolarSystemObject()
    : m_latitude(0.0)
    , m_longitude(0.0)
{
}

SolarSystemObject::~SolarSystemObject()
{
}

void SolarSystemObject::setPosition(double latitude, double longitude)
{
    m_latitude = latitude;
    m_longitude = longitude;
}

double SolarSystemObject::rev(double x)
{
    return x - floor(x / 360.0) * 360.0;
}

double SolarSystemObject::asind(double x)
{
    return asin(x) * 180.0 / M_PI;
}

double SolarSystemObject::sind(double x)
{
    return sin(x * M_PI / 180.0);
}

double SolarSystemObject::cosd(double x)
{
    return cos(x * M_PI / 180.0);
}

double SolarSystemObject::tand(double x)
{
    return tan(x * M_PI / 180.0);
}

double SolarSystemObject::atan2d(double y, double x)
{
    return atan2(y, x) * 180.0 / M_PI;
}

double SolarSystemObject::atand(double x)
{
    return atan(x) * 180.0 / M_PI;
}

void SolarSystemObject::toRectangular(double lo, double la, double r, double *x, double *y, double *z)
{
    *x = r * cosd(lo) * cosd(la);
    *y = r * sind(lo) * cosd(la);
    *z = r * sind(la);
}

void SolarSystemObject::toSpherical(double x, double y, double z, double *lo, double *la, double *r)
{
    *r = sqrt(x * x + y * y + z * z);
    *la = asind(z / *r);
    *lo = rev(atan2d(y, x));
}

QPair<double, double> SolarSystemObject::zeroPoints(QPointF p1, QPointF p2, QPointF p3)
{
    double a = ((p2.y() - p1.y()) * (p1.x() - p3.x()) + (p3.y() - p1.y()) * (p2.x() - p1.x()))
        / ((p1.x() - p3.x()) * (p2.x() * p2.x() - p1.x() * p1.x()) + (p2.x() - p1.x()) * (p3.x() * p3.x() - p1.x() * p1.x()));
    double b = ((p2.y() - p1.y()) - a * (p2.x() * p2.x() - p1.x() * p1.x())) / (p2.x() - p1.x());
    double c = p1.y() - a * p1.x() * p1.x() - b * p1.x();
    double discriminant = b * b - 4.0 * a * c;
    double z1 = -1.0, z2 = -1.0;

    if (discriminant >= 0.0) {
        z1 = (-b + sqrt(discriminant)) / (2 * a);
        z2 = (-b - sqrt(discriminant)) / (2 * a);
    }
    return QPair<double, double>(z1, z2);
}

QList<QPair<QDateTime, QDateTime>> SolarSystemObject::timesForAngles(const QList<double> &angles, const QDateTime &dt, int offset)
{
    QList<double> altitudes;
    QDate d = dt.date();
    QDateTime local(d, QTime(0, 0));
    for (int j = 0; j <= 25; ++j) {
        calcForDateTime(local, offset);
        altitudes.append(altitude());
        local = local.addSecs(60 * 60);
    }
    QList<QPair<QDateTime, QDateTime>> result;
    QTime rise, set;
    foreach (double angle, angles) {
        for (int j = 3; j <= 25; j += 2) {
            QPointF p1((j - 2) * 60 * 60, altitudes[j - 2] - angle);
            QPointF p2((j - 1) * 60 * 60, altitudes[j - 1] - angle);
            QPointF p3(j * 60 * 60, altitudes[j] - angle);
            QPair<double, double> z = zeroPoints(p1, p2, p3);
            if (z.first > p1.x() && z.first < p3.x()) {
                if (p1.y() < 0.0) {
                    rise = QTime(0, 0).addSecs(z.first);
                } else {
                    set = QTime(0, 0).addSecs(z.first);
                }
            }
            if (z.second > p1.x() && z.second < p3.x()) {
                if (p3.y() < 0.0) {
                    set = QTime(0, 0).addSecs(z.second);
                } else {
                    rise = QTime(0, 0).addSecs(z.second);
                }
            }
        }
        result.append(QPair<QDateTime, QDateTime>(QDateTime(d, rise), QDateTime(d, set)));
    }
    return result;
}

double SolarSystemObject::calcElevation()
{
    double refractionCorrection;

    if (m_altitude > 85.0) {
        refractionCorrection = 0.0;
    } else {
        double te = tand(m_altitude);
        if (m_altitude > 5.0) {
            refractionCorrection = 58.1 / te - 0.07 / (te * te * te) + 0.000086 / (te * te * te * te * te);
        } else if (m_altitude > -0.575) {
            refractionCorrection = 1735.0 + m_altitude * (-518.2 + m_altitude * (103.4 + m_altitude * (-12.79 + m_altitude * 0.711)));
        } else {
            refractionCorrection = -20.774 / te;
        }
        refractionCorrection = refractionCorrection / 3600.0;
    }
    return m_altitude + refractionCorrection;
}