GeoDataCoordinates.cpp

// SPDX-License-Identifier: LGPL-2.1-or-later
//
// SPDX-FileCopyrightText: 2004-2007 Torsten Rahn <tackat@kde.org>
// SPDX-FileCopyrightText: 2007-2008 Inge Wallin <ingwa@kde.org>
// SPDX-FileCopyrightText: 2008 Patrick Spendrin <ps_ml@gmx.de>
// SPDX-FileCopyrightText: 2011 Friedrich W. H. Kossebau <kossebau@kde.org>
// SPDX-FileCopyrightText: 2011 Bernhard Beschow <bbeschow@cs.tu-berlin.de>
// SPDX-FileCopyrightText: 2015 Alejandro Garcia Montoro <alejandro.garciamontoro@gmail.com>
//
 
#include "GeoDataCoordinates.h"
#include "GeoDataCoordinates_p.h"
#include "LonLatParser_p.h"
 
#include <QDataStream>
#include <QPointF>
#include <qmath.h>
 
#include "MarbleDebug.h"
#include "MarbleMath.h"
 
#include "Quaternion.h"
 
namespace Marble
{
 
const qreal GeoDataCoordinatesPrivate::sm_semiMajorAxis = 6378137.0;
const qreal GeoDataCoordinatesPrivate::sm_semiMinorAxis = 6356752.314;
const qreal GeoDataCoordinatesPrivate::sm_eccentricitySquared = 6.69437999013e-03;
const qreal GeoDataCoordinatesPrivate::sm_utmScaleFactor = 0.9996;
GeoDataCoordinates::Notation GeoDataCoordinates::s_notation = GeoDataCoordinates::DMS;
 
const GeoDataCoordinates GeoDataCoordinates::null = GeoDataCoordinates(0, 0, 0); // don't use default constructor!
 
GeoDataCoordinates::GeoDataCoordinates(qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit, int _detail)
    : d(new GeoDataCoordinatesPrivate(_lon, _lat, _alt, unit, _detail))
{
    d->ref.ref();
}
 
/* simply copy the d pointer
 * it will be replaced in the detach function instead
 */
GeoDataCoordinates::GeoDataCoordinates(const GeoDataCoordinates &other)
    : d(other.d)
{
    d->ref.ref();
}
 
/* simply copy null's d pointer
 * it will be replaced in the detach function
 */
GeoDataCoordinates::GeoDataCoordinates()
    : d(null.d)
{
    d->ref.ref();
}
 
/*
 * only delete the private d pointer if the number of references is 0
 * remember that all copies share the same d pointer!
 */
GeoDataCoordinates::~GeoDataCoordinates()
{
    delete d->m_q;
    d->m_q = nullptr;
 
    if (!d->ref.deref())
        delete d;
#ifdef DEBUG_GEODATA
//    mDebug() << "delete coordinates";
#endif
}
 
bool GeoDataCoordinates::isValid() const
{
    return d != null.d;
}
 
/*
 * if only one copy exists, return
 * else make a new private d pointer object and assign the values of the current
 * one to it
 * at the end, if the number of references thus reaches 0 delete it
 * this state shouldn't happen, but if it does, we have to clean up behind us.
 */
void GeoDataCoordinates::detach()
{
    delete d->m_q;
    d->m_q = nullptr;
 
    if (d->ref.fetchAndAddRelaxed(0) == 1) {
        return;
    }
 
    auto new_d = new GeoDataCoordinatesPrivate(*d);
 
    if (!d->ref.deref()) {
        delete d;
    }
 
    d = new_d;
    d->ref.ref();
}
 
/*
 * call detach() at the start of all non-static, non-const functions
 */
void GeoDataCoordinates::set(qreal _lon, qreal _lat, qreal _alt, GeoDataCoordinates::Unit unit)
{
    detach();
    d->m_altitude = _alt;
    switch (unit) {
    default:
    case Radian:
        d->m_lon = _lon;
        d->m_lat = _lat;
        break;
    case Degree:
        d->m_lon = _lon * DEG2RAD;
        d->m_lat = _lat * DEG2RAD;
        break;
    }
}
 
/*
 * call detach() at the start of all non-static, non-const functions
 */
void GeoDataCoordinates::setLongitude(qreal _lon, GeoDataCoordinates::Unit unit)
{
    detach();
    switch (unit) {
    default:
    case Radian:
        d->m_lon = _lon;
        break;
    case Degree:
        d->m_lon = _lon * DEG2RAD;
        break;
    }
}
 
/*
 * call detach() at the start of all non-static, non-const functions
 */
void GeoDataCoordinates::setLatitude(qreal _lat, GeoDataCoordinates::Unit unit)
{
    detach();
    switch (unit) {
    case Radian:
        d->m_lat = _lat;
        break;
    case Degree:
        d->m_lat = _lat * DEG2RAD;
        break;
    }
}
 
void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, GeoDataCoordinates::Unit unit) const
{
    switch (unit) {
    default:
    case Radian:
        lon = d->m_lon;
        lat = d->m_lat;
        break;
    case Degree:
        lon = d->m_lon * RAD2DEG;
        lat = d->m_lat * RAD2DEG;
        break;
    }
}
 
void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat) const
{
    lon = d->m_lon;
    lat = d->m_lat;
}
 
