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 "ellipsoid.h"
45 #include "albersequal.h"
46
47
48 #define ONE_MINUS_SQR(x) (1.0 - (x) * (x))
49 #define ALBERS_Q(slat, one_minus_sqr_e_sin, es_sin) \
50 (_one_minus_es * ((slat) / (one_minus_sqr_e_sin) - \
51 (1 / (_two_e)) * log((1 - (es_sin)) / (1 + (es_sin)))))
52 #define ALBERS_M(clat, one_minus_sqr_e_sin) \
53 ((clat) / sqrt(one_minus_sqr_e_sin))
54
55
AlbersEqual(const Ellipsoid * ellipsoid,double standardParallel1,double standardParallel2,double latitudeOrigin,double longitudeOrigin,double falseEasting,double falseNorthing)56 AlbersEqual::AlbersEqual(const Ellipsoid *ellipsoid, double standardParallel1,
57 double standardParallel2, double latitudeOrigin, double longitudeOrigin,
58 double falseEasting, double falseNorthing)
59 {
60 double sin_lat, sin_lat1, sin_lat2, cos_lat1, cos_lat2;
61 double m1, m2, sqr_m1, sqr_m2;
62 double q0, q1, q2;
63 double e_sin, e_sin1, e_sin2;
64 double one_minus_sqr_e_sin1, one_minus_sqr_e_sin2;
65 double nq0;
66 double sp1, sp2;
67
68
69 _latitudeOrigin = deg2rad(latitudeOrigin);
70 _longitudeOrigin = deg2rad(longitudeOrigin);
71 _falseEasting = falseEasting;
72 _falseNorthing = falseNorthing;
73
74 sp1 = deg2rad(standardParallel1);
75 sp2 = deg2rad(standardParallel2);
76
77 _a2 = ellipsoid->radius() * ellipsoid->radius();
78 _es = ellipsoid->es();
79 _e = sqrt(_es);
80 _one_minus_es = 1 - _es;
81 _two_e = 2 * _e;
82
83 sin_lat = sin(_latitudeOrigin);
84 e_sin = _e * sin_lat;
85 q0 = ALBERS_Q(sin_lat, ONE_MINUS_SQR(e_sin), e_sin);
86
87 sin_lat1 = sin(sp1);
88 cos_lat1 = cos(sp1);
89 e_sin1 = _e * sin_lat1;
90 one_minus_sqr_e_sin1 = ONE_MINUS_SQR(e_sin1);
91 m1 = ALBERS_M(cos_lat1, one_minus_sqr_e_sin1);
92 q1 = ALBERS_Q(sin_lat1, one_minus_sqr_e_sin1, e_sin1);
93
94 sqr_m1 = m1 * m1;
95 if (fabs(sp1 - sp2) > 1.0e-10) {
96 sin_lat2 = sin(sp2);
97 cos_lat2 = cos(sp2);
98 e_sin2 = _e * sin_lat2;
99 one_minus_sqr_e_sin2 = ONE_MINUS_SQR(e_sin2);
100 m2 = ALBERS_M(cos_lat2, one_minus_sqr_e_sin2);
101 q2 = ALBERS_Q(sin_lat2, one_minus_sqr_e_sin2, e_sin2);
102 sqr_m2 = m2 * m2;
103 _n = (sqr_m1 - sqr_m2) / (q2 - q1);
104 } else
105 _n = sin_lat1;
106
107 _c = sqr_m1 + _n * q1;
108 _a_over_n = ellipsoid->radius() / _n;
109 nq0 = _n * q0;
110 _rho0 = (_c < nq0) ? 0 : _a_over_n * sqrt(_c - nq0);
111 }
112
ll2xy(const Coordinates & c) const113 PointD AlbersEqual::ll2xy(const Coordinates &c) const
114 {
115 double dlam;
116 double sin_lat;
117 double e_sin;
118 double q;
119 double rho;
120 double theta;
121 double nq;
122
123
124 dlam = deg2rad(c.lon()) - _longitudeOrigin;
125 if (dlam > M_PI)
126 dlam -= 2 * M_PI;
127 if (dlam < -M_PI)
128 dlam += 2 * M_PI;
129
130 sin_lat = sin(deg2rad(c.lat()));
131 e_sin = _e * sin_lat;
132 q = ALBERS_Q(sin_lat, ONE_MINUS_SQR(e_sin), e_sin);
133 nq = _n * q;
134 rho = (_c < nq) ? 0 : _a_over_n * sqrt(_c - nq);
135 theta = _n * dlam;
136
137 return PointD(rho * sin(theta) + _falseEasting,
138 _rho0 - rho * cos(theta) + _falseNorthing);
139 }
140
xy2ll(const PointD & p) const141 Coordinates AlbersEqual::xy2ll(const PointD &p) const
142 {
143 double dy, dx;
144 double rho0_minus_dy;
145 double q, qc, q_over_2;
146 double rho, rho_n;
147 double phi, delta_phi = 1.0;
148 double sin_phi;
149 double e_sin, one_minus_sqr_e_sin;
150 double theta = 0.0;
151 int count = 30;
152 double tolerance = 4.85e-10;
153 double lat, lon;
154
155
156 dy = p.y() - _falseNorthing;
157 dx = p.x() - _falseEasting;
158
159 rho0_minus_dy = _rho0 - dy;
160 rho = sqrt(dx * dx + rho0_minus_dy * rho0_minus_dy);
161
162 if (_n < 0) {
163 rho *= -1.0;
164 dx *= -1.0;
165 rho0_minus_dy *= -1.0;
166 }
167
168 if (rho != 0.0)
169 theta = atan2(dx, rho0_minus_dy);
170 rho_n = rho * _n;
171 q = (_c - (rho_n * rho_n) / _a2) / _n;
172 qc = 1 - ((_one_minus_es) / (_two_e)) * log((1.0 - _e) / (1.0 + _e));
173 if (fabs(fabs(qc) - fabs(q)) > 1.0e-6) {
174 q_over_2 = q / 2.0;
175 if (q_over_2 > 1.0)
176 lat = M_PI_2;
177 else if (q_over_2 < -1.0)
178 lat = -M_PI_2;
179 else {
180 phi = asin(q_over_2);
181 if (_e < 1.0e-10)
182 lat = phi;
183 else {
184 while ((fabs(delta_phi) > tolerance) && count) {
185 sin_phi = sin(phi);
186 e_sin = _e * sin_phi;
187 one_minus_sqr_e_sin = ONE_MINUS_SQR(e_sin);
188 delta_phi = (one_minus_sqr_e_sin * one_minus_sqr_e_sin)
189 / (2.0 * cos(phi)) * (q / (_one_minus_es) - sin_phi
190 / one_minus_sqr_e_sin + (log((1.0 - e_sin)
191 / (1.0 + e_sin)) / (_two_e)));
192 phi += delta_phi;
193 count --;
194 }
195
196 lat = phi;
197 }
198
199 if (lat > M_PI_2)
200 lat = M_PI_2;
201 else if (lat < -M_PI_2)
202 lat = -M_PI_2;
203 }
204 } else {
205 if (q >= 0.0)
206 lat = M_PI_2;
207 else
208 lat = -M_PI_2;
209 }
210
211 lon = _longitudeOrigin + theta / _n;
212
213 if (lon > M_PI)
214 lon -= 2 * M_PI;
215 if (lon < -M_PI)
216 lon += 2 * M_PI;
217
218 if (lon > M_PI)
219 lon = M_PI;
220 else if (lon < -M_PI)
221 lon = -M_PI;
222
223 return Coordinates(rad2deg(lon), rad2deg(lat));
224 }
225
operator ==(const CT & ct) const226 bool AlbersEqual::operator==(const CT &ct) const
227 {
228 const AlbersEqual *other = dynamic_cast<const AlbersEqual*>(&ct);
229 return (other != 0 && _latitudeOrigin == other->_latitudeOrigin
230 && _longitudeOrigin == other->_longitudeOrigin
231 && _falseEasting == other->_falseEasting
232 && _falseNorthing == other->_falseNorthing && _a2 == other->_a2
233 && _es == other->_es && _rho0 == other->_rho0 && _c == other->_c
234 && _n == other->_n);
235 }
236