Project structure cleanup

src/map/proj/albersequal.cpp

@@ -0,0 +1,235 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include "map/ellipsoid.h"
+#include "albersequal.h"
+
+
+#define ONE_MINUS_SQR(x) (1.0 - (x) * (x))
+#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,
+  double standardParallel2, double latitudeOrigin, double longitudeOrigin,
+  double falseEasting, double falseNorthing)
+{
+	double sin_lat, sin_lat1, sin_lat2, cos_lat1, cos_lat2;
+	double m1, m2, sqr_m1, sqr_m2;
+	double q0, q1, q2;
+	double e_sin, e_sin1, e_sin2;
+	double one_minus_sqr_e_sin1, one_minus_sqr_e_sin2;
+	double nq0;
+	double sp1, sp2;
+
+
+	_latitudeOrigin = deg2rad(latitudeOrigin);
+	_longitudeOrigin = deg2rad(longitudeOrigin);
+	_falseEasting = falseEasting;
+	_falseNorthing = falseNorthing;
+
+	sp1 = deg2rad(standardParallel1);
+	sp2 = deg2rad(standardParallel2);
+
+	_a2 = ellipsoid.radius() * ellipsoid.radius();
+	_es = ellipsoid.es();
+	_e = sqrt(_es);
+	_one_minus_es = 1 - _es;
+	_two_e = 2 * _e;
+
+	sin_lat = sin(_latitudeOrigin);
+	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);
+	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);
+		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
+		_n = sin_lat1;
+
+	_c = sqr_m1 + _n * q1;
+	_a_over_n = ellipsoid.radius() / _n;
+	nq0 = _n * q0;
+	_rho0 = (_c < nq0) ? 0 : _a_over_n * sqrt(_c - nq0);
+}
+
+PointD AlbersEqual::ll2xy(const Coordinates &c) const
+{
+	double dlam;
+	double sin_lat;
+	double e_sin;
+	double q;
+	double rho;
+	double theta;
+	double nq;
+
+
+	dlam = deg2rad(c.lon()) - _longitudeOrigin;
+	if (dlam > M_PI)
+		dlam -= 2 * M_PI;
+	if (dlam < -M_PI)
+		dlam += 2 * M_PI;
+
+	sin_lat = sin(deg2rad(c.lat()));
+	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;
+
+	return PointD(rho * sin(theta) + _falseEasting,
+	  _rho0 - rho * cos(theta) + _falseNorthing);
+}
+
+Coordinates AlbersEqual::xy2ll(const PointD &p) const
+{
+	double dy, dx;
+	double rho0_minus_dy;
+	double q, qc, q_over_2;
+	double rho, rho_n;
+	double phi, delta_phi = 1.0;
+	double sin_phi;
+	double e_sin, one_minus_sqr_e_sin;
+	double theta = 0.0;
+	int count = 30;
+	double tolerance = 4.85e-10;
+	double lat, lon;
+
+
+	dy = p.y() - _falseNorthing;
+	dx = p.x() - _falseEasting;
+
+	rho0_minus_dy = _rho0 - dy;
+	rho = sqrt(dx * dx + rho0_minus_dy * rho0_minus_dy);
+
+	if (_n < 0) {
+		rho *= -1.0;
+		dx *= -1.0;
+		rho0_minus_dy *= -1.0;
+	}
+
+	if (rho != 0.0)
+		theta = atan2(dx, rho0_minus_dy);
+	rho_n = rho * _n;
+	q = (_c - (rho_n * rho_n) / _a2) / _n;
+	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)
+			lat = M_PI_2;
+		else if (q_over_2 < -1.0)
+			lat = -M_PI_2;
+		else {
+			phi = asin(q_over_2);
+			if (_e < 1.0e-10)
+				lat = phi;
+			else  {
+				while ((fabs(delta_phi) > tolerance) && count) {
+					sin_phi = sin(phi);
+					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 --;
+				}
+
+				lat = phi;
+			}
+
+			if (lat > M_PI_2)
+				lat = M_PI_2;
+			else if (lat < -M_PI_2)
+				lat = -M_PI_2;
+		}
+	} else {
+		if (q >= 0.0)
+			lat = M_PI_2;
+		else
+			lat = -M_PI_2;
+	}
+
+	lon = _longitudeOrigin + theta / _n;
+
+	if (lon > M_PI)
+		lon -= 2 * M_PI;
+	if (lon < -M_PI)
+		lon += 2 * M_PI;
+
+	if (lon > M_PI)
+		lon = M_PI;
+	else if (lon < -M_PI)
+		lon = -M_PI;
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}
+
+bool AlbersEqual::operator==(const CT &ct) const
+{
+	const AlbersEqual *other = dynamic_cast<const AlbersEqual*>(&ct);
+	return (other != 0 && _latitudeOrigin == other->_latitudeOrigin
+	  && _longitudeOrigin == other->_longitudeOrigin
+	  && _falseEasting == other->_falseEasting
+	  && _falseNorthing == other->_falseNorthing && _a2 == other->_a2
+	  && _es == other->_es && _rho0 == other->_rho0 && _c == other->_c
+	  && _n == other->_n);
+}

src/map/proj/albersequal.h