void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, qreal &alt, GeoDataCoordinates::Unit unit) const
{
    geoCoordinates(lon, lat, unit);
    alt = d->m_altitude;
}
 
void GeoDataCoordinates::geoCoordinates(qreal &lon, qreal &lat, qreal &alt) const
{
    lon = d->m_lon;
    lat = d->m_lat;
    alt = d->m_altitude;
}
 
qreal GeoDataCoordinates::longitude(GeoDataCoordinates::Unit unit) const
{
    switch (unit) {
    default:
    case Radian:
        return d->m_lon;
    case Degree:
        return d->m_lon * RAD2DEG;
    }
}
 
qreal GeoDataCoordinates::longitude() const
{
    return d->m_lon;
}
 
qreal GeoDataCoordinates::latitude(GeoDataCoordinates::Unit unit) const
{
    switch (unit) {
    default:
    case Radian:
        return d->m_lat;
    case Degree:
        return d->m_lat * RAD2DEG;
    }
}
 
qreal GeoDataCoordinates::latitude() const
{
    return d->m_lat;
}
 
// static
GeoDataCoordinates::Notation GeoDataCoordinates::defaultNotation()
{
    return s_notation;
}
 
// static
void GeoDataCoordinates::setDefaultNotation(GeoDataCoordinates::Notation notation)
{
    s_notation = notation;
}
 
// static
qreal GeoDataCoordinates::normalizeLon(qreal lon, GeoDataCoordinates::Unit unit)
{
    qreal halfCircle;
    if (unit == GeoDataCoordinates::Radian) {
        halfCircle = M_PI;
    } else {
        halfCircle = 180;
    }
 
    if (lon > halfCircle) {
        int cycles = (int)((lon + halfCircle) / (2 * halfCircle));
        return lon - (cycles * 2 * halfCircle);
    }
    if (lon < -halfCircle) {
        int cycles = (int)((lon - halfCircle) / (2 * halfCircle));
        return lon - (cycles * 2 * halfCircle);
    }
 
    return lon;
}
 
// static
qreal GeoDataCoordinates::normalizeLat(qreal lat, GeoDataCoordinates::Unit unit)
{
    qreal halfCircle;
    if (unit == GeoDataCoordinates::Radian) {
        halfCircle = M_PI;
    } else {
        halfCircle = 180;
    }
 
    if (lat > (halfCircle / 2.0)) {
        int cycles = (int)((lat + halfCircle) / (2 * halfCircle));
        qreal temp;
        if (cycles == 0) { // pi/2 < lat < pi
            temp = halfCircle - lat;
        } else {
            temp = lat - (cycles * 2 * halfCircle);
        }
        if (temp > (halfCircle / 2.0)) {
            return (halfCircle - temp);
        }
        if (temp < (-halfCircle / 2.0)) {
            return (-halfCircle - temp);
        }
        return temp;
    }
    if (lat < (-halfCircle / 2.0)) {
        int cycles = (int)((lat - halfCircle) / (2 * halfCircle));
        qreal temp;
        if (cycles == 0) {
            temp = -halfCircle - lat;
        } else {
            temp = lat - (cycles * 2 * halfCircle);
        }
        if (temp > (+halfCircle / 2.0)) {
            return (+halfCircle - temp);
        }
        if (temp < (-halfCircle / 2.0)) {
            return (-halfCircle - temp);
        }
        return temp;
    }
    return lat;
}
 
// static
void GeoDataCoordinates::normalizeLonLat(qreal &lon, qreal &lat, GeoDataCoordinates::Unit unit)
{
    qreal halfCircle;
    if (unit == GeoDataCoordinates::Radian) {
        halfCircle = M_PI;
    } else {
        halfCircle = 180;
    }
 
    if (lon > +halfCircle) {
        int cycles = (int)((lon + halfCircle) / (2 * halfCircle));
        lon = lon - (cycles * 2 * halfCircle);
    }
    if (lon < -halfCircle) {
        int cycles = (int)((lon - halfCircle) / (2 * halfCircle));
        lon = lon - (cycles * 2 * halfCircle);
    }
 
    if (lat > (+halfCircle / 2.0)) {
        int cycles = (int)((lat + halfCircle) / (2 * halfCircle));
        qreal temp;
        if (cycles == 0) { // pi/2 < lat < pi
            temp = halfCircle - lat;
        } else {
            temp = lat - (cycles * 2 * halfCircle);
        }
        if (temp > (+halfCircle / 2.0)) {
            lat = +halfCircle - temp;
        }
        if (temp < (-halfCircle / 2.0)) {
            lat = -halfCircle - temp;
        }
        lat = temp;
        if (lon > 0) {
            lon = -halfCircle + lon;
        } else {
            lon = halfCircle + lon;
        }
    }
    if (lat < (-halfCircle / 2.0)) {
        int cycles = (int)((lat - halfCircle) / (2 * halfCircle));
        qreal temp;
        if (cycles == 0) {
            temp = -halfCircle - lat;
        } else {
            temp = lat - (cycles * 2 * halfCircle);
        }
        if (temp > (+halfCircle / 2.0)) {
            lat = +halfCircle - temp;
        }
        if (temp < (-halfCircle / 2.0)) {
            lat = -halfCircle - temp;
        }
        lat = temp;
        if (lon > 0) {
            lon = -halfCircle + lon;
        } else {
            lon = halfCircle + lon;
        }
    }
    return;
}
 
