diff --git a/src/map/albersequal.cpp b/src/map/albersequal.cpp index e300e740..2aa94730 100644 --- a/src/map/albersequal.cpp +++ b/src/map/albersequal.cpp @@ -46,11 +46,11 @@ Defense. #define ONE_MINUS_SQR(x) (1.0 - (x) * (x)) -#define ALBERS_Q(slat, one_minus_sqr_es_sin, es_sin) \ - (_one_minus_es2 * ((slat) / (one_minus_sqr_es_sin) - \ - (1 / (_two_es)) * log((1 - (es_sin)) / (1 + (es_sin))))) -#define ALBERS_M(clat, one_minus_sqr_es_sin) \ - ((clat) / sqrt(one_minus_sqr_es_sin)) +#define ALBERS_Q(slat, one_minus_sqr_e_sin, es_sin) \ + (_one_minus_es * ((slat) / (one_minus_sqr_e_sin) - \ + (1 / (_two_e)) * log((1 - (es_sin)) / (1 + (es_sin))))) +#define ALBERS_M(clat, one_minus_sqr_e_sin) \ + ((clat) / sqrt(one_minus_sqr_e_sin)) 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 m1, m2, sqr_m1, sqr_m2; double q0, q1, q2; - double es_sin, es_sin1, es_sin2; - double one_minus_sqr_es_sin1, one_minus_sqr_es_sin2; + double e_sin, e_sin1, e_sin2; + double one_minus_sqr_e_sin1, one_minus_sqr_e_sin2; double nq0; double sp1, sp2; @@ -75,31 +75,30 @@ AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1, sp2 = deg2rad(standardParallel2); _a2 = ellipsoid->radius() * ellipsoid->radius(); - _es2 = 2 * ellipsoid->flattening() - ellipsoid->flattening() - * ellipsoid->flattening(); - _es = sqrt(_es2); - _one_minus_es2 = 1 - _es2; - _two_es = 2 * _es; + _es = ellipsoid->es(); + _e = sqrt(_es); + _one_minus_es = 1 - _es; + _two_e = 2 * _e; sin_lat = sin(_latitudeOrigin); - es_sin = _es * sin_lat; - q0 = ALBERS_Q(sin_lat, ONE_MINUS_SQR(es_sin), es_sin); + e_sin = _e * sin_lat; + q0 = ALBERS_Q(sin_lat, ONE_MINUS_SQR(e_sin), e_sin); sin_lat1 = sin(sp1); cos_lat1 = cos(sp1); - es_sin1 = _es * sin_lat1; - one_minus_sqr_es_sin1 = ONE_MINUS_SQR(es_sin1); - m1 = ALBERS_M(cos_lat1, one_minus_sqr_es_sin1); - q1 = ALBERS_Q(sin_lat1, one_minus_sqr_es_sin1, es_sin1); + e_sin1 = _e * sin_lat1; + one_minus_sqr_e_sin1 = ONE_MINUS_SQR(e_sin1); + m1 = ALBERS_M(cos_lat1, one_minus_sqr_e_sin1); + q1 = ALBERS_Q(sin_lat1, one_minus_sqr_e_sin1, e_sin1); sqr_m1 = m1 * m1; if (fabs(sp1 - sp2) > 1.0e-10) { sin_lat2 = sin(sp2); cos_lat2 = cos(sp2); - es_sin2 = _es * sin_lat2; - one_minus_sqr_es_sin2 = ONE_MINUS_SQR(es_sin2); - m2 = ALBERS_M(cos_lat2, one_minus_sqr_es_sin2); - q2 = ALBERS_Q(sin_lat2, one_minus_sqr_es_sin2, es_sin2); + e_sin2 = _e * sin_lat2; + one_minus_sqr_e_sin2 = ONE_MINUS_SQR(e_sin2); + m2 = ALBERS_M(cos_lat2, one_minus_sqr_e_sin2); + q2 = ALBERS_Q(sin_lat2, one_minus_sqr_e_sin2, e_sin2); sqr_m2 = m2 * m2; _n = (sqr_m1 - sqr_m2) / (q2 - q1); } else @@ -115,7 +114,7 @@ PointD AlbersEqual::ll2xy(const Coordinates &c) const { double dlam; double sin_lat; - double es_sin; + double e_sin; double q; double rho; double theta; @@ -129,8 +128,8 @@ PointD AlbersEqual::ll2xy(const Coordinates &c) const dlam += M_2_PI; sin_lat = sin(deg2rad(c.lat())); - es_sin = _es * sin_lat; - q = ALBERS_Q(sin_lat, ONE_MINUS_SQR(es_sin), es_sin); + e_sin = _e * sin_lat; + q = ALBERS_Q(sin_lat, ONE_MINUS_SQR(e_sin), e_sin); nq = _n * q; rho = (_C < nq) ? 0 : _a_over_n * sqrt(_C - nq); theta = _n * dlam; @@ -147,7 +146,7 @@ Coordinates AlbersEqual::xy2ll(const PointD &p) const double rho, rho_n; double phi, delta_phi = 1.0; double sin_phi; - double es_sin, one_minus_sqr_es_sin; + double e_sin, one_minus_sqr_e_sin; double theta = 0.0; int count = 30; double tolerance = 4.85e-10; @@ -170,7 +169,7 @@ Coordinates AlbersEqual::xy2ll(const PointD &p) const theta = atan2(dx, rho0_minus_dy); rho_n = rho * _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) { q_over_2 = q / 2.0; if (q_over_2 > 1.0) @@ -179,17 +178,17 @@ Coordinates AlbersEqual::xy2ll(const PointD &p) const lat = -M_PI_2; else { phi = asin(q_over_2); - if (_es < 1.0e-10) + if (_e < 1.0e-10) lat = phi; else { while ((fabs(delta_phi) > tolerance) && count) { sin_phi = sin(phi); - es_sin = _es * sin_phi; - one_minus_sqr_es_sin = ONE_MINUS_SQR(es_sin); - delta_phi = (one_minus_sqr_es_sin * one_minus_sqr_es_sin) - / (2.0 * cos(phi)) * (q / (_one_minus_es2) - sin_phi - / one_minus_sqr_es_sin + (log((1.0 - es_sin) - / (1.0 + es_sin)) / (_two_es))); + e_sin = _e * sin_phi; + one_minus_sqr_e_sin = ONE_MINUS_SQR(e_sin); + delta_phi = (one_minus_sqr_e_sin * one_minus_sqr_e_sin) + / (2.0 * cos(phi)) * (q / (_one_minus_es) - sin_phi + / one_minus_sqr_e_sin + (log((1.0 - e_sin) + / (1.0 + e_sin)) / (_two_e))); phi += delta_phi; count --; } diff --git a/src/map/albersequal.h b/src/map/albersequal.h index c555ef7b..ef03d2c2 100644 --- a/src/map/albersequal.h +++ b/src/map/albersequal.h @@ -27,11 +27,11 @@ private: double _rho0; double _C; double _n; + double _e; double _es; - double _es2; double _a_over_n; - double _one_minus_es2; - double _two_es; + double _one_minus_es; + double _two_e; }; #endif // ALBERSEQUAL_H diff --git a/src/map/datum.cpp b/src/map/datum.cpp index 98dab5d9..97b77f1b 100644 --- a/src/map/datum.cpp +++ b/src/map/datum.cpp @@ -113,7 +113,7 @@ Coordinates Datum::toWGS84(const Coordinates &c) const switch (_transformation) { case Helmert: return Geocentric::toGeodetic(helmert(Geocentric::fromGeodetic(c, - *this)), WGS84); + ellipsoid())), WGS84.ellipsoid()); case Molodensky: return molodensky(c, *this, WGS84); default: @@ -126,7 +126,7 @@ Coordinates Datum::fromWGS84(const Coordinates &c) const switch (_transformation) { case Helmert: return Geocentric::toGeodetic(helmertr(Geocentric::fromGeodetic(c, - WGS84)), *this); + WGS84.ellipsoid())), ellipsoid()); case Molodensky: return molodensky(c, WGS84, *this); default: diff --git a/src/map/ellipsoid.cpp b/src/map/ellipsoid.cpp index 41ee0da7..0689b4e0 100644 --- a/src/map/ellipsoid.cpp +++ b/src/map/ellipsoid.cpp @@ -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 QDebug operator<<(QDebug dbg, const Ellipsoid &ellipsoid) { diff --git a/src/map/ellipsoid.h b/src/map/ellipsoid.h index 911c934f..930c45ae 100644 --- a/src/map/ellipsoid.h +++ b/src/map/ellipsoid.h @@ -9,12 +9,13 @@ class Ellipsoid { public: - Ellipsoid() : _radius(NAN), _flattening(NAN) {} - Ellipsoid(double radius, double flattening) - : _radius(radius), _flattening(flattening) {} + Ellipsoid() : _radius(NAN), _flattening(NAN), _es(NAN), _b(NAN) {} + Ellipsoid(double radius, double flattening); double radius() const {return _radius;} double flattening() const {return _flattening;} + double es() const {return _es;} + double b() const {return _b;} bool isNull() const {return (std::isnan(_radius) && std::isnan(_flattening));} @@ -27,6 +28,7 @@ public: private: double _radius; double _flattening; + double _es, _b; static QMap WGS84(); static QMap _ellipsoids; diff --git a/src/map/geocentric.cpp b/src/map/geocentric.cpp index 545d1163..6b87e579 100644 --- a/src/map/geocentric.cpp +++ b/src/map/geocentric.cpp @@ -41,42 +41,34 @@ Defense. */ -#include "datum.h" +#include "ellipsoid.h" #include "geocentric.h" #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 lon = deg2rad(c.lon()); double slat = sin(lat); double slat2 = slat * slat; 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) lon -= M_2_PI; 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 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; - + double ep2 = (1.0 / (1.0 - e->es())) - 1.0; bool pole = false; double lat, lon; - if (p.x() == 0.0) { if (p.y() > 0) lon = M_PI_2; @@ -102,8 +94,8 @@ Coordinates Geocentric::toGeodetic(const Point3D &p, const Datum &datum) double Sin_B0 = T0 / S0; double Cos_B0 = W / S0; double Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0; - double T1 = p.z() + b * ep2 * Sin3_B0; - double Sum = W - e->radius() * e2 * Cos_B0 * Cos_B0 * Cos_B0; + double T1 = p.z() + e->b() * ep2 * Sin3_B0; + double Sum = W - e->radius() * e->es() * Cos_B0 * Cos_B0 * Cos_B0; double S1 = sqrt(T1*T1 + Sum * Sum); double Sin_p1 = T1 / S1; double Cos_p1 = Sum / S1; diff --git a/src/map/geocentric.h b/src/map/geocentric.h index 1f93e2c0..58119d39 100644 --- a/src/map/geocentric.h +++ b/src/map/geocentric.h @@ -4,7 +4,7 @@ #include #include "common/coordinates.h" -class Datum; +class Ellipsoid; class Point3D { public: @@ -22,8 +22,8 @@ private: }; namespace Geocentric { - Point3D fromGeodetic(const Coordinates &c, const Datum &datum); - Coordinates toGeodetic(const Point3D &p, const Datum &datum); + Point3D fromGeodetic(const Coordinates &c, const Ellipsoid *e); + Coordinates toGeodetic(const Point3D &p, const Ellipsoid *e); } #endif // GEOCENTRIC_H diff --git a/src/map/lambertazimuthal.cpp b/src/map/lambertazimuthal.cpp index a76082ee..6255c12f 100644 --- a/src/map/lambertazimuthal.cpp +++ b/src/map/lambertazimuthal.cpp @@ -16,19 +16,17 @@ LambertAzimuthal::LambertAzimuthal(const Ellipsoid *ellipsoid, _lon0 = deg2rad(longitudeOrigin); _a = ellipsoid->radius(); + _es = ellipsoid->es(); + _e = sqrt(_es); - _es2 = 2.0 * ellipsoid->flattening() - ellipsoid->flattening() - * ellipsoid->flattening(); - _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 + double q0 = (1.0 - _es) * ((sin(lat0) / (1.0 - _es * sin2(lat0))) + - ((1.0/(2.0*_e)) * log((1.0 - _e * sin(lat0)) / (1.0 + _e * sin(lat0))))); - _qP = (1.0 - _es2) * ((1.0 / (1.0 - _es2)) - ((1.0/(2.0*_es)) - * log((1.0 - _es) / (1.0 + _es)))); + _qP = (1.0 - _es) * ((1.0 / (1.0 - _es)) - ((1.0/(2.0*_e)) + * log((1.0 - _e) / (1.0 + _e)))); _beta0 = asin(q0 / _qP); _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 @@ -36,8 +34,8 @@ PointD LambertAzimuthal::ll2xy(const Coordinates &c) const double lon = deg2rad(c.lon()); double lat = deg2rad(c.lat()); - double q = (1.0 - _es2) * ((sin(lat) / (1.0 - _es2 * sin2(lat))) - - ((1.0/(2.0*_es)) * log((1.0 - _es * sin(lat)) / (1.0 + _es + double q = (1.0 - _es) * ((sin(lat) / (1.0 - _es * sin2(lat))) + - ((1.0/(2.0*_e)) * log((1.0 - _e * sin(lat)) / (1.0 + _e * sin(lat))))); double beta = asin(q / _qP); 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 { - double es4 = _es2 * _es2; - double es6 = _es2 * es4; + double es4 = _es * _es; + double es6 = _es * es4; double rho = sqrt(sqr((p.x() - _falseEasting) / _D) + sqr(_D * (p.y() - _falseNorthing))); @@ -64,7 +62,7 @@ Coordinates LambertAzimuthal::xy2ll(const PointD &p) const double lon = _lon0 + atan((p.x() - _falseEasting) * sin(C) / (_D * rho * cos(_beta0) * cos(C) - sqr(_D) * (p.y() - _falseNorthing) * sin(_beta0) * 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)) + ((761.0*es6/45360.0)*sin(6.0*betaS)); diff --git a/src/map/lambertazimuthal.h b/src/map/lambertazimuthal.h index 97d6016b..eecd329b 100644 --- a/src/map/lambertazimuthal.h +++ b/src/map/lambertazimuthal.h @@ -20,7 +20,7 @@ private: double _lon0; double _falseNorthing; double _falseEasting; - double _a, _es, _es2, _qP, _beta0, _Rq, _D; + double _a, _e, _es, _qP, _beta0, _Rq, _D; }; #endif // LAMBERTAZIMUTHAL_H diff --git a/src/map/lambertconic.cpp b/src/map/lambertconic.cpp index e9efc7b4..b1f114ed 100644 --- a/src/map/lambertconic.cpp +++ b/src/map/lambertconic.cpp @@ -57,8 +57,7 @@ LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin, double longitudeOrigin, double scale, double falseEasting, double falseNorthing) { - double es2; - double es_sin; + double e_sin; double m0; double lat_orig; @@ -71,16 +70,14 @@ LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin, _falseEasting = falseEasting; _falseNorthing = falseNorthing; - es2 = 2.0 * ellipsoid->flattening() - ellipsoid->flattening() - * ellipsoid->flattening(); - _es = sqrt(es2); - _es_over_2 = _es / 2.0; + _e = sqrt(ellipsoid->es()); + _e_over_2 = _e / 2.0; _n = sin(lat_orig); - es_sin = _es * sin(lat_orig); - m0 = LAMBERT_m(cos(lat_orig), es_sin); - _t0 = LAMBERT1_t(lat_orig, es_sin, _es_over_2); + e_sin = _e * sin(lat_orig); + m0 = LAMBERT_m(cos(lat_orig), e_sin); + _t0 = LAMBERT1_t(lat_orig, e_sin, _e_over_2); _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) { - 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); } else rho = 0.0; @@ -148,9 +145,9 @@ Coordinates LambertConic1::xy2ll(const PointD &p) const PHI = M_PI_2 - 2.0 * atan(t); while (fabs(PHI - tempPHI) > tolerance && count) { 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), - _es_over_2)); + _e_over_2)); count--; } diff --git a/src/map/lambertconic.h b/src/map/lambertconic.h index 32b84c26..de1e9685 100644 --- a/src/map/lambertconic.h +++ b/src/map/lambertconic.h @@ -23,8 +23,8 @@ private: double _falseEasting; double _falseNorthing; - double _es; - double _es_over_2; + double _e; + double _e_over_2; double _n; double _t0; double _rho0; diff --git a/src/map/mercator.cpp b/src/map/mercator.cpp index e4785e39..c42b3fec 100644 --- a/src/map/mercator.cpp +++ b/src/map/mercator.cpp @@ -47,28 +47,29 @@ Defense. Mercator::Mercator(const Ellipsoid *ellipsoid, double latitudeOrigin, double longitudeOrigin, double falseEasting, double falseNorthing) { + double es = ellipsoid->es(); double es2; double es3; double es4; double sin_olat; - _a = ellipsoid->radius(); _latitudeOrigin = deg2rad(latitudeOrigin); _longitudeOrigin = deg2rad(longitudeOrigin); if (_longitudeOrigin > M_PI) _longitudeOrigin -= M_2_PI; _falseNorthing = falseNorthing; _falseEasting = falseEasting; - _es = 2 * ellipsoid->flattening() - ellipsoid->flattening() - * ellipsoid->flattening(); - _e = sqrt(_es); + + _a = ellipsoid->radius(); + _e = sqrt(es); + 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)); - es2 = _es * _es; - es3 = es2 * _es; - es4 = es3 * _es; - _ab = _es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0 + 13.e0 * es4 / 360.e0; + es2 = es * es; + es3 = es2 * es; + es4 = es3 * es; + _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; _cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0; _db = 4279.e0 * es4 / 161280.e0; diff --git a/src/map/mercator.h b/src/map/mercator.h index d55db268..39d0390b 100644 --- a/src/map/mercator.h +++ b/src/map/mercator.h @@ -17,7 +17,7 @@ public: virtual Coordinates xy2ll(const PointD &p) const; private: - double _a, _es, _e; + double _a, _e; double _latitudeOrigin; double _longitudeOrigin; double _falseNorthing; diff --git a/src/map/transversemercator.cpp b/src/map/transversemercator.cpp index 884c0c1e..65f084a2 100644 --- a/src/map/transversemercator.cpp +++ b/src/map/transversemercator.cpp @@ -72,11 +72,9 @@ TransverseMercator::TransverseMercator(const Ellipsoid *ellipsoid, _falseEasting = falseEasting; _falseNorthing = falseNorthing; - _es = 2 * ellipsoid->flattening() - ellipsoid->flattening() - * ellipsoid->flattening(); + _es = ellipsoid->es(); _ebs = (1 / (1 - _es)) - 1; - - b = _a * (1 - ellipsoid->flattening()); + b = ellipsoid->b(); tn = (_a - b) / (_a + b); tn2 = tn * tn; diff --git a/src/map/transversemercator.h b/src/map/transversemercator.h index 18c2ea38..9e9e2f18 100644 --- a/src/map/transversemercator.h +++ b/src/map/transversemercator.h @@ -24,7 +24,6 @@ private: double _falseEasting; double _falseNorthing; double _a; - double _es; double _ebs; double _ap, _bp, _cp, _dp, _ep;