@@ -0,0 +1,38 @@
+#ifndef ALBERSEQUAL_H
+#define ALBERSEQUAL_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class AlbersEqual : public CT
+{
+public:
+	AlbersEqual(const Ellipsoid &ellipsoid, double standardParallel1,
+	  double standardParallel2, double latitudeOrigin, double longitudeOrigin,
+	  double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new AlbersEqual(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _latitudeOrigin;
+	double _longitudeOrigin;
+	double _falseEasting;
+	double _falseNorthing;
+
+	double _a2;
+	double _rho0;
+	double _c;
+	double _n;
+	double _e;
+	double _es;
+	double _a_over_n;
+	double _one_minus_es;
+	double _two_e;
+};
+
+#endif // ALBERSEQUAL_H

src/map/proj/krovak.cpp

@@ -0,0 +1,85 @@
+#include "map/ellipsoid.h"
+#include "krovak.h"
+
+Krovak::Krovak(const Ellipsoid &ellipsoid, double standardParallel,
+  double azimuth, double scale, double centerLatitude, double longitudeOrigin,
+  double falseEasting, double falseNorthing)
+{
+	double phiC = deg2rad(centerLatitude);
+	double sinPhiC = sin(phiC);
+	double sinPhiC2 = sinPhiC * sinPhiC;
+	double cosPhiC = cos(phiC);
+	double cosPhiC2 = cosPhiC * cosPhiC;
+	double cosPhiC4 = cosPhiC2 * cosPhiC2;
+	double alphaC = deg2rad(azimuth);
+
+	_phiP = deg2rad(standardParallel);
+	_e = sqrt(ellipsoid.es());
+	_a = ellipsoid.radius() * sqrt(1.0 - ellipsoid.es())
+	  / (1.0 - ellipsoid.es() * sinPhiC2);
+	_b = sqrt(1.0 + (ellipsoid.es() * cosPhiC4 / (1.0 - ellipsoid.es())));
+	double gamma0 = asin(sinPhiC / _b);
+	_t0 = tan(M_PI_4 + gamma0 / 2.0) * pow((1.0 + _e * sinPhiC) /
+	  (1.0 - _e * sinPhiC), _e*_b / 2.0) / pow(tan(M_PI_4 + phiC/2.0), _b);
+	_n = sin(_phiP);
+	_r0 = scale * _a / tan(_phiP);
+	_fe = falseEasting;
+	_fn = falseNorthing;
+	_cosAlphaC = cos(alphaC);
+	_sinAlphaC = sin(alphaC);
+	_lambda0 = deg2rad(longitudeOrigin);
+}
+
+PointD Krovak::ll2xy(const Coordinates &c) const
+{
+	double phi = deg2rad(c.lat());
+	double lambda = deg2rad(c.lon());
+	double eSinPhi = _e * sin(phi);
+	double U = 2.0 * (atan(_t0 * pow(tan(phi/2.0 + M_PI_4), _b)
+	  / pow((1.0 + eSinPhi) / (1.0 - eSinPhi), _e * _b/2.0)) - M_PI_4);
+	double cosU = cos(U);
+	double V = _b * (_lambda0 - lambda);
+	double T = asin(_cosAlphaC * sin(U) + _sinAlphaC * cosU * cos(V));
+	double D = asin(cosU * sin(V) / cos(T));
+	double theta = _n * D;
+	double r = _r0 * pow(tan(M_PI_4 + _phiP/2.0), _n)
+	  / pow(tan(T/2.0 + M_PI_4), _n);
+
+	return PointD(r * sin(theta) + _fe, r * cos(theta) + _fn);
+}
+
+Coordinates Krovak::xy2ll(const PointD &p) const
+{
+	double Xp = p.y() - _fn;
+	double Yp = p.x() - _fe;
+	double Xp2 = Xp * Xp;
+	double Yp2 = Yp * Yp;
+	double r = sqrt(Xp2 + Yp2);
+	double theta = atan(Yp / Xp);
+	double D = theta / sin(_phiP);
+	double T = 2.0 * (atan(pow(_r0 / r, 1.0/_n) * tan(M_PI_4 + _phiP/2.0))
+	  - M_PI_4);
+	double U = asin(_cosAlphaC * sin(T) - _sinAlphaC * cos(T) * cos(D));
+	double V = asin(cos(T) * sin(D) / cos(U));
+	double phi = U;
+	for (int i = 0; i < 3; i++)
+		phi = 2.0 * (atan(pow(_t0, -1.0/_b) * pow(tan(U/2.0 + M_PI_4), 1.0/_b)
+		  * pow((1.0 + _e * sin(phi))/(1.0 - _e * sin(phi)), _e/2.0)) - M_PI_4);
+
+	return Coordinates(rad2deg(_lambda0 - V/_b), rad2deg(phi));
+}
+
+bool Krovak::operator==(const CT &ct) const
+{
+	const Krovak *other = dynamic_cast<const Krovak*>(&ct);
+	return (other != 0 && _e == other->_e && _a == other->_a && _b == other->_b
+	  && _t0 == other->_t0 && _r0 == other->_r0 && _phiP == other->_phiP
+	  && _cosAlphaC == other->_cosAlphaC && _sinAlphaC == other->_sinAlphaC
+	  && _lambda0 == other->_lambda0 && _fe == other->_fe && _fn == other->_fn);
+}
+
+bool KrovakNE::operator==(const CT &ct) const
+{
+	const KrovakNE *other = dynamic_cast<const KrovakNE*>(&ct);
+	return (other != 0 && _k == other->_k);
+}

src/map/proj/krovak.h

@@ -0,0 +1,47 @@
+#ifndef KROVAK_H
+#define KROVAK_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class Krovak : public CT
+{
+public:
+	Krovak(const Ellipsoid &ellipsoid, double standardParallel,
+	  double azimuth, double scale, double centerLatitude,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new Krovak(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _e, _a, _b, _t0, _n, _r0, _phiP;
+	double _cosAlphaC, _sinAlphaC, _lambda0, _fe, _fn;
+};
+
+class KrovakNE : public CT
+{
+public:
+	KrovakNE(const Ellipsoid &ellipsoid, double standardParallel,
+	  double azimuth, double scale, double centerLatitude,
+	  double longitudeOrigin, double falseEasting, double falseNorthing)
+		: _k(ellipsoid, standardParallel, azimuth, scale, centerLatitude,
+		longitudeOrigin, falseEasting, falseNorthing) {}
+
+	virtual CT *clone() const {return new KrovakNE(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const
+	  {PointD p(_k.ll2xy(c)); return PointD(-p.x(), -p.y());}
+	virtual Coordinates xy2ll(const PointD &p) const
+	  {return _k.xy2ll(PointD(-p.x(), -p.y()));}
+
+private:
+	Krovak _k;
+};
+
+#endif // KROVAK_H

src/map/proj/lambertazimuthal.cpp

@@ -0,0 +1,78 @@
+#include <cmath>
+#include "map/ellipsoid.h"
+#include "lambertazimuthal.h"
+
+#define sin2(x) (sin(x) * sin(x))
+#define sqr(x) ((x) * (x))
+
+LambertAzimuthal::LambertAzimuthal(const Ellipsoid &ellipsoid,
+  double latitudeOrigin, double longitudeOrigin, double falseEasting,
+  double falseNorthing)
+{
+	double lat0 = deg2rad(latitudeOrigin);
+
+	_fe = falseEasting;
+	_fn = falseNorthing;
+	_lon0 = deg2rad(longitudeOrigin);
+
+	_a = ellipsoid.radius();
+	_es = ellipsoid.es();
+	_e = sqrt(_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 - _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 - _es * sin2(lat0))) / (_rq * cos(_beta0));
+}
+
+PointD LambertAzimuthal::ll2xy(const Coordinates &c) const
+{
+	double lon = deg2rad(c.lon());
+	double lat = deg2rad(c.lat());
+
+	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)
+	  * cos(beta) * cos(lon - _lon0))));
+
+	double x = _fe + ((B * _d) * (cos(beta) * sin(lon - _lon0)));
+	double y = _fn + (B / _d) * ((cos(_beta0) * sin(beta))
+	  - (sin(_beta0) * cos(beta) * cos(lon - _lon0)));
+
+	return PointD(x, y);
+}
+
+Coordinates LambertAzimuthal::xy2ll(const PointD &p) const
+{
+	double es4 = _es * _es;
+	double es6 = _es * es4;
+
+	double rho = sqrt(sqr((p.x() - _fe) / _d) + sqr(_d * (p.y()
+	  - _fn)));
+	double C = 2.0 * asin(rho / (2.0*_rq));
+	double betaS = asin((cos(C) * sin(_beta0)) + ((_d * (p.y() -_fn)
+	  * sin(C) * cos(_beta0)) / rho));
+
+	double lon = _lon0 + atan((p.x() - _fe) * sin(C) / (_d * rho
+	  * cos(_beta0) * cos(C) - sqr(_d) * (p.y() - _fn) * sin(_beta0)
+	  * sin(C)));
+	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));
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}
+
+bool LambertAzimuthal::operator==(const CT &ct) const
+{
+	const LambertAzimuthal *other = dynamic_cast<const LambertAzimuthal*>(&ct);
+	return (other != 0 && _lon0 == other->_lon0 && _fn == other->_fn
+	  && _fe == other->_fe && _a == other->_a && _es == other->_es
+	  && _beta0 == other->_beta0 && _rq == other->_rq && _d == other->_d);
+}

src/map/proj/lambertazimuthal.h

