From a6389f0174a07703de192cbd0323af497c4885d2 Mon Sep 17 00:00:00 2001
From: =?UTF-8?q?Martin=20T=C5=AFma?= <tumic@cbox.cz>
Date: Fri, 18 May 2018 22:10:43 +0200
Subject: [PATCH] Precompute e2s as well

---
 src/map/ellipsoid.cpp  | 1 +
 src/map/ellipsoid.h    | 3 ++-
 src/map/geocentric.cpp | 3 +--
 3 files changed, 4 insertions(+), 3 deletions(-)

diff --git a/src/map/ellipsoid.cpp b/src/map/ellipsoid.cpp
index 0689b4e0..033eee09 100644
--- a/src/map/ellipsoid.cpp
+++ b/src/map/ellipsoid.cpp
@@ -70,6 +70,7 @@ Ellipsoid::Ellipsoid(double radius, double flattening)
   : _radius(radius), _flattening(flattening)
 {
 	_es = 2.0 * flattening - flattening * flattening;
+	_e2s = (1.0 / (1.0 - _es)) - 1.0;
 	_b = radius * (1.0 - flattening);
 }
 
diff --git a/src/map/ellipsoid.h b/src/map/ellipsoid.h
index 930c45ae..bb4e6a5b 100644
--- a/src/map/ellipsoid.h
+++ b/src/map/ellipsoid.h
@@ -15,6 +15,7 @@ public:
 	double radius() const {return _radius;}
 	double flattening() const {return _flattening;}
 	double es() const {return _es;}
+	double e2s() const {return _e2s;}
 	double b() const {return _b;}
 
 	bool isNull() const
@@ -28,7 +29,7 @@ public:
 private:
 	double _radius;
 	double _flattening;
-	double _es, _b;
+	double _es, _e2s, _b;
 
 	static QMap<int, Ellipsoid> WGS84();
 	static QMap<int, Ellipsoid> _ellipsoids;
diff --git a/src/map/geocentric.cpp b/src/map/geocentric.cpp
index 6b87e579..040b85c5 100644
--- a/src/map/geocentric.cpp
+++ b/src/map/geocentric.cpp
@@ -65,7 +65,6 @@ Point3D Geocentric::fromGeodetic(const Coordinates &c, const Ellipsoid *e)
 
 Coordinates Geocentric::toGeodetic(const Point3D &p, const Ellipsoid *e)
 {
-	double ep2 = (1.0 / (1.0 - e->es())) - 1.0;
 	bool pole = false;
 	double lat, lon;
 
@@ -94,7 +93,7 @@ Coordinates Geocentric::toGeodetic(const Point3D &p, const Ellipsoid *e)
 	double Sin_B0 = T0 / S0;
 	double Cos_B0 = W / S0;
 	double Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
-	double T1 = p.z() + e->b() * ep2 * Sin3_B0;
+	double T1 = p.z() + e->b() * e->e2s() * 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;