GeoDataCoordinates GeoDataCoordinates::fromString(const QString &string, bool &successful)
{
    LonLatParser parser;
    successful = parser.parse(string);
    if (successful) {
        return {parser.lon(), parser.lat(), 0, GeoDataCoordinates::Degree};
    } else {
        return {};
    }
}
 
QString GeoDataCoordinates::toString() const
{
    return GeoDataCoordinates::toString(s_notation);
}
 
QString GeoDataCoordinates::toString(GeoDataCoordinates::Notation notation, int precision) const
{
    QString coordString;
 
    if (notation == GeoDataCoordinates::UTM) {
        int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat);
 
        // Handle lack of UTM zone number in the poles
        const QString zoneString = (zoneNumber > 0) ? QString::number(zoneNumber) : QString();
 
        QString bandString = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat);
 
        QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat), 'f', 2);
        QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat), 'f', 2);
 
        return QStringLiteral("%1%2 %3 m E, %4 m N").arg(zoneString, bandString, eastingString, northingString);
    } else {
        coordString = lonToString(d->m_lon, notation, Radian, precision) + QLatin1StringView(", ") + latToString(d->m_lat, notation, Radian, precision);
    }
 
    return coordString;
}
 
QString GeoDataCoordinates::lonToString(qreal lon, GeoDataCoordinates::Notation notation, GeoDataCoordinates::Unit unit, int precision, char format)
{
    if (notation == GeoDataCoordinates::UTM) {
        /**
         * @FIXME: UTM needs lon + lat to know zone number and easting
         * By now, this code returns the zone+easting of the point
         * (lon, equator), but this can differ a lot at different locations
         * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536
         */
 
        qreal lonRad = (unit == Radian) ? lon : lon * DEG2RAD;
 
        int zoneNumber = GeoDataCoordinatesPrivate::lonLatToZone(lonRad, 0);
 
        // Handle lack of UTM zone number in the poles
        QString result = (zoneNumber > 0) ? QString::number(zoneNumber) : QString();
 
        if (precision > 0) {
            QString eastingString = QString::number(GeoDataCoordinatesPrivate::lonLatToEasting(lonRad, 0), 'f', 2);
            result += QStringLiteral(" %1 m E").arg(eastingString);
        }
 
        return result;
    }
 
    QString weString = (lon < 0) ? tr("W") : tr("E");
 
    QString lonString;
 
    qreal lonDegF = (unit == Degree) ? fabs(lon) : fabs((qreal)(lon)*RAD2DEG);
 
    // Take care of -1 case
    precision = (precision < 0) ? 5 : precision;
 
    if (notation == DMS || notation == DM) {
        int lonDeg = (int)lonDegF;
        qreal lonMinF = 60 * (lonDegF - lonDeg);
        int lonMin = (int)lonMinF;
        qreal lonSecF = 60 * (lonMinF - lonMin);
        int lonSec = (int)lonSecF;
 
        // Adjustment for fuzziness (like 49.999999999999999999999)
        if (precision == 0) {
            lonDeg = qRound(lonDegF);
        } else if (precision <= 2) {
            lonMin = qRound(lonMinF);
        } else if (precision <= 4 && notation == DMS) {
            lonSec = qRound(lonSecF);
        } else {
            if (notation == DMS) {
                lonSec = lonSecF = qRound(lonSecF * qPow(10, precision - 4)) / qPow(10, precision - 4);
            } else {
                lonMin = lonMinF = qRound(lonMinF * qPow(10, precision - 2)) / qPow(10, precision - 2);
            }
        }
 
        if (lonSec > 59 && notation == DMS) {
            lonSec = lonSecF = 0;
            lonMin = lonMinF = lonMinF + 1;
        }
        if (lonMin > 59) {
            lonMin = lonMinF = 0;
            lonDeg = lonDegF = lonDegF + 1;
        }
 
        // Evaluate the string
        lonString = QString::fromUtf8("%1\xc2\xb0").arg(lonDeg, 3, 10, QLatin1Char(' '));
 
        if (precision == 0) {
            return lonString + weString;
        }
 
        if (notation == DMS || precision < 3) {
            lonString += QStringLiteral(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0'));
        }
 
        if (precision < 3) {
            return lonString + weString;
        }
 
        if (notation == DMS) {
            // Includes -1 case!
            if (precision < 5) {
                lonString += QStringLiteral(" %3\"").arg((double)lonSec, 2, 'f', 0, QLatin1Char('0'));
                return lonString + weString;
            }
 
            lonString += QStringLiteral(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
        } else {
            lonString += QStringLiteral(" %L3'").arg(lonMinF, precision + 1, 'f', precision - 2, QLatin1Char('0'));
        }
    } else if (notation == GeoDataCoordinates::Decimal) {
        lonString = QString::fromUtf8("%L1\xc2\xb0").arg(lonDegF, 4 + precision, format, precision, QLatin1Char(' '));
    } else if (notation == GeoDataCoordinates::Astro) {
        if (lon < 0) {
            lon += (unit == Degree) ? 360 : 2 * M_PI;
        }
 
        qreal lonHourF = (unit == Degree) ? fabs(lon / 15.0) : fabs((qreal)(lon / 15.0) * RAD2DEG);
        int lonHour = (int)lonHourF;
        qreal lonMinF = 60 * (lonHourF - lonHour);
        int lonMin = (int)lonMinF;
        qreal lonSecF = 60 * (lonMinF - lonMin);
        int lonSec = (int)lonSecF;
 
        // Adjustment for fuzziness (like 49.999999999999999999999)
        if (precision == 0) {
            lonHour = qRound(lonHourF);
        } else if (precision <= 2) {
            lonMin = qRound(lonMinF);
        } else if (precision <= 4) {
            lonSec = qRound(lonSecF);
        } else {
            lonSec = lonSecF = qRound(lonSecF * qPow(10, precision - 4)) / qPow(10, precision - 4);
        }
 
        if (lonSec > 59) {
            lonSec = lonSecF = 0;
            lonMin = lonMinF = lonMinF + 1;
        }
        if (lonMin > 59) {
            lonMin = lonMinF = 0;
            lonHour = lonHourF = lonHourF + 1;
        }
 
        // Evaluate the string
        lonString = QString::fromUtf8("%1h").arg(lonHour, 3, 10, QLatin1Char(' '));
 
        if (precision == 0) {
            return lonString;
        }
 
        lonString += QStringLiteral(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0'));
 
        if (precision < 3) {
            return lonString;
        }
 
        // Includes -1 case!
        if (precision < 5) {
            lonString += QStringLiteral(" %3\"").arg((double)lonSec, 2, 'f', 0, QLatin1Char('0'));
            return lonString;
        }
 
        lonString += QStringLiteral(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
        return lonString;
    }
 
    return lonString + weString;
}
 
QString GeoDataCoordinates::lonToString() const
{
    return GeoDataCoordinates::lonToString(d->m_lon, s_notation);
}
 
QString GeoDataCoordinates::latToString(qreal lat, GeoDataCoordinates::Notation notation, GeoDataCoordinates::Unit unit, int precision, char format)
{
    if (notation == GeoDataCoordinates::UTM) {
        /**
         * @FIXME: UTM needs lon + lat to know latitude band and northing
         * By now, this code returns the band+northing of the point
         * (meridian, lat), but this can differ a lot at different locations
         * See bug 347536 https://bugs.kde.org/show_bug.cgi?id=347536
         */
 
        qreal latRad = (unit == Radian) ? lat : lat * DEG2RAD;
 
        QString result = GeoDataCoordinatesPrivate::lonLatToLatitudeBand(0, latRad);
 
        if (precision > 0) {
            QString northingString = QString::number(GeoDataCoordinatesPrivate::lonLatToNorthing(0, latRad), 'f', 2);
            result += QStringLiteral(" %1 m N").arg(northingString);
        }
 
        return result;
    }
 
    QString pmString;
    QString nsString;
 
    if (notation == GeoDataCoordinates::Astro) {
        pmString = (lat > 0) ? QStringLiteral("+") : QStringLiteral("-");
    } else {
        nsString = (lat > 0) ? tr("N") : tr("S");
    }
 
    QString latString;
 
    qreal latDegF = (unit == Degree) ? fabs(lat) : fabs((qreal)(lat)*RAD2DEG);
 
    // Take care of -1 case
    precision = (precision < 0) ? 5 : precision;
 
    if (notation == DMS || notation == DM || notation == Astro) {
        int latDeg = (int)latDegF;
        qreal latMinF = 60 * (latDegF - latDeg);
        int latMin = (int)latMinF;
        qreal latSecF = 60 * (latMinF - latMin);
        int latSec = (int)latSecF;
 
        // Adjustment for fuzziness (like 49.999999999999999999999)
        if (precision == 0) {
            latDeg = qRound(latDegF);
        } else if (precision <= 2) {
            latMin = qRound(latMinF);
        } else if (precision <= 4 && notation == DMS) {
            latSec = qRound(latSecF);
        } else {
            if (notation == DMS || notation == Astro) {
                latSec = latSecF = qRound(latSecF * qPow(10, precision - 4)) / qPow(10, precision - 4);
            } else {
                latMin = latMinF = qRound(latMinF * qPow(10, precision - 2)) / qPow(10, precision - 2);
            }
        }
 
        if (latSec > 59 && (notation == DMS || notation == Astro)) {
            latSecF = 0;
            latSec = latSecF;
            latMin = latMin + 1;
        }
        if (latMin > 59) {
            latMinF = 0;
            latMin = latMinF;
            latDeg = latDeg + 1;
        }
 
        // Evaluate the string
        latString = QString::fromUtf8("%1\xc2\xb0").arg(latDeg, 3, 10, QLatin1Char(' '));
 
        if (precision == 0) {
            return pmString + latString + nsString;
        }
 
        if (notation == DMS || notation == Astro || precision < 3) {
            latString += QStringLiteral(" %2\'").arg(latMin, 2, 10, QLatin1Char('0'));
        }
 
        if (precision < 3) {
            return pmString + latString + nsString;
        }
 
        if (notation == DMS || notation == Astro) {
            // Includes -1 case!
            if (precision < 5) {
                latString += QStringLiteral(" %3\"").arg((double)latSec, 2, 'f', 0, QLatin1Char('0'));
                return latString + nsString;
            }
 
            latString += QStringLiteral(" %L3\"").arg(latSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
        } else {
            latString += QStringLiteral(" %L3'").arg(latMinF, precision + 1, 'f', precision - 2, QLatin1Char('0'));
        }
    } else // notation = GeoDataCoordinates::Decimal
    {
        latString = QString::fromUtf8("%L1\xc2\xb0").arg(latDegF, 4 + precision, format, precision, QLatin1Char(' '));
    }
    return pmString + latString + nsString;
}
 
QString GeoDataCoordinates::latToString() const
{
    return GeoDataCoordinates::latToString(d->m_lat, s_notation);
}
 
bool GeoDataCoordinates::operator==(const GeoDataCoordinates &rhs) const
{
    return *d == *rhs.d;
}
 
bool GeoDataCoordinates::operator!=(const GeoDataCoordinates &rhs) const
{
    return *d != *rhs.d;
}
 
void GeoDataCoordinates::setAltitude(const qreal altitude)
{
    detach();
    d->m_altitude = altitude;
}
 
qreal GeoDataCoordinates::altitude() const
{
    return d->m_altitude;
}
 
int GeoDataCoordinates::utmZone() const
{
    return GeoDataCoordinatesPrivate::lonLatToZone(d->m_lon, d->m_lat);
}
 
qreal GeoDataCoordinates::utmEasting() const
{
    return GeoDataCoordinatesPrivate::lonLatToEasting(d->m_lon, d->m_lat);
}
 
QString GeoDataCoordinates::utmLatitudeBand() const
{
    return GeoDataCoordinatesPrivate::lonLatToLatitudeBand(d->m_lon, d->m_lat);
}
 
qreal GeoDataCoordinates::utmNorthing() const
{
    return GeoDataCoordinatesPrivate::lonLatToNorthing(d->m_lon, d->m_lat);
}
 
quint8 GeoDataCoordinates::detail() const
{
    return d->m_detail;
}
 
void GeoDataCoordinates::setDetail(quint8 detail)
{
    detach();
    d->m_detail = detail;
}
 
GeoDataCoordinates GeoDataCoordinates::rotateAround(const GeoDataCoordinates &axis, qreal angle, Unit unit) const
{
    const Quaternion quatAxis = Quaternion::fromEuler(-axis.latitude(), axis.longitude(), 0);
    const Quaternion rotationAmount = Quaternion::fromEuler(0, 0, unit == Radian ? angle : angle * DEG2RAD);
    const Quaternion resultAxis = quatAxis * rotationAmount * quatAxis.inverse();
 
    return rotateAround(resultAxis);
}
 
GeoDataCoordinates GeoDataCoordinates::rotateAround(const Quaternion &rotAxis) const
{
    Quaternion rotatedQuat = quaternion();
    rotatedQuat.rotateAroundAxis(rotAxis);
    qreal rotatedLon, rotatedLat;
    rotatedQuat.getSpherical(rotatedLon, rotatedLat);
    return {rotatedLon, rotatedLat, altitude()};
}
 
qreal GeoDataCoordinates::bearing(const GeoDataCoordinates &other, Unit unit, BearingType type) const
{
    if (type == FinalBearing) {
        double const offset = unit == Degree ? 180.0 : M_PI;
        return offset + other.bearing(*this, unit, InitialBearing);
    }
 
    qreal const delta = other.d->m_lon - d->m_lon;
    double const bearing = atan2(sin(delta) * cos(other.d->m_lat), cos(d->m_lat) * sin(other.d->m_lat) - sin(d->m_lat) * cos(other.d->m_lat) * cos(delta));
    return unit == Radian ? bearing : bearing * RAD2DEG;
}
 
GeoDataCoordinates GeoDataCoordinates::moveByBearing(qreal bearing, qreal distance) const
{
    qreal newLat = asin(sin(d->m_lat) * cos(distance) + cos(d->m_lat) * sin(distance) * cos(bearing));
    qreal newLon = d->m_lon + atan2(sin(bearing) * sin(distance) * cos(d->m_lat), cos(distance) - sin(d->m_lat) * sin(newLat));
 
    return {newLon, newLat};
}
 
const Quaternion &GeoDataCoordinates::quaternion() const
{
    if (d->m_q == nullptr) {
        d->m_q = new Quaternion(Quaternion::fromSpherical(d->m_lon, d->m_lat));
    }
    return *d->m_q;
}
 
GeoDataCoordinates GeoDataCoordinates::interpolate(const GeoDataCoordinates &target, double t_) const
{
    double const t = qBound(0.0, t_, 1.0);
    Quaternion const quat = Quaternion::slerp(quaternion(), target.quaternion(), t);
    qreal lon, lat;
    quat.getSpherical(lon, lat);
    double const alt = (1.0 - t) * d->m_altitude + t * target.d->m_altitude;
    return {lon, lat, alt};
}
 
GeoDataCoordinates GeoDataCoordinates::nlerp(const GeoDataCoordinates &target, double t) const
{
    qreal lon = 0.0;
    qreal lat = 0.0;
 
    const Quaternion itpos = Quaternion::nlerp(quaternion(), target.quaternion(), t);
    itpos.getSpherical(lon, lat);
 
    const qreal altitude = 0.5 * (d->m_altitude + target.altitude());
 
    return {lon, lat, altitude};
}
 
GeoDataCoordinates
GeoDataCoordinates::interpolate(const GeoDataCoordinates &before, const GeoDataCoordinates &target, const GeoDataCoordinates &after, double t_) const
{
    double const t = qBound(0.0, t_, 1.0);
    Quaternion const b1 = GeoDataCoordinatesPrivate::basePoint(before.quaternion(), quaternion(), target.quaternion());
    Quaternion const a2 = GeoDataCoordinatesPrivate::basePoint(quaternion(), target.quaternion(), after.quaternion());
    Quaternion const a = Quaternion::slerp(quaternion(), target.quaternion(), t);
    Quaternion const b = Quaternion::slerp(b1, a2, t);
    Quaternion c = Quaternion::slerp(a, b, 2 * t * (1.0 - t));
    qreal lon, lat;
    c.getSpherical(lon, lat);
    // @todo spline interpolation of altitude?
    double const alt = (1.0 - t) * d->m_altitude + t * target.d->m_altitude;
    return {lon, lat, alt};
}
 
bool GeoDataCoordinates::isPole(Pole pole) const
{
    // Evaluate the most likely case first:
    // The case where we haven't hit the pole and where our latitude is normalized
    // to the range of 90 deg S ... 90 deg N
    if (fabs((qreal)2.0 * d->m_lat) < M_PI) {
        return false;
    } else {
        if (fabs((qreal)2.0 * d->m_lat) == M_PI) {
            // Ok, we have hit a pole. Now let's check whether it's the one we've asked for:
            if (pole == AnyPole) {
                return true;
            } else {
                if (pole == NorthPole && 2.0 * d->m_lat == +M_PI) {
                    return true;
                }
                if (pole == SouthPole && 2.0 * d->m_lat == -M_PI) {
                    return true;
                }
                return false;
            }
        }
        //
        else {
            // FIXME: Should we just normalize latitude and longitude and be done?
            //        While this might work well for persistent data it would create some
            //        possible overhead for temporary data, so this needs careful thinking.
            mDebug() << "GeoDataCoordinates not normalized!";
 
            // Only as a last resort we cover the unlikely case where
            // the latitude is not normalized to the range of
            // 90 deg S ... 90 deg N
            if (fabs((qreal)2.0 * normalizeLat(d->m_lat)) < M_PI) {
                return false;
            } else {
                // Ok, we have hit a pole. Now let's check whether it's the one we've asked for:
                if (pole == AnyPole) {
                    return true;
                } else {
                    if (pole == NorthPole && 2.0 * d->m_lat == +M_PI) {
                        return true;
                    }
                    if (pole == SouthPole && 2.0 * d->m_lat == -M_PI) {
                        return true;
                    }
                    return false;
                }
            }
        }
    }
}
 
qreal GeoDataCoordinates::sphericalDistanceTo(const GeoDataCoordinates &other) const
{
    qreal lon2, lat2;
    other.geoCoordinates(lon2, lat2);
 
    // FIXME: Take the altitude into account!
 
    return distanceSphere(d->m_lon, d->m_lat, lon2, lat2);
}
 
GeoDataCoordinates &GeoDataCoordinates::operator=(const GeoDataCoordinates &other)
{
    qAtomicAssign(d, other.d);
    return *this;
}
 
void GeoDataCoordinates::pack(QDataStream &stream) const
{
    stream << d->m_lon;
    stream << d->m_lat;
    stream << d->m_altitude;
}
 
void GeoDataCoordinates::unpack(QDataStream &stream)
{
    // call detach even though it shouldn't be needed - one never knows
    detach();
    stream >> d->m_lon;
    stream >> d->m_lat;
    stream >> d->m_altitude;
}
 
Quaternion GeoDataCoordinatesPrivate::basePoint(const Quaternion &q1, const Quaternion &q2, const Quaternion &q3)
{
    Quaternion const a = (q2.inverse() * q3).log();
    Quaternion const b = (q2.inverse() * q1).log();
    return q2 * ((a + b) * -0.25).exp();
}
 
qreal GeoDataCoordinatesPrivate::arcLengthOfMeridian(qreal phi)
{
    // Precalculate n
    qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis)
        / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis);
 
    // Precalculate alpha
    qreal const alpha = ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0)
        * (1.0 + (qPow(n, 2.0) / 4.0) + (qPow(n, 4.0) / 64.0));
 
    // Precalculate beta
    qreal const beta = (-3.0 * n / 2.0) + (9.0 * qPow(n, 3.0) / 16.0) + (-3.0 * qPow(n, 5.0) / 32.0);
 
    // Precalculate gamma
    qreal const gamma = (15.0 * qPow(n, 2.0) / 16.0) + (-15.0 * qPow(n, 4.0) / 32.0);
 
    // Precalculate delta
    qreal const delta = (-35.0 * qPow(n, 3.0) / 48.0) + (105.0 * qPow(n, 5.0) / 256.0);
 
    // Precalculate epsilon
    qreal const epsilon = (315.0 * qPow(n, 4.0) / 512.0);
 
    // Now calculate the sum of the series and return
    qreal const result = alpha * (phi + (beta * qSin(2.0 * phi)) + (gamma * qSin(4.0 * phi)) + (delta * qSin(6.0 * phi)) + (epsilon * qSin(8.0 * phi)));
 
    return result;
}
 
qreal GeoDataCoordinatesPrivate::centralMeridianUTM(qreal zone)
{
    return DEG2RAD * (-183.0 + (zone * 6.0));
}
 
qreal GeoDataCoordinatesPrivate::footpointLatitude(qreal northing)
{
    // Precalculate n (Eq. 10.18)
    qreal const n = (GeoDataCoordinatesPrivate::sm_semiMajorAxis - GeoDataCoordinatesPrivate::sm_semiMinorAxis)
        / (GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis);
 
    // Precalculate alpha (Eq. 10.22)
    // (Same as alpha in Eq. 10.17)
    qreal const alpha =
        ((GeoDataCoordinatesPrivate::sm_semiMajorAxis + GeoDataCoordinatesPrivate::sm_semiMinorAxis) / 2.0) * (1 + (qPow(n, 2.0) / 4) + (qPow(n, 4.0) / 64));
 
    // Precalculate y (Eq. 10.23)
    qreal const y = northing / alpha;
 
    // Precalculate beta (Eq. 10.22)
    qreal const beta = (3.0 * n / 2.0) + (-27.0 * qPow(n, 3.0) / 32.0) + (269.0 * qPow(n, 5.0) / 512.0);
 
    // Precalculate gamma (Eq. 10.22)
    qreal const gamma = (21.0 * qPow(n, 2.0) / 16.0) + (-55.0 * qPow(n, 4.0) / 32.0);
 
    // Precalculate delta (Eq. 10.22)
    qreal const delta = (151.0 * qPow(n, 3.0) / 96.0) + (-417.0 * qPow(n, 5.0) / 128.0);
 
    // Precalculate epsilon (Eq. 10.22)
    qreal const epsilon = (1097.0 * qPow(n, 4.0) / 512.0);
 
    // Now calculate the sum of the series (Eq. 10.21)
    qreal const result = y + (beta * qSin(2.0 * y)) + (gamma * qSin(4.0 * y)) + (delta * qSin(6.0 * y)) + (epsilon * qSin(8.0 * y));
 
    return result;
}
 
QPointF GeoDataCoordinatesPrivate::mapLonLatToXY(qreal lambda, qreal phi, qreal lambda0)
{
    // Equation (10.15)
 
    // Precalculate second numerical eccentricity
    const qreal ep2 = (qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) - qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0))
        / qPow(GeoDataCoordinatesPrivate::sm_semiMinorAxis, 2.0);
 
    // Precalculate the square of nu, just an auxilar quantity
    const qreal nu2 = ep2 * qPow(qCos(phi), 2.0);
 
    // Precalculate the radius of curvature in prime vertical
    const qreal N = qPow(GeoDataCoordinatesPrivate::sm_semiMajorAxis, 2.0) / (GeoDataCoordinatesPrivate::sm_semiMinorAxis * qSqrt(1 + nu2));
 
    // Precalculate the tangent of phi and its square
    const qreal t = qTan(phi);
    const qreal t2 = t * t;
 
    // Precalculate longitude difference
    const qreal l = lambda - lambda0;
 
    /*
     * Precalculate coefficients for l**n in the equations below
     * so a normal human being can read the expressions for easting
     * and northing
     * -- l**1 and l**2 have coefficients of 1.0
     *
     * The actual used coefficients starts at coef[1], just to
     * follow the meaningful nomenclature in equation 10.15
     * (coef[n] corresponds to qPow(l,n) factor)
     */
 
    const qreal coef1 = 1;
    const qreal coef2 = 1;
    const qreal coef3 = 1.0 - t2 + nu2;
    const qreal coef4 = 5.0 - t2 + 9 * nu2 + 4.0 * (nu2 * nu2);
    const qreal coef5 = 5.0 - 18.0 * t2 + (t2 * t2) + 14.0 * nu2 - 58.0 * t2 * nu2;
    const qreal coef6 = 61.0 - 58.0 * t2 + (t2 * t2) + 270.0 * nu2 - 330.0 * t2 * nu2;
    const qreal coef7 = 61.0 - 479.0 * t2 + 179.0 * (t2 * t2) - (t2 * t2 * t2);
    const qreal coef8 = 1385.0 - 3111.0 * t2 + 543.0 * (t2 * t2) - (t2 * t2 * t2);
 
    // Calculate easting (x)
    const qreal easting = N * qCos(phi) * coef1 * l + (N / 6.0 * qPow(qCos(phi), 3.0) * coef3 * qPow(l, 3.0))
        + (N / 120.0 * qPow(qCos(phi), 5.0) * coef5 * qPow(l, 5.0)) + (N / 5040.0 * qPow(qCos(phi), 7.0) * coef7 * qPow(l, 7.0));
 
    // Calculate northing (y)
    const qreal northing = arcLengthOfMeridian(phi) + (t / 2.0 * N * qPow(qCos(phi), 2.0) * coef2 * qPow(l, 2.0))
        + (t / 24.0 * N * qPow(qCos(phi), 4.0) * coef4 * qPow(l, 4.0)) + (t / 720.0 * N * qPow(qCos(phi), 6.0) * coef6 * qPow(l, 6.0))
        + (t / 40320.0 * N * qPow(qCos(phi), 8.0) * coef8 * qPow(l, 8.0));
 
    return {easting, northing};
}
 
