GeoDataCoordinates.cpp

// SPDX-License-Identifier: LGPL-2.1-or-later
//
// SPDX-FileCopyrightText: 2004-2007 Torsten Rahn <[email protected]>
// SPDX-FileCopyrightText: 2007-2008 Inge Wallin <[email protected]>
// SPDX-FileCopyrightText: 2008 Patrick Spendrin <[email protected]>
// SPDX-FileCopyrightText: 2011 Friedrich W. H. Kossebau <[email protected]>
// SPDX-FileCopyrightText: 2011 Bernhard Beschow <[email protected]>
// SPDX-FileCopyrightText: 2015 Alejandro Garcia Montoro <[email protected]>
//
 
 
#include "GeoDataCoordinates.h"
#include "GeoDataCoordinates_p.h"
#include "LonLatParser_p.h"
 
#include <qmath.h>
#include <QDataStream>
#include <QPointF>
 
#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.load() == 1) {
        return;
    }
 
    GeoDataCoordinatesPrivate *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 GeoDataCoordinates( parser.lon(), parser.lat(), 0, GeoDataCoordinates::Degree );
    } else {
        return GeoDataCoordinates();
    }
}
 
 
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 QString("%1%2 %3 m E, %4 m N").arg(zoneString, bandString, eastingString, northingString);
        }
        else{
            coordString = lonToString( d->m_lon, notation, Radian, precision )
                        + QLatin1String(", ")
                        + 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 += QString(" %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 += QString(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0'));
        }
 
        if ( precision < 3 ) {
            return lonString + weString;
        }
 
        if ( notation == DMS ) {
            // Includes -1 case!
            if ( precision < 5 ) {
                lonString += QString(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0'));
                return lonString + weString;
            }
 
            lonString += QString(" %L3\"").arg(lonSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
        }
        else {
            lonString += QString(" %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 += QString(" %2\'").arg(lonMin, 2, 10, QLatin1Char('0'));
 
        if ( precision < 3 ) {
            return lonString;
        }
 
        // Includes -1 case!
        if ( precision < 5 ) {
            lonString += QString(" %3\"").arg(lonSec, 2, 'f', 0, QLatin1Char('0'));
            return lonString;
        }
 
        lonString += QString(" %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 += QString(" %1 m N").arg(northingString);
        }
 
        return result;
    }
 
    QString pmString;
    QString nsString;
 
    if (notation == GeoDataCoordinates::Astro){
        pmString = ( lat > 0 ) ? "+" : "-";
    }
    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 += QString(" %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 += QString(" %3\"").arg(latSec, 2, 'f', 0, QLatin1Char('0'));
                return latString + nsString;
            }
 
            latString += QString(" %L3\"").arg(latSecF, precision - 1, 'f', precision - 4, QLatin1Char('0'));
        }
        else {
            latString += QString(" %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 GeoDataCoordinates(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 GeoDataCoordinates( 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 GeoDataCoordinates( 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 GeoDataCoordinates(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 GeoDataCoordinates( 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 QPointF(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 QString( "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;
}
 
uint qHash(const GeoDataCoordinates &coordinates)
{
    uint seed = ::qHash(coordinates.altitude());
    seed = ::qHash(coordinates.latitude(), seed);
 
    return ::qHash(coordinates.longitude(), seed);
}
 
}