2 * Based on libgeotrans with the following Source Code Disclaimer:
4 1. The GEOTRANS source code ("the software") is provided free of charge by
5 the National Imagery and Mapping Agency (NIMA) of the United States
6 Department of Defense. Although NIMA makes no copyright claim under Title 17
7 U.S.C., NIMA claims copyrights in the source code under other legal regimes.
8 NIMA hereby grants to each user of the software a license to use and
9 distribute the software, and develop derivative works.
11 2. Warranty Disclaimer: The software was developed to meet only the internal
12 requirements of the U.S. National Imagery and Mapping Agency. The software
13 is provided "as is," and no warranty, express or implied, including but not
14 limited to the implied warranties of merchantability and fitness for
15 particular purpose or arising by statute or otherwise in law or from a
16 course of dealing or usage in trade, is made by NIMA as to the accuracy and
17 functioning of the software.
19 3. NIMA and its personnel are not required to provide technical support or
20 general assistance with respect to the software.
22 4. Neither NIMA nor its personnel will be liable for any claims, losses, or
23 damages arising from or connected with the use of the software. The user
24 agrees to hold harmless the United States National Imagery and Mapping
25 Agency. The user's sole and exclusive remedy is to stop using the software.
27 5. NIMA requests that products developed using the software credit the
28 source of the software with the following statement, "The product was
29 developed using GEOTRANS, a product of the National Imagery and Mapping
30 Agency and U.S. Army Engineering Research and Development Center."
32 6. For any products developed using the software, NIMA requires a disclaimer
33 that use of the software does not indicate endorsement or approval of the
34 product by the Secretary of Defense or the National Imagery and Mapping
35 Agency. Pursuant to the United States Code, 10 U.S.C. Sec. 2797, the name of
36 the National Imagery and Mapping Agency, the initials "NIMA", the seal of
37 the National Imagery and Mapping Agency, or any colorable imitation thereof
38 shall not be used to imply approval, endorsement, or authorization of a
39 product without prior written permission from United States Secretary of
44 #include "map/ellipsoid.h"
45 #include "polarstereographic.h"
48 #define POLAR_POW(EsSin) pow((1.0 - EsSin) / (1.0 + EsSin), _es_OVER_2)
50 PolarStereographic::PolarStereographic(const Ellipsoid
&ellipsoid
,
51 double latitudeOrigin
, double longitudeOrigin
,
52 double falseEasting
, double falseNorthing
)
54 _two_a
= 2.0 * ellipsoid
.radius();
56 if (longitudeOrigin
> M_PI
)
57 longitudeOrigin
-= 2*M_PI
;
59 if (latitudeOrigin
< 0) {
60 _southernHemisphere
= 1;
61 _originLatitude
= -deg2rad(latitudeOrigin
);
62 _originLongitude
= -deg2rad(longitudeOrigin
);
64 _southernHemisphere
= 0;
65 _originLatitude
= deg2rad(latitudeOrigin
);
66 _originLongitude
= deg2rad(longitudeOrigin
);
68 _falseEasting
= falseEasting
;
69 _falseNorthing
= falseNorthing
;
71 double es2
= ellipsoid
.es();
73 _es_OVER_2
= _es
/ 2.0;
75 if (fabs(fabs(_originLatitude
) - M_PI_2
) > 1.0e-10) {
76 double slat
= sin(_originLatitude
);
77 double essin
= _es
* slat
;
78 double pow_es
= POLAR_POW(essin
);
79 double clat
= cos(_originLatitude
);
80 _mc
= clat
/ sqrt(1.0 - essin
* essin
);
81 _a_mc
= ellipsoid
.radius() * _mc
;
82 _tc
= tan(M_PI_4
- _originLatitude
/ 2.0) / pow_es
;
84 double one_PLUS_es
= 1.0 + _es
;
85 double one_MINUS_es
= 1.0 - _es
;
86 _e4
= sqrt(pow(one_PLUS_es
, one_PLUS_es
) * pow(one_MINUS_es
,
91 PointD
PolarStereographic::ll2xy(const Coordinates
&c
) const
93 double Easting
, Northing
;
94 double Longitude
= deg2rad(c
.lon());
95 double Latitude
= deg2rad(c
.lat());
98 if (fabs(fabs(Latitude
) - M_PI_2
) < 1.0e-10) {
102 if (_southernHemisphere
!= 0) {
107 double dlam
= Longitude
- _originLongitude
;
113 double slat
= sin(Latitude
);
114 double essin
= _es
* slat
;
115 double pow_es
= POLAR_POW(essin
);
116 double t
= tan(M_PI_4
- Latitude
/ 2.0) / pow_es
;
119 if (fabs(fabs(_originLatitude
) - M_PI_2
) > 1.0e-10)
120 rho
= _a_mc
* t
/ _tc
;
122 rho
= _two_a
* t
/ _e4
;
124 Easting
= rho
* sin(dlam
) + _falseEasting
;
126 if (_southernHemisphere
!= 0) {
128 Northing
= rho
* cos(dlam
) + _falseNorthing
;
130 Northing
= -rho
* cos(dlam
) + _falseNorthing
;
133 return PointD(Easting
, Northing
);
136 Coordinates
PolarStereographic::xy2ll(const PointD
&p
) const
138 double Latitude
, Longitude
;
139 double dy
= p
.y() - _falseNorthing
;
140 double dx
= p
.x() - _falseEasting
;
142 if ((dy
== 0.0) && (dx
== 0.0)) {
144 Longitude
= _originLongitude
;
146 if (_southernHemisphere
!= 0) {
151 double rho
= sqrt(dx
* dx
+ dy
* dy
);
153 if (fabs(fabs(_originLatitude
) - M_PI_2
) > 1.0e-10)
154 t
= rho
* _tc
/ (_a_mc
);
156 t
= rho
* _e4
/ (_two_a
);
157 double PHI
= M_PI_2
- 2.0 * atan(t
);
159 double tempPHI
= 0.0;
160 while (fabs(PHI
- tempPHI
) > 1.0e-10) {
162 double sin_PHI
= sin(PHI
);
163 double essin
= _es
* sin_PHI
;
164 double pow_es
= POLAR_POW(essin
);
165 PHI
= M_PI_2
- 2.0 * atan(t
* pow_es
);
168 Longitude
= _originLongitude
+ atan2(dx
, -dy
);
170 if (Longitude
> M_PI
)
172 else if (Longitude
< -M_PI
)
175 if (Latitude
> M_PI_2
)
177 else if (Latitude
< -M_PI_2
)
180 if (Longitude
> M_PI
)
182 else if (Longitude
< -M_PI
)
186 if (_southernHemisphere
!= 0) {
191 return Coordinates(rad2deg(Longitude
), rad2deg(Latitude
));
194 bool PolarStereographic::operator==(const CT
&ct
) const
196 const PolarStereographic
*other
197 = dynamic_cast<const PolarStereographic
*>(&ct
);
198 return (other
!= 0 && _originLatitude
== other
->_originLatitude
199 && _originLongitude
== other
->_originLongitude
200 && _falseEasting
== other
->_falseEasting
201 && _falseNorthing
== other
->_falseNorthing
&& _two_a
== other
->_two_a
202 && _es
== other
->_es
&& _southernHemisphere
== other
->_southernHemisphere
);