int GeoDataCoordinatesPrivate::lonLatToZone(qreal lon, qreal lat)
{
    // Converts lon and lat to degrees
    qreal lonDeg = lon * RAD2DEG;
    qreal latDeg = lat * RAD2DEG;
 
    /* Round the value of the longitude when the distance to the nearest integer
     * is less than 0.0000001. This avoids fuzzy values such as -114.0000000001, which
     * can produce a misbehaviour when calculating the zone associated at the borders
     * of the zone intervals (for example, the interval [-114, -108[ is associated with
     * zone number 12; if the following rounding is not done, the value returned by
     * lonLatToZone(114,0) is 11 instead of 12, as the function actually receives
     * -114.0000000001, which is in the interval [-120,-114[, associated to zone 11
     */
    qreal precision = 0.0000001;
 
    if (qAbs(lonDeg - qFloor(lonDeg)) < precision || qAbs(lonDeg - qCeil(lonDeg)) < precision) {
        lonDeg = qRound(lonDeg);
    }
 
    // There is no numbering associated to the poles, special value 0 is returned.
    if (latDeg < -80 || latDeg > 84) {
        return 0;
    }
 
    // Obtains the zone number handling all the so called "exceptions"
    // See problem: https://en.wikipedia.org/wiki/Universal_Transverse_Mercator_coordinate_system#Exceptions
    // See solution: https://gis.stackexchange.com/questions/13291/computing-utm-zone-from-lat-long-point
 
    // General
    int zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1;
 
    // Southwest Norway
    if (latDeg >= 56 && latDeg < 64 && lonDeg >= 3 && lonDeg < 12) {
        zoneNumber = 32;
    }
 
    // Svalbard
    if (latDeg >= 72 && latDeg < 84) {
        if (lonDeg >= 0 && lonDeg < 9) {
            zoneNumber = 31;
        } else if (lonDeg >= 9 && lonDeg < 21) {
            zoneNumber = 33;
        } else if (lonDeg >= 21 && lonDeg < 33) {
            zoneNumber = 35;
        } else if (lonDeg >= 33 && lonDeg < 42) {
            zoneNumber = 37;
        }
    }
 
    return zoneNumber;
}
 
