1 /*
2 * Based on libgeotrans with the following Source Code Disclaimer:
3
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.
10
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.
18
19 3. NIMA and its personnel are not required to provide technical support or
20 general assistance with respect to the software.
21
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.
26
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."
31
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
40 Defense.
41
42 */
43
44 #include <cmath>
45 #include "ellipsoid.h"
46 #include "lambertconic.h"
47
48
49 #define LAMBERT_m(clat, essin) (clat / sqrt(1.0 - essin * essin))
50 #define LAMBERT2_t(lat, essin, es_over_2) \
51 (tan(M_PI_4 - lat / 2) * pow((1.0 + essin) / (1.0 - essin), es_over_2))
52 #define LAMBERT1_t(lat, essin, es_over_2) \
53 (tan(M_PI_4 - lat / 2) / pow((1.0 - essin) / (1.0 + essin), es_over_2))
54
55
LambertConic1(const Ellipsoid * ellipsoid,double latitudeOrigin,double longitudeOrigin,double scale,double falseEasting,double falseNorthing)56 LambertConic1::LambertConic1(const Ellipsoid *ellipsoid, double latitudeOrigin,
57 double longitudeOrigin, double scale, double falseEasting,
58 double falseNorthing)
59 {
60 double e_sin;
61 double m0;
62 double lat_orig;
63
64
65 lat_orig = deg2rad(latitudeOrigin);
66 _longitudeOrigin = deg2rad(longitudeOrigin);
67 if (_longitudeOrigin > M_PI)
68 _longitudeOrigin -= 2 * M_PI;
69
70 _falseEasting = falseEasting;
71 _falseNorthing = falseNorthing;
72
73 _e = sqrt(ellipsoid->es());
74 _e_over_2 = _e / 2.0;
75
76 _n = sin(lat_orig);
77
78 e_sin = _e * sin(lat_orig);
79 m0 = LAMBERT_m(cos(lat_orig), e_sin);
80 _t0 = LAMBERT1_t(lat_orig, e_sin, _e_over_2);
81
82 _rho0 = ellipsoid->radius() * scale * m0 / _n;
83
84 _rho_olat = _rho0;
85 }
86
ll2xy(const Coordinates & c) const87 PointD LambertConic1::ll2xy(const Coordinates &c) const
88 {
89 double t;
90 double rho;
91 double dlam;
92 double theta;
93 double lat = deg2rad(c.lat());
94
95
96 if (fabs(fabs(lat) - M_PI_2) > 1.0e-10) {
97 t = LAMBERT1_t(lat, _e * sin(lat), _e_over_2);
98 rho = _rho0 * pow(t / _t0, _n);
99 } else
100 rho = 0.0;
101
102 dlam = deg2rad(c.lon()) - _longitudeOrigin;
103
104 if (dlam > M_PI)
105 dlam -= 2 * M_PI;
106 if (dlam < -M_PI)
107 dlam += 2 * M_PI;
108
109 theta = _n * dlam;
110
111 return PointD(rho * sin(theta) + _falseEasting, _rho_olat - rho
112 * cos(theta) + _falseNorthing);
113 }
114
xy2ll(const PointD & p) const115 Coordinates LambertConic1::xy2ll(const PointD &p) const
116 {
117 double dx;
118 double dy;
119 double rho;
120 double rho_olat_minus_dy;
121 double t;
122 double PHI;
123 double es_sin;
124 double tempPHI = 0.0;
125 double theta = 0.0;
126 double tolerance = 4.85e-10;
127 int count = 30;
128 double lat, lon;
129
130
131 dy = p.y() - _falseNorthing;
132 dx = p.x() - _falseEasting;
133 rho_olat_minus_dy = _rho_olat - dy;
134 rho = sqrt(dx * dx + (rho_olat_minus_dy) * (rho_olat_minus_dy));
135
136 if (_n < 0.0) {
137 rho *= -1.0;
138 dx *= -1.0;
139 rho_olat_minus_dy *= -1.0;
140 }
141
142 if (rho != 0.0) {
143 theta = atan2(dx, rho_olat_minus_dy) / _n;
144 t = _t0 * pow(rho / _rho0, 1 / _n);
145 PHI = M_PI_2 - 2.0 * atan(t);
146 while (fabs(PHI - tempPHI) > tolerance && count) {
147 tempPHI = PHI;
148 es_sin = _e * sin(PHI);
149 PHI = M_PI_2 - 2.0 * atan(t * pow((1.0 - es_sin) / (1.0 + es_sin),
150 _e_over_2));
151 count--;
152 }
153
154 if (!count)
155 return Coordinates();
156
157 lat = PHI;
158 lon = theta + _longitudeOrigin;
159
160 if (fabs(lat) < 2.0e-7)
161 lat = 0.0;
162 if (lat > M_PI_2)
163 lat = M_PI_2;
164 else if (lat < -M_PI_2)
165 lat = -M_PI_2;
166
167 if (lon > M_PI) {
168 if (lon - M_PI < 3.5e-6)
169 lon = M_PI;
170 else
171 lon -= 2 * M_PI;
172 }
173 if (lon < -M_PI) {
174 if (fabs(lon + M_PI) < 3.5e-6)
175 lon = -M_PI;
176 else
177 lon += 2 * M_PI;
178 }
179
180 if (fabs(lon) < 2.0e-7)
181 lon = 0.0;
182 if (lon > M_PI)
183 lon = M_PI;
184 else if (lon < -M_PI)
185 lon = -M_PI;
186 } else {
187 if (_n > 0.0)
188 lat = M_PI_2;
189 else
190 lat = -M_PI_2;
191 lon = _longitudeOrigin;
192 }
193
194 return Coordinates(rad2deg(lon), rad2deg(lat));
195 }
196
operator ==(const CT & ct) const197 bool LambertConic1::operator==(const CT &ct) const
198 {
199 const LambertConic1 *other = dynamic_cast<const LambertConic1*>(&ct);
200 return (other != 0 && _longitudeOrigin == other->_longitudeOrigin
201 && _falseEasting == other->_falseEasting
202 && _falseNorthing == other->_falseNorthing && _e == other->_e
203 && _n == other->_n && _rho0 == other->_rho0);
204 }
205
206
LambertConic2(const Ellipsoid * ellipsoid,double standardParallel1,double standardParallel2,double latitudeOrigin,double longitudeOrigin,double falseEasting,double falseNorthing)207 LambertConic2::LambertConic2(const Ellipsoid *ellipsoid,
208 double standardParallel1, double standardParallel2, double latitudeOrigin,
209 double longitudeOrigin, double falseEasting, double falseNorthing)
210 {
211 double e, e_over_2, e_sin;
212 double lat0;
213 double k0;
214 double t0;
215 double t1, t2;
216 double t_olat;
217 double m0;
218 double m1;
219 double m2;
220 double n;
221 double const_value;
222 double sp1, sp2;
223 double lat_orig;
224
225
226 lat_orig = deg2rad(latitudeOrigin);
227 sp1 = deg2rad(standardParallel1);
228 sp2 = deg2rad(standardParallel2);
229
230 if (fabs(sp1 - sp2) > 1.0e-10) {
231 e = sqrt(ellipsoid->es());
232 e_over_2 = e / 2.0;
233
234 e_sin = e * sin(lat_orig);
235 t_olat = LAMBERT2_t(lat_orig, e_sin, e_over_2);
236
237 e_sin = e * sin(sp1);
238 m1 = LAMBERT_m(cos(sp1), e_sin);
239 t1 = LAMBERT2_t(sp1, e_sin, e_over_2);
240
241 e_sin = e * sin(sp2);
242 m2 = LAMBERT_m(cos(sp2), e_sin);
243 t2 = LAMBERT2_t(sp2, e_sin, e_over_2);
244
245 n = log(m1 / m2) / log(t1 / t2);
246
247 lat0 = asin(n);
248
249 e_sin = e * sin(lat0);
250 m0 = LAMBERT_m(cos(lat0), e_sin);
251 t0 = LAMBERT2_t(lat0, e_sin, e_over_2);
252
253 k0 = (m1 / m0) * (pow(t0 / t1, n));
254
255 const_value = ((ellipsoid->radius() * m2) / (n * pow(t2, n)));
256
257 falseNorthing += (const_value * pow(t_olat, n)) - (const_value
258 * pow(t0, n));
259 } else {
260 lat0 = sp1;
261 k0 = 1.0;
262 }
263
264 _lc1 = LambertConic1(ellipsoid, rad2deg(lat0), longitudeOrigin, k0,
265 falseEasting, falseNorthing);
266 }
267
ll2xy(const Coordinates & c) const268 PointD LambertConic2::ll2xy(const Coordinates &c) const
269 {
270 return _lc1.ll2xy(c);
271 }
272
xy2ll(const PointD & p) const273 Coordinates LambertConic2::xy2ll(const PointD &p) const
274 {
275 return _lc1.xy2ll(p);
276 }
277
operator ==(const CT & ct) const278 bool LambertConic2::operator==(const CT &ct) const
279 {
280 const LambertConic2 *other = dynamic_cast<const LambertConic2*>(&ct);
281 return (other != 0 && _lc1 == other->_lc1);
282 }
283