@@ -0,0 +1,26 @@
+#ifndef LAMBERTAZIMUTHAL_H
+#define LAMBERTAZIMUTHAL_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class LambertAzimuthal : public CT
+{
+public:
+	LambertAzimuthal(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new LambertAzimuthal(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _lon0;
+	double _fn, _fe;
+	double _a, _e, _es, _qP, _beta0, _rq, _d;
+};
+
+#endif // LAMBERTAZIMUTHAL_H

src/map/proj/lambertconic.cpp

@@ -0,0 +1,282 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include <cmath>
+#include "map/ellipsoid.h"
+#include "lambertconic.h"
+
+
+#define LAMBERT_m(clat, essin) (clat / sqrt(1.0 - essin * essin))
+#define LAMBERT2_t(lat, essin, es_over_2) \
+	(tan(M_PI_4 - lat / 2) * pow((1.0 + essin) / (1.0 - essin), es_over_2))
+#define LAMBERT1_t(lat, essin, es_over_2) \
+	(tan(M_PI_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_over_2))
+
+
+LambertConic1::LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
+  double longitudeOrigin, double scale, double falseEasting,
+  double falseNorthing)
+{
+	double e_sin;
+	double m0;
+	double lat_orig;
+
+
+	lat_orig = deg2rad(latitudeOrigin);
+	_longitudeOrigin = deg2rad(longitudeOrigin);
+	if (_longitudeOrigin > M_PI)
+		_longitudeOrigin -= 2 * M_PI;
+
+	_falseEasting = falseEasting;
+	_falseNorthing = falseNorthing;
+
+	_e = sqrt(ellipsoid.es());
+	_e_over_2 = _e / 2.0;
+
+	_n = sin(lat_orig);
+
+	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;
+
+	_rho_olat = _rho0;
+}
+
+PointD LambertConic1::ll2xy(const Coordinates &c) const
+{
+	double t;
+	double rho;
+	double dlam;
+	double theta;
+	double lat = deg2rad(c.lat());
+
+
+	if (fabs(fabs(lat) - M_PI_2) > 1.0e-10) {
+		t = LAMBERT1_t(lat, _e * sin(lat), _e_over_2);
+		rho = _rho0 * pow(t / _t0, _n);
+	} else
+		rho = 0.0;
+
+	dlam = deg2rad(c.lon()) - _longitudeOrigin;
+
+	if (dlam > M_PI)
+		dlam -= 2 * M_PI;
+	if (dlam < -M_PI)
+		dlam += 2 * M_PI;
+
+	theta = _n * dlam;
+
+	return PointD(rho * sin(theta) + _falseEasting, _rho_olat - rho
+	  * cos(theta) + _falseNorthing);
+}
+
+Coordinates LambertConic1::xy2ll(const PointD &p) const
+{
+	double dx;
+	double dy;
+	double rho;
+	double rho_olat_minus_dy;
+	double t;
+	double PHI;
+	double es_sin;
+	double tempPHI = 0.0;
+	double theta = 0.0;
+	double tolerance = 4.85e-10;
+	int count = 30;
+	double lat, lon;
+
+
+	dy = p.y() - _falseNorthing;
+	dx = p.x() - _falseEasting;
+	rho_olat_minus_dy = _rho_olat - dy;
+	rho = sqrt(dx * dx + (rho_olat_minus_dy) * (rho_olat_minus_dy));
+
+	if (_n < 0.0) {
+		rho *= -1.0;
+		dx *= -1.0;
+		rho_olat_minus_dy *= -1.0;
+	}
+
+	if (rho != 0.0) {
+		theta = atan2(dx, rho_olat_minus_dy) / _n;
+		t = _t0 * pow(rho / _rho0, 1 / _n);
+		PHI = M_PI_2 - 2.0 * atan(t);
+		while (fabs(PHI - tempPHI) > tolerance && count) {
+			tempPHI = PHI;
+			es_sin = _e * sin(PHI);
+			PHI = M_PI_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin),
+			  _e_over_2));
+			count--;
+		}
+
+		if (!count)
+			return Coordinates();
+
+		lat = PHI;
+		lon = theta + _longitudeOrigin;
+
+		if (fabs(lat) < 2.0e-7)
+			lat = 0.0;
+		if (lat > M_PI_2)
+			lat = M_PI_2;
+		else if (lat < -M_PI_2)
+			lat = -M_PI_2;
+
+		if (lon > M_PI) {
+			if (lon - M_PI < 3.5e-6)
+				lon = M_PI;
+			else
+				lon -= 2 * M_PI;
+		}
+		if (lon < -M_PI) {
+			if (fabs(lon + M_PI) < 3.5e-6)
+				lon = -M_PI;
+			else
+				lon += 2 * M_PI;
+		}
+
+		if (fabs(lon) < 2.0e-7)
+			lon = 0.0;
+		if (lon > M_PI)
+			lon = M_PI;
+		else if (lon < -M_PI)
+			lon = -M_PI;
+	} else {
+		if (_n > 0.0)
+			lat = M_PI_2;
+		else
+			lat = -M_PI_2;
+		lon = _longitudeOrigin;
+	}
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}
+
+bool LambertConic1::operator==(const CT &ct) const
+{
+	const LambertConic1 *other = dynamic_cast<const LambertConic1*>(&ct);
+	return (other != 0 && _longitudeOrigin == other->_longitudeOrigin
+	  && _falseEasting == other->_falseEasting
+	  && _falseNorthing == other->_falseNorthing && _e == other->_e
+	  && _n == other->_n && _rho0 == other->_rho0);
+}
+
+
+LambertConic2::LambertConic2(const Ellipsoid &ellipsoid,
+  double standardParallel1, double standardParallel2, double latitudeOrigin,
+  double longitudeOrigin, double falseEasting, double falseNorthing)
+{
+	double e, e_over_2, e_sin;
+	double lat0;
+	double k0;
+	double t0;
+	double t1, t2;
+	double t_olat;
+	double m0;
+	double m1;
+	double m2;
+	double n;
+	double const_value;
+	double sp1, sp2;
+	double lat_orig;
+
+
+	lat_orig = deg2rad(latitudeOrigin);
+	sp1 = deg2rad(standardParallel1);
+	sp2 = deg2rad(standardParallel2);
+
+	if (fabs(sp1 - sp2) > 1.0e-10) {
+		e = sqrt(ellipsoid.es());
+		e_over_2 = e / 2.0;
+
+		e_sin = e * sin(lat_orig);
+		t_olat = LAMBERT2_t(lat_orig, e_sin, e_over_2);
+
+		e_sin = e * sin(sp1);
+		m1 = LAMBERT_m(cos(sp1), e_sin);
+		t1 = LAMBERT2_t(sp1, e_sin, e_over_2);
+
+		e_sin = e * sin(sp2);
+		m2 = LAMBERT_m(cos(sp2), e_sin);
+		t2 = LAMBERT2_t(sp2, e_sin, e_over_2);
+
+		n = log(m1 / m2) / log(t1 / t2);
+
+		lat0 = asin(n);
+
+		e_sin = e * sin(lat0);
+		m0 = LAMBERT_m(cos(lat0), e_sin);
+		t0 = LAMBERT2_t(lat0, e_sin, e_over_2);
+
+		k0 = (m1 / m0) * (pow(t0 / t1, n));
+
+		const_value = ((ellipsoid.radius() * m2) / (n * pow(t2, n)));
+
+		falseNorthing += (const_value * pow(t_olat, n)) - (const_value
+		  * pow(t0, n));
+	} else {
+		lat0 = sp1;
+		k0 = 1.0;
+	}
+
+	_lc1 = LambertConic1(ellipsoid, rad2deg(lat0), longitudeOrigin, k0,
+	  falseEasting, falseNorthing);
+}
+
+PointD LambertConic2::ll2xy(const Coordinates &c) const
+{
+	return _lc1.ll2xy(c);
+}
+
+Coordinates LambertConic2::xy2ll(const PointD &p) const
+{
+	return _lc1.xy2ll(p);
+}
+
+bool LambertConic2::operator==(const CT &ct) const
+{
+	const LambertConic2 *other = dynamic_cast<const LambertConic2*>(&ct);
+	return (other != 0 && _lc1 == other->_lc1);
+}

src/map/proj/lambertconic.h