qreal GeoDataCoordinatesPrivate::lonLatToEasting(qreal lon, qreal lat)
{
    int zoneNumber = lonLatToZone(lon, lat);
 
    if (zoneNumber == 0) {
        qreal lonDeg = lon * RAD2DEG;
        zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1;
    }
 
    QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY(lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber));
 
    // Adjust easting and northing for UTM system
    qreal easting = coordinates.x() * GeoDataCoordinatesPrivate::sm_utmScaleFactor + 500000.0;
 
    return easting;
}
 
QString GeoDataCoordinatesPrivate::lonLatToLatitudeBand(qreal lon, qreal lat)
{
    // Obtains the latitude bands handling all the so called "exceptions"
 
    // Converts lon and lat to degrees
    qreal lonDeg = lon * RAD2DEG;
    qreal latDeg = lat * RAD2DEG;
 
    // Regular latitude bands between 80 S and 80 N (that is, between 10 and 170 in the [0,180] interval)
    int bandLetterIndex = 24; // Avoids "may be used uninitialized" warning
 
    if (latDeg < -80) {
        // South pole (A for zones 1-30, B for zones 31-60)
        bandLetterIndex = ((lonDeg + 180) < 6 * 31) ? 0 : 1;
    } else if (latDeg >= -80 && latDeg <= 80) {
        // General (+2 because the general lettering starts in C)
        bandLetterIndex = static_cast<int>((latDeg + 80.0) / 8.0) + 2;
    } else if (latDeg >= 80 && latDeg < 84) {
        // Band X is extended 4 more degrees
        bandLetterIndex = 21;
    } else if (latDeg >= 84) {
        // North pole (Y for zones 1-30, Z for zones 31-60)
        bandLetterIndex = ((lonDeg + 180) < 6 * 31) ? 22 : 23;
    }
 
    return QStringLiteral("ABCDEFGHJKLMNPQRSTUVWXYZ?").at(bandLetterIndex);
}
 
qreal GeoDataCoordinatesPrivate::lonLatToNorthing(qreal lon, qreal lat)
{
    int zoneNumber = lonLatToZone(lon, lat);
 
    if (zoneNumber == 0) {
        qreal lonDeg = lon * RAD2DEG;
        zoneNumber = static_cast<int>((lonDeg + 180) / 6.0) + 1;
    }
 
    QPointF coordinates = GeoDataCoordinatesPrivate::mapLonLatToXY(lon, lat, GeoDataCoordinatesPrivate::centralMeridianUTM(zoneNumber));
 
    qreal northing = coordinates.y() * GeoDataCoordinatesPrivate::sm_utmScaleFactor;
 
    if (northing < 0.0) {
        northing += 10000000.0;
    }
 
    return northing;
}
 
size_t qHash(const GeoDataCoordinates &coordinates)
{
    uint seed = ::qHash(coordinates.altitude());
    seed = ::qHash(coordinates.latitude(), seed);
 
    return ::qHash(coordinates.longitude(), seed);
}
 
}