Code cleanup

src/datum.cpp

@@ -1,5 +1,7 @@
+#include <cmath>
 #include <QFile>
 #include "wgs84.h"
+#include "rd.h"
 #include "datum.h"
 
 
@@ -81,3 +83,39 @@ bool Datum::loadList(const QString &path)
 
 	return true;
 }
+
+// Abridged Molodensky transformation
+Coordinates Datum::toWGS84(const Coordinates &c) const
+{
+	if (_ellipsoid.radius() == WGS84_RADIUS
+	  && _ellipsoid.flattening() == WGS84_FLATTENING
+	  && _dx == 0 && _dy == 0 && _dz == 0)
+		return c;
+
+	double rlat = deg2rad(c.lat());
+	double rlon = deg2rad(c.lon());
+
+	double slat = sin(rlat);
+	double clat = cos(rlat);
+	double slon = sin(rlon);
+	double clon = cos(rlon);
+	double ssqlat = slat * slat;
+
+	double from_f = ellipsoid().flattening();
+	double df = WGS84_FLATTENING - from_f;
+	double from_a = ellipsoid().radius();
+	double da = WGS84_RADIUS - from_a;
+	double from_esq = ellipsoid().flattening()
+	  * (2.0 - ellipsoid().flattening());
+	double adb = 1.0 / (1.0 - from_f);
+	double rn = from_a / sqrt(1 - from_esq * ssqlat);
+	double rm = from_a * (1 - from_esq) / pow((1 - from_esq * ssqlat), 1.5);
+	double from_h = 0.0;
+
+	double dlat = (-dx() * slat * clon - dy() * slat * slon + dz() * clat + da
+	  * rn * from_esq * slat * clat / from_a + df * (rm * adb + rn / adb) * slat
+	  * clat) / (rm + from_h);
+	double dlon = (-dx() * slon + dy() * clon) / ((rn + from_h) * clat);
+
+	return Coordinates(c.lon() + rad2deg(dlon), c.lat() + rad2deg(dlat));
+}

src/datum.h

@@ -3,6 +3,7 @@
 
 #include <QMap>
 #include "ellipsoid.h"
+#include "coordinates.h"
 
 class Datum
 {
@@ -17,8 +18,8 @@ public:
 	double dz() const {return _dz;}
 
 	bool isNull() const {return _ellipsoid.isNull();}
-	bool isWGS84() const
+
-	  {return _ellipsoid.isWGS84() && _dx == 0 && _dy == 0 && _dz == 0;}
+	Coordinates toWGS84(const Coordinates &c) const;
 
 	static bool loadList(const QString &path);
 	static const QString &errorString() {return _errorString;}

src/ellipsoid.h

@@ -3,7 +3,6 @@
 
 #include <QString>
 #include <QMap>
-#include "wgs84.h"
 
 class Ellipsoid
 {
@@ -16,8 +15,6 @@ public:
 	double flattening() const {return _flattening;}
 
 	bool isNull() const {return _radius < 0 || _flattening < 0;}
-	bool isWGS84() const
-	  {return _radius == WGS84_RADIUS && _flattening == WGS84_FLATTENING;}
 
 	static bool loadList(const QString &path);
 	static const QString &errorString() {return _errorString;}

src/offlinemap.cpp

@@ -23,39 +23,6 @@
 #include "offlinemap.h"
 
 
-// Abridged Molodensky transformation
-static Coordinates toWGS84(Coordinates c, const Datum &datum)
-{
-	double rlat = deg2rad(c.lat());
-	double rlon = deg2rad(c.lon());
-
-	double slat = sin(rlat);
-	double clat = cos(rlat);
-	double slon = sin(rlon);
-	double clon = cos(rlon);
-	double ssqlat = slat * slat;
-
-	double from_f = datum.ellipsoid().flattening();
-	double df = WGS84_FLATTENING - from_f;
-	double from_a = datum.ellipsoid().radius();
-	double da = WGS84_RADIUS - from_a;
-	double from_esq = datum.ellipsoid().flattening()
-	  * (2.0 - datum.ellipsoid().flattening());
-	double adb = 1.0 / (1.0 - from_f);
-	double rn = from_a / sqrt(1 - from_esq * ssqlat);
-	double rm = from_a * (1 - from_esq) / pow((1 - from_esq * ssqlat), 1.5);
-	double from_h = 0.0;
-
-	double dlat = (-datum.dx() * slat * clon - datum.dy() * slat * slon
-	  + datum.dz() * clat + da * rn * from_esq * slat * clat / from_a + df
-	  * (rm * adb + rn / adb) * slat * clat) / (rm + from_h);
-	double dlon = (-datum.dx() * slon + datum.dy() * clon) / ((rn + from_h)
-	  * clat);
-
-	return Coordinates(c.lon() + rad2deg(dlon), c.lat() + rad2deg(dlat));
-}
-
-
 int OfflineMap::parse(QIODevice &device, QList<ReferencePoint> &points,
   QString &projection, ProjectionSetup &setup, QString &datum)
 {
@@ -263,15 +230,10 @@ bool OfflineMap::createProjection(const QString &datum,
 	}
 
 	for (int i = 0; i < points.size(); i++) {
-		if (points.at(i).ll.isNull()) {
+		if (points.at(i).ll.isNull())
-			if (d.isWGS84())
+			points[i].ll = d.toWGS84(_projection->xy2ll(points.at(i).pp));
-				points[i].ll = _projection->xy2ll(points.at(i).pp);
+		else
-			else
+			points[i].ll = d.toWGS84(points[i].ll);
-				points[i].ll = toWGS84(_projection->xy2ll(points.at(i).pp), d);
-		} else {
-			if (!d.isWGS84())
-				points[i].ll = toWGS84(points[i].ll, d);
-		}
 	}
 
 	return true;