From a486f9e78d06b30ce7951c1036f42a6958990481 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Martin=20T=C5=AFma?= <tumic@cbox.cz>
Date: Thu, 17 May 2018 22:41:56 +0200
Subject: [PATCH] Datum shift optimization (precompute ellipsoid constants) +
 projections code unification

---
 src/map/albersequal.cpp        | 67 +++++++++++++++++-----------------
 src/map/albersequal.h          |  6 +--
 src/map/datum.cpp              |  4 +-
 src/map/ellipsoid.cpp          |  7 ++++
 src/map/ellipsoid.h            |  8 ++--
 src/map/geocentric.cpp         | 24 ++++--------
 src/map/geocentric.h           |  6 +--
 src/map/lambertazimuthal.cpp   | 26 ++++++-------
 src/map/lambertazimuthal.h     |  2 +-
 src/map/lambertconic.cpp       | 21 +++++------
 src/map/lambertconic.h         |  4 +-
 src/map/mercator.cpp           | 19 +++++-----
 src/map/mercator.h             |  2 +-
 src/map/transversemercator.cpp |  6 +--
 src/map/transversemercator.h   |  1 -
 15 files changed, 98 insertions(+), 105 deletions(-)

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<int, Ellipsoid> WGS84();
 	static QMap<int, Ellipsoid> _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 <cmath>
 #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;