@@ -0,0 +1,52 @@
+#ifndef LAMBERTCONIC_H
+#define LAMBERTCONIC_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class LambertConic1 : public CT
+{
+public:
+	LambertConic1() {}
+	LambertConic1(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double scale, double falseEasting,
+	  double falseNorthing);
+
+	virtual CT *clone() const {return new LambertConic1(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _longitudeOrigin;
+	double _falseEasting;
+	double _falseNorthing;
+
+	double _e;
+	double _e_over_2;
+	double _n;
+	double _t0;
+	double _rho0;
+	double _rho_olat;
+};
+
+class LambertConic2 : public CT
+{
+public:
+	LambertConic2(const Ellipsoid &ellipsoid, double standardParallel1,
+	  double standardParallel2, double latitudeOrigin, double longitudeOrigin,
+	   double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new LambertConic2(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	LambertConic1 _lc1;
+};
+
+#endif // LAMBERTCONIC_H

src/map/proj/latlon.h

@@ -0,0 +1,28 @@
+#ifndef LATLON_H
+#define LATLON_H
+
+#include "map/ct.h"
+#include "map/angularunits.h"
+
+class LatLon : public CT
+{
+public:
+	LatLon(const AngularUnits &au) : _au(au) {}
+
+	virtual CT *clone() const {return new LatLon(*this);}
+	virtual bool operator==(const CT &ct) const
+	{
+		const LatLon *other = dynamic_cast<const LatLon*>(&ct);
+		return (other != 0 && _au == other->_au);
+	}
+
+	virtual PointD ll2xy(const Coordinates &c) const
+	  {return PointD(_au.fromDegrees(c.lon()), _au.fromDegrees(c.lat()));}
+	virtual Coordinates xy2ll(const PointD &p) const
+	  {return Coordinates(_au.toDegrees(p.x()), _au.toDegrees(p.y()));}
+
+private:
+	AngularUnits _au;
+};
+
+#endif // LATLON_H

src/map/proj/mercator.cpp

@@ -0,0 +1,133 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include "map/ellipsoid.h"
+#include "mercator.h"
+
+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;
+
+	_latitudeOrigin = deg2rad(latitudeOrigin);
+	_longitudeOrigin = deg2rad(longitudeOrigin);
+	if (_longitudeOrigin > M_PI)
+		_longitudeOrigin -= 2 * M_PI;
+	_falseNorthing = falseNorthing;
+	_falseEasting = falseEasting;
+
+	_a = ellipsoid.radius();
+	_e = sqrt(es);
+
+	sin_olat = sin(_latitudeOrigin);
+	_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;
+	_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;
+}
+
+PointD Mercator::ll2xy(const Coordinates &c) const
+{
+	double lon = deg2rad(c.lon());
+	double lat = deg2rad(c.lat());
+	double ctanz2;
+	double e_x_sinlat;
+	double delta_lon;
+	double tan_temp;
+	double pow_temp;
+
+	if (lon > M_PI)
+		lon -= 2 * M_PI;
+	e_x_sinlat = _e * sin(lat);
+	tan_temp = tan(M_PI_4 + lat / 2.e0);
+	pow_temp = pow((1.e0 - e_x_sinlat) / (1.e0 + e_x_sinlat), _e / 2.e0);
+	ctanz2 = tan_temp * pow_temp;
+	delta_lon = lon - _longitudeOrigin;
+	if (delta_lon > M_PI)
+	  delta_lon -= 2 * M_PI;
+	if (delta_lon < -M_PI)
+	  delta_lon += 2 * M_PI;
+
+	return PointD(_scaleFactor * _a * delta_lon + _falseEasting,
+	  _scaleFactor * _a * log(ctanz2) + _falseNorthing);
+}
+
+Coordinates Mercator::xy2ll(const PointD &p) const
+{
+	double dx;
+	double dy;
+	double xphi;
+	double lat, lon;
+
+	dy = p.y() - _falseNorthing;
+	dx = p.x() - _falseEasting;
+	lon = _longitudeOrigin + dx / (_scaleFactor * _a);
+	xphi = M_PI_2 - 2.e0 * atan(1.e0 / exp(dy / (_scaleFactor * _a)));
+	lat = xphi + _ab * sin(2.e0 * xphi) + _bb * sin(4.e0 * xphi)
+	  + _cb * sin(6.e0 * xphi) + _db * sin(8.e0 * xphi);
+	if (lon > M_PI)
+		lon -= 2 * M_PI;
+	if (lon < -M_PI)
+		lon += 2 * M_PI;
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}
+
+bool Mercator::operator==(const CT &ct) const
+{
+	const Mercator *other = dynamic_cast<const Mercator*>(&ct);
+	return (other != 0 && _a == other->_a && _e == other->_e
+	  && _latitudeOrigin == other->_latitudeOrigin
+	  && _longitudeOrigin == other->_longitudeOrigin
+	  && _falseNorthing == other->_falseNorthing
+	  && _falseEasting == other->_falseEasting);
+}

src/map/proj/mercator.h

@@ -0,0 +1,30 @@
+#ifndef MERCATOR_H
+#define MERCATOR_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class Mercator : public CT
+{
+public:
+	Mercator(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new Mercator(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _a, _e;
+	double _latitudeOrigin;
+	double _longitudeOrigin;
+	double _falseNorthing;
+	double _falseEasting;
+	double _scaleFactor;
+	double _ab, _bb, _cb, _db;
+};
+
+#endif // MERCATOR_H

src/map/proj/obliquestereographic.cpp

@@ -0,0 +1,92 @@
+#include "map/ellipsoid.h"
+#include "obliquestereographic.h"
+
+
+#define S1(x) ((1.0 + x) / (1.0 - x))
+#define S2(x) ((1.0 - _e * x) / (1.0 + _e * x))
+
+ObliqueStereographic::ObliqueStereographic(const Ellipsoid &ellipsoid,
+  double latitudeOrigin, double longitudeOrigin, double scale,
+  double falseEasting, double falseNorthing)
+  : _fe(falseEasting), _fn(falseNorthing)
+{
+	double lat0 = deg2rad(latitudeOrigin);
+	double sinPhi0 = sin(lat0);
+	double cosPhi0 = cos(lat0);
+	double sinLat0s = sinPhi0 * sinPhi0;
+
+	double a = ellipsoid.radius();
+	_es = ellipsoid.es();
+	_e = sqrt(_es);
+
+	double rho0 = a * (1.0 - _es) / pow(1.0 - _es * sinLat0s, 1.5);
+	double nu0 = a / sqrt(1.0 - _es * sinLat0s);
+
+	_n = sqrt(1.0 + ((_es * pow(cosPhi0, 4)) / (1.0 - _es)));
+	double w1 = pow(S1(sinPhi0) * pow(S2(sinPhi0), _e), _n);
+	double sinChi0 = (w1 - 1) / (w1 + 1);
+	_c = (_n + sinPhi0) * (1.0 - sinChi0) / ((_n - sinPhi0) * (1.0 + sinChi0));
+	double w2 = _c * w1;
+
+	_chi0 = asin((w2 - 1.0)/(w2 + 1.0));
+	_sinChi0 = sin(_chi0); _cosChi0 = cos(_chi0);
+	_lambda0 = deg2rad(longitudeOrigin);
+	_twoRk0 = 2.0 * sqrt(rho0 * nu0) * scale;
+	_g = _twoRk0 * tan(M_PI_4 - _chi0 / 2.0);
+	_h = 2.0 * _twoRk0 * tan(_chi0) + _g;
+}
+
+PointD ObliqueStereographic::ll2xy(const Coordinates &c) const
+{
+	double lon = deg2rad(c.lon());
+	double lat = deg2rad(c.lat());
+
+	double sinPhi = sin(lat);
+	double w = _c * pow(S1(sinPhi) * pow(S2(sinPhi), _e), _n);
+	double chi = asin((w - 1.0) / (w + 1.0));
+	double lambda = _n * (lon - _lambda0) + _lambda0;
+	double B = (1.0 + sin(chi) * _sinChi0 + cos(chi) * _cosChi0
+	  * cos(lambda - _lambda0));
+
+	return PointD(_fe + _twoRk0 * cos(chi) * sin(lambda - _lambda0) / B,
+	  _fn + _twoRk0 * (sin(chi) * _cosChi0 - cos(chi) * _sinChi0
+	  * cos(lambda - _lambda0)) / B);
+}
+
+Coordinates ObliqueStereographic::xy2ll(const PointD &p) const
+{
+	double i = atan((p.x() - _fe) / (_h + (p.y() - _fn)));
+	double j = atan((p.x() - _fe) / (_g - (p.y() - _fn))) - i;
+
+	double chi = _chi0 + 2.0 * atan(((p.y() - _fn) - (p.x() - _fe) * tan(j/2.0))
+	  / _twoRk0);
+	double lambda = j + 2.0 * i + _lambda0;
+
+	double psi0 = 0.5 * log((1.0 + sin(chi)) / (_c * (1.0 - sin(chi)))) / _n;
+	double phi1 = 2.0 * atan(pow(M_E, psi0)) - M_PI_2;
+
+	double psi, phi = phi1, prev = phi;
+	for (int i = 0; i < 8; i++) {
+		double sinPhi = sin(phi);
+		psi = log((tan(phi/2.0 + M_PI_4)) * pow((1.0 - _e * sinPhi)
+		  / (1.0 + _e * sinPhi), _e/2.0));
+		phi = phi - (psi - psi0) * cos(phi) * (1.0 - _es * sinPhi * sinPhi)
+		  / (1.0 - _es);
+		if (fabs(phi - prev) < 1e-10)
+			break;
+		prev = phi;
+	}
+
+	return Coordinates(rad2deg((lambda - _lambda0) / _n + _lambda0),
+	  rad2deg(phi));
+}
+
+bool ObliqueStereographic::operator==(const CT &ct) const
+{
+	const ObliqueStereographic *other
+	  = dynamic_cast<const ObliqueStereographic*>(&ct);
+	return (other != 0 && _es == other->_es && _chi0 == other->_chi0
+	  && _sinChi0 == other->_sinChi0 && _cosChi0 == other->_cosChi0
+	  && _lambda0 == other->_lambda0 && _n == other->_n && _c == other->_c
+	  && _fe == other->_fe && _fn == other->_fn && _twoRk0 == other->_twoRk0);
+}

src/map/proj/obliquestereographic.h

@@ -0,0 +1,31 @@
+#ifndef OBLIQUESTEREOGRAPHIC_H
+#define OBLIQUESTEREOGRAPHIC_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class ObliqueStereographic : public CT
+{
+public:
+	ObliqueStereographic(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double scale, double falseEasting,
+	  double falseNorthing);
+
+	virtual CT *clone() const {return new ObliqueStereographic(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _e, _es;
+	double _chi0, _sinChi0, _cosChi0;
+	double _lambda0;
+	double _n;
+	double _c;
+	double _fe, _fn;
+	double _twoRk0,	_g, _h;
+};
+
+#endif // OBLIQUESTEREOGRAPHIC_H

src/map/proj/polarstereographic.cpp

@@ -0,0 +1,203 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include "map/ellipsoid.h"
+#include "polarstereographic.h"
+
+
+#define POLAR_POW(EsSin) pow((1.0 - EsSin) / (1.0 + EsSin), _es_OVER_2)
+
+PolarStereographic::PolarStereographic(const Ellipsoid &ellipsoid,
+  double latitudeOrigin, double longitudeOrigin,
+  double falseEasting, double falseNorthing)
+{
+	_two_a = 2.0 * ellipsoid.radius();
+
+	if (longitudeOrigin > M_PI)
+		longitudeOrigin -= 2*M_PI;
+
+	if (latitudeOrigin < 0) {
+		_southernHemisphere = 1;
+		_originLatitude = -deg2rad(latitudeOrigin);
+		_originLongitude = -deg2rad(longitudeOrigin);
+	} else {
+		_southernHemisphere = 0;
+		_originLatitude = deg2rad(latitudeOrigin);
+		_originLongitude = deg2rad(longitudeOrigin);
+	}
+	_falseEasting = falseEasting;
+	_falseNorthing = falseNorthing;
+
+	double es2 = ellipsoid.es();
+	_es = sqrt(es2);
+	_es_OVER_2 = _es / 2.0;
+
+	if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10) {
+		double slat = sin(_originLatitude);
+		double essin = _es * slat;
+		double pow_es = POLAR_POW(essin);
+		double clat = cos(_originLatitude);
+		_mc = clat / sqrt(1.0 - essin * essin);
+		_a_mc = ellipsoid.radius() * _mc;
+		_tc = tan(M_PI_4 - _originLatitude / 2.0) / pow_es;
+	} else {
+		double one_PLUS_es = 1.0 + _es;
+		double one_MINUS_es = 1.0 - _es;
+		_e4 = sqrt(pow(one_PLUS_es, one_PLUS_es) * pow(one_MINUS_es,
+		  one_MINUS_es));
+	}
+}
+
+PointD PolarStereographic::ll2xy(const Coordinates &c) const
+{
+	double Easting, Northing;
+	double Longitude = deg2rad(c.lon());
+	double Latitude = deg2rad(c.lat());
+
+
+	if (fabs(fabs(Latitude) - M_PI_2) < 1.0e-10) {
+		Easting = 0.0;
+		Northing = 0.0;
+	} else {
+		if (_southernHemisphere != 0) {
+			Longitude *= -1.0;
+			Latitude *= -1.0;
+		}
+
+		double dlam = Longitude - _originLongitude;
+		if (dlam > M_PI)
+			dlam -= 2*M_PI;
+		if (dlam < -M_PI)
+			dlam += 2*M_PI;
+
+		double slat = sin(Latitude);
+		double essin = _es * slat;
+		double pow_es = POLAR_POW(essin);
+		double t = tan(M_PI_4 - Latitude / 2.0) / pow_es;
+
+		double rho;
+		if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10)
+			rho = _a_mc * t / _tc;
+		else
+			rho = _two_a * t / _e4;
+
+		Easting = rho * sin(dlam) + _falseEasting;
+
+		if (_southernHemisphere != 0) {
+			Easting *= -1.0;
+			Northing = rho * cos(dlam) + _falseNorthing;
+		} else
+			Northing = -rho * cos(dlam) + _falseNorthing;
+	}
+
+	return PointD(Easting, Northing);
+}
+
+Coordinates PolarStereographic::xy2ll(const PointD &p) const
+{
+	double Latitude, Longitude;
+	double dy = p.y() - _falseNorthing;
+	double dx = p.x() - _falseEasting;
+
+	if ((dy == 0.0) && (dx == 0.0)) {
+		Latitude = M_PI_2;
+		Longitude = _originLongitude;
+	} else {
+		if (_southernHemisphere != 0) {
+			dy *= -1.0;
+			dx *= -1.0;
+		}
+
+		double rho = sqrt(dx * dx + dy * dy);
+		double t;
+		if (fabs(fabs(_originLatitude) - M_PI_2) > 1.0e-10)
+			t = rho * _tc / (_a_mc);
+		else
+			t = rho * _e4 / (_two_a);
+		double PHI = M_PI_2 - 2.0 * atan(t);
+
+		double tempPHI = 0.0;
+		while (fabs(PHI - tempPHI) > 1.0e-10) {
+			tempPHI = PHI;
+			double sin_PHI = sin(PHI);
+			double essin =  _es * sin_PHI;
+			double pow_es = POLAR_POW(essin);
+			PHI = M_PI_2 - 2.0 * atan(t * pow_es);
+		}
+		Latitude = PHI;
+		Longitude = _originLongitude + atan2(dx, -dy);
+
+		if (Longitude > M_PI)
+			Longitude -= 2*M_PI;
+		else if (Longitude < -M_PI)
+			Longitude += 2*M_PI;
+
+		if (Latitude > M_PI_2)
+			Latitude = M_PI_2;
+		else if (Latitude < -M_PI_2)
+			Latitude = -M_PI_2;
+
+		if (Longitude > M_PI)
+			Longitude = M_PI;
+		else if (Longitude < -M_PI)
+			Longitude = -M_PI;
+	}
+
+	if (_southernHemisphere != 0) {
+		Latitude *= -1.0;
+		Longitude *= -1.0;
+	}
+
+	return Coordinates(rad2deg(Longitude), rad2deg(Latitude));
+}
+
+bool PolarStereographic::operator==(const CT &ct) const
+{
+	const PolarStereographic *other
+	  = dynamic_cast<const PolarStereographic*>(&ct);
+	return (other != 0 && _originLatitude == other->_originLatitude
+	  && _originLongitude == other->_originLongitude
+	  && _falseEasting == other->_falseEasting
+	  && _falseNorthing == other->_falseNorthing && _two_a == other->_two_a
+	  && _es == other->_es && _southernHemisphere == other->_southernHemisphere);
+}

src/map/proj/polarstereographic.h

@@ -0,0 +1,36 @@
+#ifndef POLARSTEREOGRAPHIC_H
+#define POLARSTEREOGRAPHIC_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class PolarStereographic  : public CT
+{
+public:
+	PolarStereographic(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new PolarStereographic(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _originLatitude;
+	double _originLongitude;
+	double _falseEasting;
+	double _falseNorthing;
+
+	double _a_mc;
+	double _es;
+	double _es_OVER_2;
+	double _two_a;
+	double _mc;
+	double _tc;
+	double _e4;
+	int _southernHemisphere;
+};
+
+#endif // POLARSTEREOGRAPHIC_H

src/map/proj/polyconic.cpp

@@ -0,0 +1,211 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include "map/ellipsoid.h"
+#include "polyconic.h"
+
+#define POLY_COEFF_TIMES_SIN(coeff, x, latit) \
+	(coeff * (sin (x * latit)))
+#define POLY_M(c0lat, c1s2lat, c2s4lat, c3s6lat) \
+	(_a * (c0lat - c1s2lat + c2s4lat - c3s6lat))
+#define FLOAT_EQ(x, v, epsilon) \
+	(((v - epsilon) < x) && (x < (v + epsilon)))
+
+Polyconic::Polyconic(const Ellipsoid &ellipsoid, double latitudeOrigin,
+  double longitudeOrigin, double falseEasting, double falseNorthing)
+{
+	double j, three_es4;
+	double lat, sin2lat, sin4lat, sin6lat;
+	double a2;
+	double b2;
+
+	_longitudeOrigin = deg2rad(longitudeOrigin);
+	_latitudeOrigin = deg2rad(latitudeOrigin);
+	_a = ellipsoid.radius();
+	_b = ellipsoid.b();
+	if (_longitudeOrigin > M_PI)
+		_longitudeOrigin -= 2 * M_PI;
+	_falseNorthing = falseNorthing;
+	_falseEasting = falseEasting;
+	a2 = _a * _a;
+	b2 = _b * _b;
+	_es2 = (a2 - b2) / a2;
+	_es4 = _es2 * _es2;
+	_es6 = _es4 * _es2;
+
+	j = 45.0 * _es6 / 1024.0;
+	three_es4 = 3.0 * _es4;
+	_c0 = 1.0 - _es2 / 4.0 - three_es4 / 64.0 - 5.0 * _es6 / 256.0;
+	_c1 = 3.0 * _es2 / 8.0 + three_es4 / 32.0 + j;
+	_c2 = 15.0 * _es4 / 256.0 + j;
+	_c3 = 35.0 * _es6 / 3072.0;
+
+	lat = _c0 * _latitudeOrigin;
+	sin2lat = POLY_COEFF_TIMES_SIN(_c1, 2.0, _latitudeOrigin);
+	sin4lat = POLY_COEFF_TIMES_SIN(_c2, 4.0, _latitudeOrigin);
+	sin6lat = POLY_COEFF_TIMES_SIN(_c3, 6.0, _latitudeOrigin);
+	_M0 = POLY_M(lat, sin2lat, sin4lat, sin6lat);
+}
+
+PointD Polyconic::ll2xy(const Coordinates &c) const
+{
+	double Longitude = deg2rad(c.lon());
+	double Latitude = deg2rad(c.lat());
+	double slat = sin(Latitude);
+	double lat, sin2lat, sin4lat, sin6lat;
+	double dlam;
+	double NN;
+	double NN_OVER_tlat;
+	double MM;
+	double EE;
+
+	dlam = Longitude - _longitudeOrigin;
+	if (dlam > M_PI)
+		dlam -= 2 * M_PI;
+	if (dlam < -M_PI)
+		dlam += 2 * M_PI;
+
+	if (Latitude == 0.0) {
+		return PointD(_a * dlam + _falseEasting,
+		  -_M0 + _falseNorthing);
+	} else {
+		NN = _a / sqrt(1.0 - _es2 * (slat * slat));
+		NN_OVER_tlat = NN  / tan(Latitude);
+		lat = _c0 * Latitude;
+		sin2lat = POLY_COEFF_TIMES_SIN(_c1, 2.0, Latitude);
+		sin4lat = POLY_COEFF_TIMES_SIN(_c2, 4.0, Latitude);
+		sin6lat = POLY_COEFF_TIMES_SIN(_c3, 6.0, Latitude);
+		MM = POLY_M(lat, sin2lat, sin4lat, sin6lat);
+		EE = dlam * slat;
+		return PointD(NN_OVER_tlat * sin(EE) + _falseEasting,
+		  MM - _M0 + NN_OVER_tlat * (1.0 - cos(EE)) + _falseNorthing);
+	}
+}
+
+Coordinates Polyconic::xy2ll(const PointD &p) const
+{
+	double dx;
+	double dy;
+	double dx_OVER_Poly_a;
+	double AA;
+	double BB;
+	double CC = 0.0;
+	double PHIn, Delta_PHI = 1.0;
+	double sin_PHIn;
+	double PHI, sin2PHI, sin4PHI, sin6PHI;
+	double Mn, Mn_prime, Ma;
+	double AA_Ma;
+	double Ma2_PLUS_BB;
+	double AA_MINUS_Ma;
+	double tolerance = 1.0e-12;
+	double Latitude;
+	double Longitude;
+
+
+	dy = p.y() - _falseNorthing;
+	dx = p.x() - _falseEasting;
+	dx_OVER_Poly_a = dx / _a;
+
+	if (FLOAT_EQ(dy,-_M0,1)) {
+		Latitude = 0.0;
+		Longitude = dx_OVER_Poly_a + _longitudeOrigin;
+	} else {
+		AA = (_M0 + dy) / _a;
+		BB = dx_OVER_Poly_a * dx_OVER_Poly_a + (AA * AA);
+		PHIn = AA;
+
+		while (fabs(Delta_PHI) > tolerance) {
+			sin_PHIn = sin(PHIn);
+			CC = sqrt(1.0 - _es2 * sin_PHIn * sin_PHIn) * tan(PHIn);
+			PHI = _c0 * PHIn;
+			sin2PHI = POLY_COEFF_TIMES_SIN(_c1, 2.0, PHIn);
+			sin4PHI = POLY_COEFF_TIMES_SIN(_c2, 4.0, PHIn);
+			sin6PHI = POLY_COEFF_TIMES_SIN(_c3, 6.0, PHIn);
+			Mn = POLY_M(PHI, sin2PHI, sin4PHI, sin6PHI);
+			Mn_prime = _c0 - 2.0 * _c1 * cos(2.0 * PHIn) + 4.0 * _c2
+			  * cos(4.0 * PHIn) - 6.0 * _c3 * cos(6.0 * PHIn);
+			Ma = Mn / _a;
+			AA_Ma = AA * Ma;
+			Ma2_PLUS_BB = Ma * Ma + BB;
+			AA_MINUS_Ma = AA - Ma;
+			Delta_PHI = (AA_Ma * CC + AA_MINUS_Ma - 0.5 * (Ma2_PLUS_BB) * CC) /
+			  (_es2 * sin2PHI * (Ma2_PLUS_BB - 2.0 * AA_Ma) / 4.0 * CC
+			   + (AA_MINUS_Ma) * (CC * Mn_prime - 2.0 / sin2PHI) - Mn_prime);
+			PHIn -= Delta_PHI;
+		}
+		Latitude = PHIn;
+
+		if (Latitude > M_PI_2)
+			Latitude = M_PI_2;
+		else if (Latitude < -M_PI_2)
+			Latitude = -M_PI_2;
+
+		if (FLOAT_EQ(fabs(Latitude), M_PI_2, 0.00001) || (Latitude == 0))
+			Longitude = _longitudeOrigin;
+		else
+			Longitude = (asin(dx_OVER_Poly_a * CC)) / sin(Latitude)
+			  + _longitudeOrigin;
+	}
+
+	if (Longitude > M_PI)
+		Longitude -= 2 * M_PI;
+	if (Longitude < -M_PI)
+		Longitude += 2 *M_PI;
+
+	if (Longitude > M_PI)
+		Longitude = M_PI;
+	else if (Longitude < -M_PI)
+		Longitude = -M_PI;
+
+	return Coordinates(rad2deg(Longitude), rad2deg(Latitude));
+}
+
+bool Polyconic::operator==(const CT &ct) const
+{
+	const Polyconic *other = dynamic_cast<const Polyconic*>(&ct);
+	return (other != 0 && _a == other->_a && _b == other->_b
+	  && _latitudeOrigin == other->_latitudeOrigin
+	  && _longitudeOrigin == other->_longitudeOrigin
+	  && _falseNorthing == other->_falseNorthing
+	  && _falseEasting == other->_falseEasting);
+}

src/map/proj/polyconic.h

@@ -0,0 +1,37 @@
+#ifndef POLYCONIC_H
+#define POLYCONIC_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class Polyconic : public CT
+{
+public:
+	Polyconic(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double falseEasting, double falseNorthing);
+
+	virtual CT *clone() const {return new Polyconic(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _a;
+	double _b;
+	double _es2;
+	double _es4;
+	double _es6;
+	double _M0;
+	double _c0;
+	double _c1;
+	double _c2;
+	double _c3;
+	double _longitudeOrigin;
+	double _latitudeOrigin;
+	double _falseEasting;
+	double _falseNorthing;
+};
+
+#endif // POLYCONIC_H

src/map/proj/transversemercator.cpp

@@ -0,0 +1,265 @@
+/*
+ * Based on libgeotrans with the following Source Code Disclaimer:
+
+1. The GEOTRANS source code ("the software") is provided free of charge by
+the National Imagery and Mapping Agency (NIMA) of the United States
+Department of Defense. Although NIMA makes no copyright claim under Title 17
+U.S.C., NIMA claims copyrights in the source code under other legal regimes.
+NIMA hereby grants to each user of the software a license to use and
+distribute the software, and develop derivative works.
+
+2. Warranty Disclaimer: The software was developed to meet only the internal
+requirements of the U.S. National Imagery and Mapping Agency. The software
+is provided "as is," and no warranty, express or implied, including but not
+limited to the implied warranties of merchantability and fitness for
+particular purpose or arising by statute or otherwise in law or from a
+course of dealing or usage in trade, is made by NIMA as to the accuracy and
+functioning of the software.
+
+3. NIMA and its personnel are not required to provide technical support or
+general assistance with respect to the software.
+
+4. Neither NIMA nor its personnel will be liable for any claims, losses, or
+damages arising from or connected with the use of the software. The user
+agrees to hold harmless the United States National Imagery and Mapping
+Agency. The user's sole and exclusive remedy is to stop using the software.
+
+5. NIMA requests that products developed using the software credit the
+source of the software with the following statement, "The product was
+developed using GEOTRANS, a product of the National Imagery and Mapping
+Agency and U.S. Army Engineering Research and Development Center."
+
+6. For any products developed using the software, NIMA requires a disclaimer
+that use of the software does not indicate endorsement or approval of the
+product by the Secretary of Defense or the National Imagery and Mapping
+Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
+the National Imagery and Mapping Agency, the initials "NIMA", the seal of
+the National Imagery and Mapping Agency, or any colorable imitation thereof
+shall not be used to imply approval, endorsement, or authorization of a
+product without prior written permission from United States Secretary of
+Defense.
+
+*/
+
+#include <cmath>
+#include "map/ellipsoid.h"
+#include "transversemercator.h"
+
+
+#define SPHSN(lat) \
+	((double)(_a / sqrt(1.e0 - _es * pow(sin(lat), 2))))
+#define SPHTMD(lat) \
+	((double)(_ap * lat - _bp * sin(2.e0 * lat) + _cp * sin(4.e0 * lat) \
+	  - _dp * sin(6.e0 * lat) + _ep * sin(8.e0 * lat)))
+#define DENOM(lat) \
+	((double)(sqrt(1.e0 - _es * pow(sin(lat),2))))
+#define SPHSR(lat) \
+	((double)(_a * (1.e0 - _es) / pow(DENOM(lat), 3)))
+
+
+TransverseMercator::TransverseMercator(const Ellipsoid &ellipsoid,
+  double latitudeOrigin, double longitudeOrigin, double scale,
+  double falseEasting, double falseNorthing)
+{
+	double tn, tn2, tn3, tn4, tn5;
+	double b;
+
+
+	_a = ellipsoid.radius();
+	_longitudeOrigin = deg2rad(longitudeOrigin);
+	_latitudeOrigin = deg2rad(latitudeOrigin);
+	_scale = scale;
+	_falseEasting = falseEasting;
+	_falseNorthing = falseNorthing;
+
+	_es = ellipsoid.es();
+	_ebs = (1 / (1 - _es)) - 1;
+	b = ellipsoid.b();
+
+	tn = (_a - b) / (_a + b);
+	tn2 = tn * tn;
+	tn3 = tn2 * tn;
+	tn4 = tn3 * tn;
+	tn5 = tn4 * tn;
+
+	_ap = _a * (1.e0 - tn + 5.e0 * (tn2 - tn3) / 4.e0 + 81.e0
+	  * (tn4 - tn5) / 64.e0);
+	_bp = 3.e0 * _a * (tn - tn2 + 7.e0 * (tn3 - tn4) / 8.e0 + 55.e0
+	  * tn5 / 64.e0 ) / 2.e0;
+	_cp = 15.e0 * _a * (tn2 - tn3 + 3.e0 * (tn4 - tn5 ) / 4.e0) / 16.0;
+	_dp = 35.e0 * _a * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
+	_ep = 315.e0 * _a * (tn4 - tn5) / 512.e0;
+}
+
+PointD TransverseMercator::ll2xy(const Coordinates &c) const
+{
+	double rl;
+	double cl, c2, c3, c5, c7;
+	double dlam;
+	double eta, eta2, eta3, eta4;
+	double sl, sn;
+	double t, tan2, tan3, tan4, tan5, tan6;
+	double t1, t2, t3, t4, t5, t6, t7, t8, t9;
+	double tmd, tmdo;
+	double x, y;
+
+
+	dlam = deg2rad(c.lon()) - _longitudeOrigin;
+
+	if (dlam > M_PI)
+		dlam -= 2 * M_PI;
+	if (dlam < -M_PI)
+		dlam += 2 * M_PI;
+	if (fabs(dlam) < 2.e-10)
+		dlam = 0.0;
+
+	rl = deg2rad(c.lat());
+	sl = sin(rl);
+	cl = cos(rl);
+	c2 = cl * cl;
+	c3 = c2 * cl;
+	c5 = c3 * c2;
+	c7 = c5 * c2;
+	t = sl / cl;
+	tan2 = t * t;
+	tan3 = tan2 * t;
+	tan4 = tan3 * t;
+	tan5 = tan4 * t;
+	tan6 = tan5 * t;
+	eta = _ebs * c2;
+	eta2 = eta * eta;
+	eta3 = eta2 * eta;
+	eta4 = eta3 * eta;
+
+	sn = SPHSN(rl);
+	tmd = SPHTMD(rl);
+	tmdo = SPHTMD (_latitudeOrigin);
+
+
+	t1 = (tmd - tmdo) * _scale;
+	t2 = sn * sl * cl * _scale / 2.e0;
+	t3 = sn * sl * c3 * _scale * (5.e0 - tan2 + 9.e0 * eta + 4.e0 * eta2)
+	  / 24.e0;
+	t4 = sn * sl * c5 * _scale * (61.e0 - 58.e0 * tan2 + tan4 + 270.e0 * eta
+	  - 330.e0 * tan2 * eta + 445.e0 * eta2 + 324.e0 * eta3 - 680.e0 * tan2
+	  * eta2 + 88.e0 * eta4 - 600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4)
+	  / 720.e0;
+	t5 = sn * sl * c7 * _scale * (1385.e0 - 3111.e0 * tan2 + 543.e0 * tan4
+	  - tan6) / 40320.e0;
+
+	y = _falseNorthing + t1 + pow(dlam, 2.e0) * t2 + pow(dlam, 4.e0) * t3
+	  + pow(dlam, 6.e0) * t4 + pow(dlam, 8.e0) * t5;
+
+
+	t6 = sn * cl * _scale;
+	t7 = sn * c3 * _scale * (1.e0 - tan2 + eta) /6.e0;
+	t8 = sn * c5 * _scale * (5.e0 - 18.e0 * tan2 + tan4 + 14.e0 * eta - 58.e0
+	  * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3 - 64.e0 * tan2 * eta2 - 24.e0
+	  * tan2 * eta3) / 120.e0;
+	t9 = sn * c7 * _scale * (61.e0 - 479.e0 * tan2 + 179.e0 * tan4 - tan6)
+	  / 5040.e0;
+
+	x = _falseEasting + dlam * t6 + pow(dlam, 3.e0) * t7 + pow(dlam, 5.e0)
+	  * t8 + pow(dlam, 7.e0) * t9;
+
+	return PointD(x, y);
+}
+
+Coordinates TransverseMercator::xy2ll(const PointD &p) const
+{
+	double cl;
+	double de;
+	double dlam;
+	double eta, eta2, eta3, eta4;
+	double ftphi;
+	double sn;
+	double sr;
+	double t, tan2, tan4;
+	double t10, t11, t12, t13, t14, t15, t16, t17;
+	double tmd, tmdo;
+	double lat, lon;
+
+
+	tmdo = SPHTMD(_latitudeOrigin);
+	tmd = tmdo + (p.y() - _falseNorthing) / _scale;
+
+	sr = SPHSR(0.e0);
+	ftphi = tmd / sr;
+
+	for (int i = 0; i < 5 ; i++) {
+		t10 = SPHTMD(ftphi);
+		sr = SPHSR(ftphi);
+		ftphi = ftphi + (tmd - t10) / sr;
+	}
+
+	sr = SPHSR(ftphi);
+	sn = SPHSN(ftphi);
+
+	cl = cos(ftphi);
+
+	t = tan(ftphi);
+	tan2 = t * t;
+	tan4 = tan2 * tan2;
+	eta = _ebs * pow(cl, 2);
+	eta2 = eta * eta;
+	eta3 = eta2 * eta;
+	eta4 = eta3 * eta;
+	de = p.x() - _falseEasting;
+	if (fabs(de) < 0.0001)
+		de = 0.0;
+
+	t10 = t / (2.e0 * sr * sn * pow(_scale, 2));
+	t11 = t * (5.e0  + 3.e0 * tan2 + eta - 4.e0 * pow(eta, 2) - 9.e0 * tan2
+	  * eta) / (24.e0 * sr * pow(sn, 3) * pow(_scale, 4));
+	t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4 - 252.e0 * tan2
+	  * eta - 3.e0 * eta2 + 100.e0 * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
+	  * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2 + 84.e0 * tan2 * eta3 - 192.e0
+	  * tan2 * eta4) / (720.e0 * sr * pow(sn, 5) * pow(_scale, 6));
+	t13 = t * (1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0 * pow(t,6))
+	  / (40320.e0 * sr * pow(sn, 7) * pow(_scale, 8));
+	lat = ftphi - pow(de, 2) * t10 + pow(de, 4) * t11 - pow(de, 6) * t12
+	  + pow(de, 8) * t13;
+
+	t14 = 1.e0 / (sn * cl * _scale);
+	t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn, 3) * cl * pow(_scale, 3));
+	t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2 + 8.e0 * tan2 * eta
+	  + 24.e0 * tan4 - 4.e0 * eta3 + 4.e0 * tan2 * eta2 + 24.e0 * tan2 * eta3)
+	  / (120.e0 * pow(sn, 5) * cl * pow(_scale, 5));
+	t17 = (61.e0 +  662.e0 * tan2 + 1320.e0 * tan4 + 720.e0 * pow(t,6))
+	  / (5040.e0 * pow(sn, 7) * cl * pow(_scale, 7));
+
+	dlam = de * t14 - pow(de, 3) * t15 + pow(de, 5) * t16 - pow(de, 7) * t17;
+
+	lon = _longitudeOrigin + dlam;
+	while (lat > deg2rad(90.0)) {
+		lat = M_PI - lat;
+		lon += M_PI;
+		if (lon > M_PI)
+			lon -= 2 * M_PI;
+	}
+
+	while (lat < deg2rad(-90.0)) {
+		lat = - (lat + M_PI);
+		lon += M_PI;
+		if (lon > M_PI)
+			lon -= 2 * M_PI;
+	}
+
+	if (lon > 2 * M_PI)
+		lon -= 2 * M_PI;
+	if (lon < -M_PI)
+		lon += 2 * M_PI;
+
+	return Coordinates(rad2deg(lon), rad2deg(lat));
+}
+
+bool TransverseMercator::operator==(const CT &ct) const
+{
+	const TransverseMercator *other
+	  = dynamic_cast<const TransverseMercator*>(&ct);
+	return (other != 0 && _longitudeOrigin == other->_longitudeOrigin
+	  && _latitudeOrigin == other->_latitudeOrigin && _scale == other->_scale
+	  && _falseEasting == other->_falseEasting
+	  && _falseNorthing == other->_falseNorthing && _a == other->_a
+	  && _es == other->_es);
+}

src/map/proj/transversemercator.h

@@ -0,0 +1,33 @@
+#ifndef TRANSVERSEMERCATOR_H
+#define TRANSVERSEMERCATOR_H
+
+#include "map/ct.h"
+
+class Ellipsoid;
+
+class TransverseMercator : public CT
+{
+public:
+	TransverseMercator(const Ellipsoid &ellipsoid, double latitudeOrigin,
+	  double longitudeOrigin, double scale, double falseEasting,
+	  double falseNorthing);
+
+	virtual CT *clone() const {return new TransverseMercator(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+
+private:
+	double _longitudeOrigin;
+	double _latitudeOrigin;
+	double _scale;
+	double _falseEasting;
+	double _falseNorthing;
+	double _a;
+	double _es;
+	double _ebs;
+	double _ap, _bp, _cp, _dp, _ep;
+};
+
+#endif // TRANSVERSEMERCATOR_H

src/map/proj/webmercator.cpp

@@ -0,0 +1,22 @@
+#include <cmath>
+#include "common/coordinates.h"
+#include "common/wgs84.h"
+#include "webmercator.h"
+
+PointD WebMercator::ll2xy(const Coordinates &c) const
+{
+	return PointD(deg2rad(c.lon()) * WGS84_RADIUS,
+	  log(tan(M_PI_4 + deg2rad(c.lat())/2.0)) * WGS84_RADIUS);
+}
+
+Coordinates WebMercator::xy2ll(const PointD &p) const
+{
+	return Coordinates(rad2deg(p.x() / WGS84_RADIUS),
+	  rad2deg(2.0 * atan(exp(p.y() / WGS84_RADIUS)) - M_PI_2));
+}
+
+bool WebMercator::operator==(const CT &ct) const
+{
+	const WebMercator *other = dynamic_cast<const WebMercator*>(&ct);
+	return (other != 0);
+}

src/map/proj/webmercator.h

@@ -0,0 +1,16 @@
+#ifndef WEBMERCATOR_H
+#define WEBMERCATOR_H
+
+#include "map/ct.h"
+
+class WebMercator : public CT
+{
+public:
+	virtual CT *clone() const {return new WebMercator(*this);}
+	virtual bool operator==(const CT &ct) const;
+
+	virtual PointD ll2xy(const Coordinates &c) const;
+	virtual Coordinates xy2ll(const PointD &p) const;
+};
+
+#endif // WEBMERCATOR_H