Logo Search packages:      
Sourcecode: qlandkarte version File versions  Download package

IProjection.cpp

/**********************************************************************************************
    Copyright (C) 2006, 2007 Oliver Eichler oliver.eichler@gmx.de

    This program is free software; you can redistribute it and/or modify
    it under the terms of the GNU General Public License as published by
    the Free Software Foundation; either version 2 of the License, or
    (at your option) any later version.

    This program is distributed in the hope that it will be useful,
    but WITHOUT ANY WARRANTY; without even the implied warranty of
    MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
    GNU General Public License for more details.

    You should have received a copy of the GNU General Public License
    along with this program; if not, write to the Free Software
    Foundation, Inc., 59 Temple Place - Suite 330, Boston, MA 02111 USA

**********************************************************************************************/

#include "IProjection.h"
#include "CProjMerc.h"
#include "CProjMercDynamic.h"
#include "CProjEckII.h"
#include "CProjUTM.h"

#include <QtCore>

static CProjMerc projMerc;
static CProjMercDynamic projMercDyn;
static CProjEckII projEckII;
static CProjUTM projUTM;

const proj_t projections[] =
{
    {0,QObject::tr("Mercator"),&projMerc}
    ,{1,QObject::tr("Mercator (dynamic)"),&projMercDyn}
    ,{2,QObject::tr("Eckert II"),&projEckII}
    ,{3,QObject::tr("UTM"),&projUTM}
    ,{0,0,0}
};

IProjection * gpProj = &projMercDyn;

void setupProjection(quint32 key)
{
    if(key < (sizeof(projections) / sizeof(proj_t) - 1)) {
        gpProj = projections[key].proj;
    }
}


IProjection::IProjection()
{

}


IProjection::~IProjection()
{

}


// from http://www.movable-type.co.uk/scripts/LatLongVincenty.html
double IProjection::distance(const XY& p1, const XY& p2, double& a1, double& a2)
{
    double cosSigma = 0.0;
    double sigma = 0.0;
    double sinAlpha = 0.0;
    double cosSqAlpha = 0.0;
    double cos2SigmaM = 0.0;
    double sinSigma = 0.0;
    double sinLambda = 0.0;
    double cosLambda = 0.0;

                                 // WGS-84 ellipsiod
    double a = 6378137.0, b = 6356752.3142,  f = 1.0/298.257223563;
    double L = p2.u - p1.u;

    double U1 = atan((1-f) * tan(p1.v));
    double U2 = atan((1-f) * tan(p2.v));
    double sinU1 = sin(U1), cosU1 = cos(U1);
    double sinU2 = sin(U2), cosU2 = cos(U2);
    double lambda = L, lambdaP = 2*PI;
    unsigned iterLimit = 20;

    while ((fabs(lambda - lambdaP) > 1e-12) && (--iterLimit > 0)) {
        sinLambda = sin(lambda);
        cosLambda = cos(lambda);
        sinSigma = sqrt((cosU2*sinLambda) * (cosU2*sinLambda) + (cosU1*sinU2-sinU1*cosU2*cosLambda) * (cosU1*sinU2-sinU1*cosU2*cosLambda));

        if (sinSigma==0) {
            return 0;            // co-incident points
        }

        cosSigma = sinU1 * sinU2 + cosU1 * cosU2 * cosLambda;
        sigma = atan2(sinSigma, cosSigma);
        sinAlpha = cosU1 * cosU2 * sinLambda / sinSigma;
        cosSqAlpha = 1 - sinAlpha * sinAlpha;
        cos2SigmaM = cosSigma - 2 * sinU1 * sinU2 / cosSqAlpha;

        if (isnan(cos2SigmaM)) {
            cos2SigmaM = 0;      // equatorial line: cosSqAlpha=0 (6)
        }

        double C = f/16 * cosSqAlpha * (4 + f * (4 - 3 * cosSqAlpha));
        lambdaP = lambda;
        lambda = L + (1-C) * f * sinAlpha * (sigma + C*sinSigma*(cos2SigmaM + C * cosSigma * (-1 + 2 * cos2SigmaM * cos2SigmaM)));
    }
                                 // formula failed to converge
    if (iterLimit==0) return FP_NAN;

    double uSq = cosSqAlpha * (a*a - b*b) / (b*b);
    double A = 1 + uSq/16384*(4096+uSq*(-768+uSq*(320-175*uSq)));
    double B = uSq/1024 * (256+uSq*(-128+uSq*(74-47*uSq)));
    double deltaSigma = B*sinSigma*(cos2SigmaM+B/4*(cosSigma*(-1+2*cos2SigmaM*cos2SigmaM)-B/6*cos2SigmaM*(-3+4*sinSigma*sinSigma)*(-3+4*cos2SigmaM*cos2SigmaM)));
    double s = b*A*(sigma-deltaSigma);

    a1 = atan2(cosU2 * sinLambda, cosU1 * sinU2 - sinU1 * cosU2 * cosLambda) * 360 / TWOPI;
    a2 = atan2(cosU1 * sinLambda, -sinU1 * cosU2 + cosU1 * sinU2 * cosLambda) * 360 / TWOPI;
    return s;
}

Generated by  Doxygen 1.6.0   Back to index