root/GPXSee
1
0
Fork 0
mirror of https://github.com/tumic0/GPXSee.git synced 2024-11-28 13:41:16 +01:00
Code Issues Projects Releases Activity

Datum shift optimization (precompute ellipsoid constants)

Browse Source 
+ projections code unification
This commit is contained in:
Martin Tůma 2018-05-17 22:41:56 +02:00
parent b3940283a8
commit a486f9e78d
15 changed files with 98 additions and 105 deletions
Show Stats Download Patch File Download Diff File Expand all files Collapse all files

67
src/map/albersequal.cpp
View File

@ -46,11 +46,11 @@ Defense.
#define ONE_MINUS_SQR(x) (1.0 - (x) * (x)) #define ONE_MINUS_SQR(x) (1.0 - (x) * (x))
#define ALBERS_Q(slat, one_minus_sqr_es_sin, es_sin) \ #define ALBERS_Q(slat, one_minus_sqr_e_sin, es_sin) \
(_one_minus_es2 * ((slat) / (one_minus_sqr_es_sin) - \ (_one_minus_es * ((slat) / (one_minus_sqr_e_sin) - \
(1 / (_two_es)) * log((1 - (es_sin)) / (1 + (es_sin))))) (1 / (_two_e)) * log((1 - (es_sin)) / (1 + (es_sin)))))
#define ALBERS_M(clat, one_minus_sqr_es_sin) \ #define ALBERS_M(clat, one_minus_sqr_e_sin) \
((clat) / sqrt(one_minus_sqr_es_sin)) ((clat) / sqrt(one_minus_sqr_e_sin))
AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1, AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
@ -60,8 +60,8 @@ AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
double sin_lat, sin_lat1, sin_lat2, cos_lat1, cos_lat2; double sin_lat, sin_lat1, sin_lat2, cos_lat1, cos_lat2;
double m1, m2, sqr_m1, sqr_m2; double m1, m2, sqr_m1, sqr_m2;
double q0, q1, q2; double q0, q1, q2;
double es_sin, es_sin1, es_sin2; double e_sin, e_sin1, e_sin2;
double one_minus_sqr_es_sin1, one_minus_sqr_es_sin2; double one_minus_sqr_e_sin1, one_minus_sqr_e_sin2;
double nq0; double nq0;
double sp1, sp2; double sp1, sp2;
@ -75,31 +75,30 @@ AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
sp2 = deg2rad(standardParallel2); sp2 = deg2rad(standardParallel2);
_a2 = ellipsoid->radius() * ellipsoid->radius(); _a2 = ellipsoid->radius() * ellipsoid->radius();
_es2 = 2 * ellipsoid->flattening() - ellipsoid->flattening() _es = ellipsoid->es();
* ellipsoid->flattening(); _e = sqrt(_es);
_es = sqrt(_es2); _one_minus_es = 1 - _es;
_one_minus_es2 = 1 - _es2; _two_e = 2 * _e;
_two_es = 2 * _es;
sin_lat = sin(_latitudeOrigin); sin_lat = sin(_latitudeOrigin);
es_sin = _es * sin_lat; e_sin = _e * sin_lat;
q0 = ALBERS_Q(sin_lat, ONE_MINUS_SQR(es_sin), es_sin); q0 = ALBERS_Q(sin_lat, ONE_MINUS_SQR(e_sin), e_sin);
sin_lat1 = sin(sp1); sin_lat1 = sin(sp1);
cos_lat1 = cos(sp1); cos_lat1 = cos(sp1);
es_sin1 = _es * sin_lat1; e_sin1 = _e * sin_lat1;
one_minus_sqr_es_sin1 = ONE_MINUS_SQR(es_sin1); one_minus_sqr_e_sin1 = ONE_MINUS_SQR(e_sin1);
m1 = ALBERS_M(cos_lat1, one_minus_sqr_es_sin1); m1 = ALBERS_M(cos_lat1, one_minus_sqr_e_sin1);
q1 = ALBERS_Q(sin_lat1, one_minus_sqr_es_sin1, es_sin1); q1 = ALBERS_Q(sin_lat1, one_minus_sqr_e_sin1, e_sin1);
sqr_m1 = m1 * m1; sqr_m1 = m1 * m1;
if (fabs(sp1 - sp2) > 1.0e-10) { if (fabs(sp1 - sp2) > 1.0e-10) {
sin_lat2 = sin(sp2); sin_lat2 = sin(sp2);
cos_lat2 = cos(sp2); cos_lat2 = cos(sp2);
es_sin2 = _es * sin_lat2; e_sin2 = _e * sin_lat2;
one_minus_sqr_es_sin2 = ONE_MINUS_SQR(es_sin2); one_minus_sqr_e_sin2 = ONE_MINUS_SQR(e_sin2);
m2 = ALBERS_M(cos_lat2, one_minus_sqr_es_sin2); m2 = ALBERS_M(cos_lat2, one_minus_sqr_e_sin2);
q2 = ALBERS_Q(sin_lat2, one_minus_sqr_es_sin2, es_sin2); q2 = ALBERS_Q(sin_lat2, one_minus_sqr_e_sin2, e_sin2);
sqr_m2 = m2 * m2; sqr_m2 = m2 * m2;
_n = (sqr_m1 - sqr_m2) / (q2 - q1); _n = (sqr_m1 - sqr_m2) / (q2 - q1);
} else } else
@ -115,7 +114,7 @@ PointD AlbersEqual::ll2xy(const Coordinates &c) const
{ {
double dlam; double dlam;
double sin_lat; double sin_lat;
double es_sin; double e_sin;
double q; double q;
double rho; double rho;
double theta; double theta;
@ -129,8 +128,8 @@ PointD AlbersEqual::ll2xy(const Coordinates &c) const
dlam += M_2_PI; dlam += M_2_PI;
sin_lat = sin(deg2rad(c.lat())); sin_lat = sin(deg2rad(c.lat()));
es_sin = _es * sin_lat; e_sin = _e * sin_lat;
q = ALBERS_Q(sin_lat, ONE_MINUS_SQR(es_sin), es_sin); q = ALBERS_Q(sin_lat, ONE_MINUS_SQR(e_sin), e_sin);
nq = _n * q; nq = _n * q;
rho = (_C < nq) ? 0 : _a_over_n * sqrt(_C - nq); rho = (_C < nq) ? 0 : _a_over_n * sqrt(_C - nq);
theta = _n * dlam; theta = _n * dlam;
@ -147,7 +146,7 @@ Coordinates AlbersEqual::xy2ll(const PointD &p) const
double rho, rho_n; double rho, rho_n;
double phi, delta_phi = 1.0; double phi, delta_phi = 1.0;
double sin_phi; double sin_phi;
double es_sin, one_minus_sqr_es_sin; double e_sin, one_minus_sqr_e_sin;
double theta = 0.0; double theta = 0.0;
int count = 30; int count = 30;
double tolerance = 4.85e-10; double tolerance = 4.85e-10;
@ -170,7 +169,7 @@ Coordinates AlbersEqual::xy2ll(const PointD &p) const
theta = atan2(dx, rho0_minus_dy); theta = atan2(dx, rho0_minus_dy);
rho_n = rho * _n; rho_n = rho * _n;
q = (_C - (rho_n * rho_n) / _a2) / _n; q = (_C - (rho_n * rho_n) / _a2) / _n;
qc = 1 - ((_one_minus_es2) / (_two_es)) * log((1.0 - _es) / (1.0 + _es)); qc = 1 - ((_one_minus_es) / (_two_e)) * log((1.0 - _e) / (1.0 + _e));
if (fabs(fabs(qc) - fabs(q)) > 1.0e-6) { if (fabs(fabs(qc) - fabs(q)) > 1.0e-6) {
q_over_2 = q / 2.0; q_over_2 = q / 2.0;
if (q_over_2 > 1.0) if (q_over_2 > 1.0)
@ -179,17 +178,17 @@ Coordinates AlbersEqual::xy2ll(const PointD &p) const
lat = -M_PI_2; lat = -M_PI_2;
else { else {
phi = asin(q_over_2); phi = asin(q_over_2);
if (_es < 1.0e-10) if (_e < 1.0e-10)
lat = phi; lat = phi;
else { else {
while ((fabs(delta_phi) > tolerance) && count) { while ((fabs(delta_phi) > tolerance) && count) {
sin_phi = sin(phi); sin_phi = sin(phi);
es_sin = _es * sin_phi; e_sin = _e * sin_phi;
one_minus_sqr_es_sin = ONE_MINUS_SQR(es_sin); one_minus_sqr_e_sin = ONE_MINUS_SQR(e_sin);
delta_phi = (one_minus_sqr_es_sin * one_minus_sqr_es_sin) delta_phi = (one_minus_sqr_e_sin * one_minus_sqr_e_sin)
/ (2.0 * cos(phi)) * (q / (_one_minus_es2) - sin_phi / (2.0 * cos(phi)) * (q / (_one_minus_es) - sin_phi
/ one_minus_sqr_es_sin + (log((1.0 - es_sin) / one_minus_sqr_e_sin + (log((1.0 - e_sin)
/ (1.0 + es_sin)) / (_two_es))); / (1.0 + e_sin)) / (_two_e)));
phi += delta_phi; phi += delta_phi;
count --; count --;
} }

6
src/map/albersequal.h
View File

@ -27,11 +27,11 @@ private:
double _rho0; double _rho0;
double _C; double _C;
double _n; double _n;
double _e;
double _es; double _es;
double _es2;
double _a_over_n; double _a_over_n;
double _one_minus_es2; double _one_minus_es;
double _two_es; double _two_e;
}; };
#endif // ALBERSEQUAL_H #endif // ALBERSEQUAL_H

4
src/map/datum.cpp
View File

@ -113,7 +113,7 @@ Coordinates Datum::toWGS84(const Coordinates &c) const
switch (_transformation) { switch (_transformation) {
case Helmert: case Helmert:
return Geocentric::toGeodetic(helmert(Geocentric::fromGeodetic(c, return Geocentric::toGeodetic(helmert(Geocentric::fromGeodetic(c,
*this)), WGS84); ellipsoid())), WGS84.ellipsoid());
case Molodensky: case Molodensky:
return molodensky(c, *this, WGS84); return molodensky(c, *this, WGS84);
default: default:
@ -126,7 +126,7 @@ Coordinates Datum::fromWGS84(const Coordinates &c) const
switch (_transformation) { switch (_transformation) {
case Helmert: case Helmert:
return Geocentric::toGeodetic(helmertr(Geocentric::fromGeodetic(c, return Geocentric::toGeodetic(helmertr(Geocentric::fromGeodetic(c,
WGS84)), *this); WGS84.ellipsoid())), ellipsoid());
case Molodensky: case Molodensky:
return molodensky(c, WGS84, *this); return molodensky(c, WGS84, *this);
default: default:

7
src/map/ellipsoid.cpp
View File

@ -66,6 +66,13 @@ void Ellipsoid::loadList(const QString &path)
} }
} }
Ellipsoid::Ellipsoid(double radius, double flattening)
: _radius(radius), _flattening(flattening)
{
_es = 2.0 * flattening - flattening * flattening;
_b = radius * (1.0 - flattening);
}
#ifndef QT_NO_DEBUG #ifndef QT_NO_DEBUG
QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid) QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid)
{ {

8
src/map/ellipsoid.h
View File

@ -9,12 +9,13 @@
class Ellipsoid class Ellipsoid
{ {
public: public:
Ellipsoid() : _radius(NAN), _flattening(NAN) {} Ellipsoid() : _radius(NAN), _flattening(NAN), _es(NAN), _b(NAN) {}
Ellipsoid(double radius, double flattening) Ellipsoid(double radius, double flattening);
: _radius(radius), _flattening(flattening) {}
double radius() const {return _radius;} double radius() const {return _radius;}
double flattening() const {return _flattening;} double flattening() const {return _flattening;}
double es() const {return _es;}
double b() const {return _b;}
bool isNull() const bool isNull() const
{return (std::isnan(_radius) && std::isnan(_flattening));} {return (std::isnan(_radius) && std::isnan(_flattening));}
@ -27,6 +28,7 @@ public:
private: private:
double _radius; double _radius;
double _flattening; double _flattening;
double _es, _b;
static QMap<int, Ellipsoid> WGS84(); static QMap<int, Ellipsoid> WGS84();
static QMap<int, Ellipsoid> _ellipsoids; static QMap<int, Ellipsoid> _ellipsoids;

24
src/map/geocentric.cpp
View File

@ -41,42 +41,34 @@ Defense.
*/ */
#include "datum.h" #include "ellipsoid.h"
#include "geocentric.h" #include "geocentric.h"
#define AD_C 1.0026000 #define AD_C 1.0026000
Point3D Geocentric::fromGeodetic(const Coordinates &c, const Datum &datum) Point3D Geocentric::fromGeodetic(const Coordinates &c, const Ellipsoid *e)
{ {
const Ellipsoid *e = datum.ellipsoid();
double e2 = 2.0 * e->flattening() - e->flattening() * e->flattening();
double lat = deg2rad(c.lat()); double lat = deg2rad(c.lat());
double lon = deg2rad(c.lon()); double lon = deg2rad(c.lon());
double slat = sin(lat); double slat = sin(lat);
double slat2 = slat * slat; double slat2 = slat * slat;
double clat = cos(lat); double clat = cos(lat);
double Rn = e->radius() / (sqrt(1.0 - e2 * slat2)); double Rn = e->radius() / (sqrt(1.0 - e->es() * slat2));
if (lon > M_PI) if (lon > M_PI)
lon -= M_2_PI; lon -= M_2_PI;
return Point3D(Rn * clat * cos(lon), Rn * clat * sin(lon), return Point3D(Rn * clat * cos(lon), Rn * clat * sin(lon),
(Rn * (1 - e2)) * slat); (Rn * (1 - e->es())) * slat);
} }
Coordinates Geocentric::toGeodetic(const Point3D &p, const Datum &datum) Coordinates Geocentric::toGeodetic(const Point3D &p, const Ellipsoid *e)
{ {
const Ellipsoid *e = datum.ellipsoid(); double ep2 = (1.0 / (1.0 - e->es())) - 1.0;
double b = e->radius() * (1.0 - e->flattening());
double e2 = 2.0 * e->flattening() - e->flattening() * e->flattening();
double ep2 = (1.0 / (1.0 - e2)) - 1.0;
bool pole = false; bool pole = false;
double lat, lon; double lat, lon;
if (p.x() == 0.0) { if (p.x() == 0.0) {
if (p.y() > 0) if (p.y() > 0)
lon = M_PI_2; lon = M_PI_2;
@ -102,8 +94,8 @@ Coordinates Geocentric::toGeodetic(const Point3D &p, const Datum &datum)
double Sin_B0 = T0 / S0; double Sin_B0 = T0 / S0;
double Cos_B0 = W / S0; double Cos_B0 = W / S0;
double Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0; double Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
double T1 = p.z() + b * ep2 * Sin3_B0; double T1 = p.z() + e->b() * ep2 * Sin3_B0;
double Sum = W - e->radius() * e2 * Cos_B0 * Cos_B0 * Cos_B0; double Sum = W - e->radius() * e->es() * Cos_B0 * Cos_B0 * Cos_B0;
double S1 = sqrt(T1*T1 + Sum * Sum); double S1 = sqrt(T1*T1 + Sum * Sum);
double Sin_p1 = T1 / S1; double Sin_p1 = T1 / S1;
double Cos_p1 = Sum / S1; double Cos_p1 = Sum / S1;

6
src/map/geocentric.h
View File

@ -4,7 +4,7 @@
#include <cmath> #include <cmath>
#include "common/coordinates.h" #include "common/coordinates.h"
class Datum; class Ellipsoid;
class Point3D { class Point3D {
public: public:
@ -22,8 +22,8 @@ private:
}; };
namespace Geocentric { namespace Geocentric {
Point3D fromGeodetic(const Coordinates &c, const Datum &datum); Point3D fromGeodetic(const Coordinates &c, const Ellipsoid *e);
Coordinates toGeodetic(const Point3D &p, const Datum &datum); Coordinates toGeodetic(const Point3D &p, const Ellipsoid *e);
} }
#endif // GEOCENTRIC_H #endif // GEOCENTRIC_H

26
src/map/lambertazimuthal.cpp
View File

@ -16,19 +16,17 @@ LambertAzimuthal::LambertAzimuthal(const Ellipsoid *ellipsoid,
_lon0 = deg2rad(longitudeOrigin); _lon0 = deg2rad(longitudeOrigin);
_a = ellipsoid->radius(); _a = ellipsoid->radius();
_es = ellipsoid->es();
_e = sqrt(_es);
_es2 = 2.0 * ellipsoid->flattening() - ellipsoid->flattening() double q0 = (1.0 - _es) * ((sin(lat0) / (1.0 - _es * sin2(lat0)))
* ellipsoid->flattening(); - ((1.0/(2.0*_e)) * log((1.0 - _e * sin(lat0)) / (1.0 + _e
_es = sqrt(_es2);
double q0 = (1.0 - _es2) * ((sin(lat0) / (1.0 - _es2 * sin2(lat0)))
- ((1.0/(2.0*_es)) * log((1.0 - _es * sin(lat0)) / (1.0 + _es
* sin(lat0))))); * sin(lat0)))));
_qP = (1.0 - _es2) * ((1.0 / (1.0 - _es2)) - ((1.0/(2.0*_es)) _qP = (1.0 - _es) * ((1.0 / (1.0 - _es)) - ((1.0/(2.0*_e))
* log((1.0 - _es) / (1.0 + _es)))); * log((1.0 - _e) / (1.0 + _e))));
_beta0 = asin(q0 / _qP); _beta0 = asin(q0 / _qP);
_Rq = _a * sqrt(_qP / 2.0); _Rq = _a * sqrt(_qP / 2.0);
_D = _a * (cos(lat0) / sqrt(1.0 - _es2 * sin2(lat0))) / (_Rq * cos(_beta0)); _D = _a * (cos(lat0) / sqrt(1.0 - _es * sin2(lat0))) / (_Rq * cos(_beta0));
} }
PointD LambertAzimuthal::ll2xy(const Coordinates &c) const PointD LambertAzimuthal::ll2xy(const Coordinates &c) const
@ -36,8 +34,8 @@ PointD LambertAzimuthal::ll2xy(const Coordinates &c) const
double lon = deg2rad(c.lon()); double lon = deg2rad(c.lon());
double lat = deg2rad(c.lat()); double lat = deg2rad(c.lat());
double q = (1.0 - _es2) * ((sin(lat) / (1.0 - _es2 * sin2(lat))) double q = (1.0 - _es) * ((sin(lat) / (1.0 - _es * sin2(lat)))
- ((1.0/(2.0*_es)) * log((1.0 - _es * sin(lat)) / (1.0 + _es - ((1.0/(2.0*_e)) * log((1.0 - _e * sin(lat)) / (1.0 + _e
* sin(lat))))); * sin(lat)))));
double beta = asin(q / _qP); double beta = asin(q / _qP);
double B = _Rq * sqrt(2.0 / (1.0 + sin(_beta0) * sin(beta) + (cos(_beta0) double B = _Rq * sqrt(2.0 / (1.0 + sin(_beta0) * sin(beta) + (cos(_beta0)
@ -52,8 +50,8 @@ PointD LambertAzimuthal::ll2xy(const Coordinates &c) const
Coordinates LambertAzimuthal::xy2ll(const PointD &p) const Coordinates LambertAzimuthal::xy2ll(const PointD &p) const
{ {
double es4 = _es2 * _es2; double es4 = _es * _es;
double es6 = _es2 * es4; double es6 = _es * es4;
double rho = sqrt(sqr((p.x() - _falseEasting) / _D) + sqr(_D * (p.y() double rho = sqrt(sqr((p.x() - _falseEasting) / _D) + sqr(_D * (p.y()
- _falseNorthing))); - _falseNorthing)));
@ -64,7 +62,7 @@ Coordinates LambertAzimuthal::xy2ll(const PointD &p) const
double lon = _lon0 + atan((p.x() - _falseEasting) * sin(C) / (_D * rho double lon = _lon0 + atan((p.x() - _falseEasting) * sin(C) / (_D * rho
* cos(_beta0) * cos(C) - sqr(_D) * (p.y() - _falseNorthing) * sin(_beta0) * cos(_beta0) * cos(C) - sqr(_D) * (p.y() - _falseNorthing) * sin(_beta0)
* sin(C))); * sin(C)));
double lat = betaS + ((_es2/3.0 + 31.0*es4/180.0 + 517.0*es6/5040.0) double lat = betaS + ((_es/3.0 + 31.0*es4/180.0 + 517.0*es6/5040.0)
* sin(2.0*betaS)) + ((23.0*es4/360.0 + 251.0*es6/3780.0) * sin(4.0*betaS)) * sin(2.0*betaS)) + ((23.0*es4/360.0 + 251.0*es6/3780.0) * sin(4.0*betaS))
+ ((761.0*es6/45360.0)*sin(6.0*betaS)); + ((761.0*es6/45360.0)*sin(6.0*betaS));

2
src/map/lambertazimuthal.h
View File

@ -20,7 +20,7 @@ private:
double _lon0; double _lon0;
double _falseNorthing; double _falseNorthing;
double _falseEasting; double _falseEasting;
double _a, _es, _es2, _qP, _beta0, _Rq, _D; double _a, _e, _es, _qP, _beta0, _Rq, _D;
}; };
#endif // LAMBERTAZIMUTHAL_H #endif // LAMBERTAZIMUTHAL_H

21
src/map/lambertconic.cpp
View File

@ -57,8 +57,7 @@ LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin,
double longitudeOrigin, double scale, double falseEasting, double longitudeOrigin, double scale, double falseEasting,
double falseNorthing) double falseNorthing)
{ {
double es2; double e_sin;
double es_sin;
double m0; double m0;
double lat_orig; double lat_orig;
@ -71,16 +70,14 @@ LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin,
_falseEasting = falseEasting; _falseEasting = falseEasting;
_falseNorthing = falseNorthing; _falseNorthing = falseNorthing;
es2 = 2.0 * ellipsoid->flattening() - ellipsoid->flattening() _e = sqrt(ellipsoid->es());
* ellipsoid->flattening(); _e_over_2 = _e / 2.0;
_es = sqrt(es2);
_es_over_2 = _es / 2.0;
_n = sin(lat_orig); _n = sin(lat_orig);
es_sin = _es * sin(lat_orig); e_sin = _e * sin(lat_orig);
m0 = LAMBERT_m(cos(lat_orig), es_sin); m0 = LAMBERT_m(cos(lat_orig), e_sin);
_t0 = LAMBERT1_t(lat_orig, es_sin, _es_over_2); _t0 = LAMBERT1_t(lat_orig, e_sin, _e_over_2);
_rho0 = ellipsoid->radius() * scale * m0 / _n; _rho0 = ellipsoid->radius() * scale * m0 / _n;
@ -97,7 +94,7 @@ PointD LambertConic1::ll2xy(const Coordinates &c) const
if (fabs(fabs(lat) - M_PI_2) > 1.0e-10) { if (fabs(fabs(lat) - M_PI_2) > 1.0e-10) {
t = LAMBERT1_t(lat, _es * sin(lat), _es_over_2); t = LAMBERT1_t(lat, _e * sin(lat), _e_over_2);
rho = _rho0 * pow(t / _t0, _n); rho = _rho0 * pow(t / _t0, _n);
} else } else
rho = 0.0; rho = 0.0;
@ -148,9 +145,9 @@ Coordinates LambertConic1::xy2ll(const PointD &p) const
PHI = M_PI_2 - 2.0 * atan(t); PHI = M_PI_2 - 2.0 * atan(t);
while (fabs(PHI - tempPHI) > tolerance && count) { while (fabs(PHI - tempPHI) > tolerance && count) {
tempPHI = PHI; tempPHI = PHI;
es_sin = _es * sin(PHI); es_sin = _e * sin(PHI);
PHI = M_PI_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin), PHI = M_PI_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin),
_es_over_2)); _e_over_2));
count--; count--;
} }

4
src/map/lambertconic.h
View File

@ -23,8 +23,8 @@ private:
double _falseEasting; double _falseEasting;
double _falseNorthing; double _falseNorthing;
double _es; double _e;
double _es_over_2; double _e_over_2;
double _n; double _n;
double _t0; double _t0;
double _rho0; double _rho0;

19
src/map/mercator.cpp
View File

@ -47,28 +47,29 @@ Defense.
Mercator::Mercator(const Ellipsoid *ellipsoid, double latitudeOrigin, Mercator::Mercator(const Ellipsoid *ellipsoid, double latitudeOrigin,
double longitudeOrigin, double falseEasting, double falseNorthing) double longitudeOrigin, double falseEasting, double falseNorthing)
{ {
double es = ellipsoid->es();
double es2; double es2;
double es3; double es3;
double es4; double es4;
double sin_olat; double sin_olat;
_a = ellipsoid->radius();
_latitudeOrigin = deg2rad(latitudeOrigin); _latitudeOrigin = deg2rad(latitudeOrigin);
_longitudeOrigin = deg2rad(longitudeOrigin); _longitudeOrigin = deg2rad(longitudeOrigin);
if (_longitudeOrigin > M_PI) if (_longitudeOrigin > M_PI)
_longitudeOrigin -= M_2_PI; _longitudeOrigin -= M_2_PI;
_falseNorthing = falseNorthing; _falseNorthing = falseNorthing;
_falseEasting = falseEasting; _falseEasting = falseEasting;
_es = 2 * ellipsoid->flattening() - ellipsoid->flattening()
* ellipsoid->flattening(); _a = ellipsoid->radius();
_e = sqrt(_es); _e = sqrt(es);
sin_olat = sin(_latitudeOrigin); sin_olat = sin(_latitudeOrigin);
_scaleFactor = 1.0 / (sqrt(1.e0 - _es * sin_olat * sin_olat) _scaleFactor = 1.0 / (sqrt(1.e0 - es * sin_olat * sin_olat)
/ cos(_latitudeOrigin)); / cos(_latitudeOrigin));
es2 = _es * _es; es2 = es * es;
es3 = es2 * _es; es3 = es2 * es;
es4 = es3 * _es; es4 = es3 * es;
_ab = _es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0 + 13.e0 * es4 / 360.e0; _ab = es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0 + 13.e0 * es4 / 360.e0;
_bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0 + 811.e0 * es4 / 11520.e0; _bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0 + 811.e0 * es4 / 11520.e0;
_cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0; _cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
_db = 4279.e0 * es4 / 161280.e0; _db = 4279.e0 * es4 / 161280.e0;

2
src/map/mercator.h
View File

@ -17,7 +17,7 @@ public:
virtual Coordinates xy2ll(const PointD &p) const; virtual Coordinates xy2ll(const PointD &p) const;
private: private:
double _a, _es, _e; double _a, _e;
double _latitudeOrigin; double _latitudeOrigin;
double _longitudeOrigin; double _longitudeOrigin;
double _falseNorthing; double _falseNorthing;

6
src/map/transversemercator.cpp
View File

@ -72,11 +72,9 @@ TransverseMercator::TransverseMercator(const Ellipsoid *ellipsoid,
_falseEasting = falseEasting; _falseEasting = falseEasting;
_falseNorthing = falseNorthing; _falseNorthing = falseNorthing;
_es = 2 * ellipsoid->flattening() - ellipsoid->flattening() _es = ellipsoid->es();
* ellipsoid->flattening();
_ebs = (1 / (1 - _es)) - 1; _ebs = (1 / (1 - _es)) - 1;
b = ellipsoid->b();
b = _a * (1 - ellipsoid->flattening());
tn = (_a - b) / (_a + b); tn = (_a - b) / (_a + b);
tn2 = tn * tn; tn2 = tn * tn;

1
src/map/transversemercator.h
View File

@ -24,7 +24,6 @@ private:
double _falseEasting; double _falseEasting;
double _falseNorthing; double _falseNorthing;
double _a; double _a;
double _es; double _es;
double _ebs; double _ebs;
double _ap, _bp, _cp, _dp, _ep; double _ap, _bp, _cp, _dp, _ep;