GeographicLib 2.5
AlbersEqualArea.cpp
Go to the documentation of this file.
1/**
2 * \file AlbersEqualArea.cpp
3 * \brief Implementation for GeographicLib::AlbersEqualArea class
4 *
5 * Copyright (c) Charles Karney (2010-2022) <karney@alum.mit.edu> and licensed
6 * under the MIT/X11 License. For more information, see
7 * https://geographiclib.sourceforge.io/
8 **********************************************************************/
9
10#include <GeographicLib/AlbersEqualArea.hpp>
11
12namespace GeographicLib {
13
14 using namespace std;
15
16 AlbersEqualArea::AlbersEqualArea(real a, real f, real stdlat, real k0)
17 : eps_(numeric_limits<real>::epsilon())
18 , epsx_(Math::sq(eps_))
19 , epsx2_(Math::sq(epsx_))
20 , tol_(sqrt(eps_))
21 , tol0_(tol_ * sqrt(sqrt(eps_)))
22 , _a(a)
23 , _f(f)
24 , _fm(1 - _f)
25 , _e2(_f * (2 - _f))
26 , _e(sqrt(fabs(_e2)))
27 , _e2m(1 - _e2)
28 , _qZ(1 + _e2m * atanhee(real(1)))
29 , _qx(_qZ / ( 2 * _e2m ))
30 {
31 if (!(isfinite(_a) && _a > 0))
32 throw GeographicErr("Equatorial radius is not positive");
33 if (!(isfinite(_f) && _f < 1))
34 throw GeographicErr("Polar semi-axis is not positive");
35 if (!(isfinite(k0) && k0 > 0))
36 throw GeographicErr("Scale is not positive");
37 if (!(fabs(stdlat) <= Math::qd))
38 throw GeographicErr("Standard latitude not in [-" + to_string(Math::qd)
39 + "d, " + to_string(Math::qd) + "d]");
40 real sphi, cphi;
41 Math::sincosd(stdlat, sphi, cphi);
42 Init(sphi, cphi, sphi, cphi, k0);
43 }
44
45 AlbersEqualArea::AlbersEqualArea(real a, real f, real stdlat1, real stdlat2,
46 real k1)
47 : eps_(numeric_limits<real>::epsilon())
48 , epsx_(Math::sq(eps_))
49 , epsx2_(Math::sq(epsx_))
50 , tol_(sqrt(eps_))
51 , tol0_(tol_ * sqrt(sqrt(eps_)))
52 , _a(a)
53 , _f(f)
54 , _fm(1 - _f)
55 , _e2(_f * (2 - _f))
56 , _e(sqrt(fabs(_e2)))
57 , _e2m(1 - _e2)
58 , _qZ(1 + _e2m * atanhee(real(1)))
59 , _qx(_qZ / ( 2 * _e2m ))
60 {
61 if (!(isfinite(_a) && _a > 0))
62 throw GeographicErr("Equatorial radius is not positive");
63 if (!(isfinite(_f) && _f < 1))
64 throw GeographicErr("Polar semi-axis is not positive");
65 if (!(isfinite(k1) && k1 > 0))
66 throw GeographicErr("Scale is not positive");
67 if (!(fabs(stdlat1) <= Math::qd))
68 throw GeographicErr("Standard latitude 1 not in [-"
69 + to_string(Math::qd) + "d, "
70 + to_string(Math::qd) + "d]");
71 if (!(fabs(stdlat2) <= Math::qd))
72 throw GeographicErr("Standard latitude 2 not in [-"
73 + to_string(Math::qd) + "d, "
74 + to_string(Math::qd) + "d]");
75 real sphi1, cphi1, sphi2, cphi2;
76 Math::sincosd(stdlat1, sphi1, cphi1);
77 Math::sincosd(stdlat2, sphi2, cphi2);
78 Init(sphi1, cphi1, sphi2, cphi2, k1);
79 }
80
81 AlbersEqualArea::AlbersEqualArea(real a, real f,
82 real sinlat1, real coslat1,
83 real sinlat2, real coslat2,
84 real k1)
85 : eps_(numeric_limits<real>::epsilon())
86 , epsx_(Math::sq(eps_))
87 , epsx2_(Math::sq(epsx_))
88 , tol_(sqrt(eps_))
89 , tol0_(tol_ * sqrt(sqrt(eps_)))
90 , _a(a)
91 , _f(f)
92 , _fm(1 - _f)
93 , _e2(_f * (2 - _f))
94 , _e(sqrt(fabs(_e2)))
95 , _e2m(1 - _e2)
96 , _qZ(1 + _e2m * atanhee(real(1)))
97 , _qx(_qZ / ( 2 * _e2m ))
98 {
99 if (!(isfinite(_a) && _a > 0))
100 throw GeographicErr("Equatorial radius is not positive");
101 if (!(isfinite(_f) && _f < 1))
102 throw GeographicErr("Polar semi-axis is not positive");
103 if (!(isfinite(k1) && k1 > 0))
104 throw GeographicErr("Scale is not positive");
105 if (signbit(coslat1))
106 throw GeographicErr("Standard latitude 1 not in [-"
107 + to_string(Math::qd) + "d, "
108 + to_string(Math::qd) + "d]");
109 if (signbit(coslat2))
110 throw GeographicErr("Standard latitude 2 not in [-"
111 + to_string(Math::qd) + "d, "
112 + to_string(Math::qd) + "d]");
113 if (!(fabs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
114 throw GeographicErr("Bad sine/cosine of standard latitude 1");
115 if (!(fabs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
116 throw GeographicErr("Bad sine/cosine of standard latitude 2");
117 if (coslat1 == 0 && coslat2 == 0 && sinlat1 * sinlat2 <= 0)
118 throw GeographicErr
119 ("Standard latitudes cannot be opposite poles");
120 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
121 }
122
123 void AlbersEqualArea::Init(real sphi1, real cphi1,
124 real sphi2, real cphi2, real k1) {
125 {
126 real r;
127 r = hypot(sphi1, cphi1);
128 sphi1 /= r; cphi1 /= r;
129 r = hypot(sphi2, cphi2);
130 sphi2 /= r; cphi2 /= r;
131 }
132 bool polar = (cphi1 == 0);
133 cphi1 = fmax(epsx_, cphi1); // Avoid singularities at poles
134 cphi2 = fmax(epsx_, cphi2);
135 // Determine hemisphere of tangent latitude
136 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
137 // Internally work with tangent latitude positive
138 sphi1 *= _sign; sphi2 *= _sign;
139 if (sphi1 > sphi2) {
140 swap(sphi1, sphi2); swap(cphi1, cphi2); // Make phi1 < phi2
141 }
142 real
143 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2;
144
145 // q = (1-e^2)*(sphi/(1-e^2*sphi^2) - atanhee(sphi))
146 // qZ = q(pi/2) = (1 + (1-e^2)*atanhee(1))
147 // atanhee(x) = atanh(e*x)/e
148 // q = sxi * qZ
149 // dq/dphi = 2*(1-e^2)*cphi/(1-e^2*sphi^2)^2
150 //
151 // n = (m1^2-m2^2)/(q2-q1) -> sin(phi0) for phi1, phi2 -> phi0
152 // C = m1^2 + n*q1 = (m1^2*q2-m2^2*q1)/(q2-q1)
153 // let
154 // rho(pi/2)/rho(-pi/2) = (1-s)/(1+s)
155 // s = n*qZ/C
156 // = qZ * (m1^2-m2^2)/(m1^2*q2-m2^2*q1)
157 // = qZ * (scbet2^2 - scbet1^2)/(scbet2^2*q2 - scbet1^2*q1)
158 // = (scbet2^2 - scbet1^2)/(scbet2^2*sxi2 - scbet1^2*sxi1)
159 // = (tbet2^2 - tbet1^2)/(scbet2^2*sxi2 - scbet1^2*sxi1)
160 // 1-s = -((1-sxi2)*scbet2^2 - (1-sxi1)*scbet1^2)/
161 // (scbet2^2*sxi2 - scbet1^2*sxi1)
162 //
163 // Define phi0 to give same value of s, i.e.,
164 // s = sphi0 * qZ / (m0^2 + sphi0*q0)
165 // = sphi0 * scbet0^2 / (1/qZ + sphi0 * scbet0^2 * sxi0)
166
167 real tphi0, C;
168 if (polar || tphi1 == tphi2) {
169 tphi0 = tphi2;
170 C = 1; // ignored
171 } else {
172 real
173 tbet1 = _fm * tphi1, scbet12 = 1 + Math::sq(tbet1),
174 tbet2 = _fm * tphi2, scbet22 = 1 + Math::sq(tbet2),
175 txi1 = txif(tphi1), cxi1 = 1/hyp(txi1), sxi1 = txi1 * cxi1,
176 txi2 = txif(tphi2), cxi2 = 1/hyp(txi2), sxi2 = txi2 * cxi2,
177 dtbet2 = _fm * (tbet1 + tbet2),
178 es1 = 1 - _e2 * Math::sq(sphi1), es2 = 1 - _e2 * Math::sq(sphi2),
179 /*
180 dsxi = ( (_e2 * sq(sphi2 + sphi1) + es2 + es1) / (2 * es2 * es1) +
181 Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) /
182 ( 2 * _qx ),
183 */
184 dsxi = ( (1 + _e2 * sphi1 * sphi2) / (es2 * es1) +
185 Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) /
186 ( 2 * _qx ),
187 den = (sxi2 + sxi1) * dtbet2 + (scbet22 + scbet12) * dsxi,
188 // s = (sq(tbet2) - sq(tbet1)) / (scbet22*sxi2 - scbet12*sxi1)
189 s = 2 * dtbet2 / den,
190 // 1-s = -(sq(scbet2)*(1-sxi2) - sq(scbet1)*(1-sxi1)) /
191 // (scbet22*sxi2 - scbet12*sxi1)
192 // Write
193 // sq(scbet)*(1-sxi) = sq(scbet)*(1-sphi) * (1-sxi)/(1-sphi)
194 sm1 = -Dsn(tphi2, tphi1, sphi2, sphi1) *
195 ( -( ((sphi2 <= 0 ? (1 - sxi2) / (1 - sphi2) :
196 Math::sq(cxi2/cphi2) * (1 + sphi2) / (1 + sxi2)) +
197 (sphi1 <= 0 ? (1 - sxi1) / (1 - sphi1) :
198 Math::sq(cxi1/cphi1) * (1 + sphi1) / (1 + sxi1))) ) *
199 (1 + _e2 * (sphi1 + sphi2 + sphi1 * sphi2)) /
200 (1 + (sphi1 + sphi2 + sphi1 * sphi2)) +
201 (scbet22 * (sphi2 <= 0 ? 1 - sphi2 :
202 Math::sq(cphi2) / ( 1 + sphi2)) +
203 scbet12 * (sphi1 <= 0 ? 1 - sphi1 : Math::sq(cphi1) / ( 1 + sphi1)))
204 * (_e2 * (1 + sphi1 + sphi2 + _e2 * sphi1 * sphi2)/(es1 * es2)
205 +_e2m * DDatanhee(sphi1, sphi2) ) / _qZ ) / den;
206 // C = (scbet22*sxi2 - scbet12*sxi1) / (scbet22 * scbet12 * (sx2 - sx1))
207 C = den / (2 * scbet12 * scbet22 * dsxi);
208 tphi0 = (tphi2 + tphi1)/2;
209 real stol = tol0_ * fmax(real(1), fabs(tphi0));
210 for (int i = 0;
211 i < 2*numit0_ ||
212 GEOGRAPHICLIB_PANIC("Convergence failure in AlbersEqualArea");
213 ++i) {
214 // Solve (scbet0^2 * sphi0) / (1/qZ + scbet0^2 * sphi0 * sxi0) = s
215 // for tphi0 by Newton's method on
216 // v(tphi0) = (scbet0^2 * sphi0) - s * (1/qZ + scbet0^2 * sphi0 * sxi0)
217 // = 0
218 // Alt:
219 // (scbet0^2 * sphi0) / (1/qZ - scbet0^2 * sphi0 * (1-sxi0))
220 // = s / (1-s)
221 // w(tphi0) = (1-s) * (scbet0^2 * sphi0)
222 // - s * (1/qZ - scbet0^2 * sphi0 * (1-sxi0))
223 // = (1-s) * (scbet0^2 * sphi0)
224 // - S/qZ * (1 - scbet0^2 * sphi0 * (qZ-q0))
225 // Now
226 // qZ-q0 = (1+e2*sphi0)*(1-sphi0)/(1-e2*sphi0^2) +
227 // (1-e2)*atanhee((1-sphi0)/(1-e2*sphi0))
228 // In limit sphi0 -> 1, qZ-q0 -> 2*(1-sphi0)/(1-e2), so wrte
229 // qZ-q0 = 2*(1-sphi0)/(1-e2) + A + B
230 // A = (1-sphi0)*( (1+e2*sphi0)/(1-e2*sphi0^2) - (1+e2)/(1-e2) )
231 // = -e2 *(1-sphi0)^2 * (2+(1+e2)*sphi0) / ((1-e2)*(1-e2*sphi0^2))
232 // B = (1-e2)*atanhee((1-sphi0)/(1-e2*sphi0)) - (1-sphi0)
233 // = (1-sphi0)*(1-e2)/(1-e2*sphi0)*
234 // ((atanhee(x)/x-1) - e2*(1-sphi0)/(1-e2))
235 // x = (1-sphi0)/(1-e2*sphi0), atanhee(x)/x = atanh(e*x)/(e*x)
236 //
237 // 1 - scbet0^2 * sphi0 * (qZ-q0)
238 // = 1 - scbet0^2 * sphi0 * (2*(1-sphi0)/(1-e2) + A + B)
239 // = D - scbet0^2 * sphi0 * (A + B)
240 // D = 1 - scbet0^2 * sphi0 * 2*(1-sphi0)/(1-e2)
241 // = (1-sphi0)*(1-e2*(1+2*sphi0*(1+sphi0)))/((1-e2)*(1+sphi0))
242 // dD/dsphi0 = -2*(1-e2*sphi0^2*(2*sphi0+3))/((1-e2)*(1+sphi0)^2)
243 // d(A+B)/dsphi0 = 2*(1-sphi0^2)*e2*(2-e2*(1+sphi0^2))/
244 // ((1-e2)*(1-e2*sphi0^2)^2)
245
246 real
247 scphi02 = 1 + Math::sq(tphi0), scphi0 = sqrt(scphi02),
248 // sphi0m = 1-sin(phi0) = 1/( sec(phi0) * (tan(phi0) + sec(phi0)) )
249 sphi0 = tphi0 / scphi0, sphi0m = 1/(scphi0 * (tphi0 + scphi0)),
250 // scbet0^2 * sphi0
251 g = (1 + Math::sq( _fm * tphi0 )) * sphi0,
252 // dg/dsphi0 = dg/dtphi0 * scphi0^3
253 dg = _e2m * scphi02 * (1 + 2 * Math::sq(tphi0)) + _e2,
254 D = sphi0m * (1 - _e2*(1 + 2*sphi0*(1+sphi0))) / (_e2m * (1+sphi0)),
255 // dD/dsphi0
256 dD = -2 * (1 - _e2*Math::sq(sphi0) * (2*sphi0+3)) /
257 (_e2m * Math::sq(1+sphi0)),
258 A = -_e2 * Math::sq(sphi0m) * (2+(1+_e2)*sphi0) /
259 (_e2m*(1-_e2*Math::sq(sphi0))),
260 B = (sphi0m * _e2m / (1 - _e2*sphi0) *
261 (atanhxm1(_e2 *
262 Math::sq(sphi0m / (1-_e2*sphi0))) - _e2*sphi0m/_e2m)),
263 // d(A+B)/dsphi0
264 dAB = (2 * _e2 * (2 - _e2 * (1 + Math::sq(sphi0))) /
265 (_e2m * Math::sq(1 - _e2*Math::sq(sphi0)) * scphi02)),
266 u = sm1 * g - s/_qZ * ( D - g * (A + B) ),
267 // du/dsphi0
268 du = sm1 * dg - s/_qZ * (dD - dg * (A + B) - g * dAB),
269 dtu = -u/du * (scphi0 * scphi02);
270 tphi0 += dtu;
271 if (!(fabs(dtu) >= stol))
272 break;
273 }
274 }
275 _txi0 = txif(tphi0); _scxi0 = hyp(_txi0); _sxi0 = _txi0 / _scxi0;
276 _n0 = tphi0/hyp(tphi0);
277 _m02 = 1 / (1 + Math::sq(_fm * tphi0));
278 _nrho0 = polar ? 0 : _a * sqrt(_m02);
279 _k0 = sqrt(tphi1 == tphi2 ? 1 : C / (_m02 + _n0 * _qZ * _sxi0)) * k1;
280 _k2 = Math::sq(_k0);
281 _lat0 = _sign * atan(tphi0)/Math::degree();
282 }
283
284 const AlbersEqualArea& AlbersEqualArea::CylindricalEqualArea() {
285 static const AlbersEqualArea
286 cylindricalequalarea(Constants::WGS84_a(), Constants::WGS84_f(),
287 real(0), real(1), real(0), real(1), real(1));
288 return cylindricalequalarea;
289 }
290
291 const AlbersEqualArea& AlbersEqualArea::AzimuthalEqualAreaNorth() {
292 static const AlbersEqualArea
293 azimuthalequalareanorth(Constants::WGS84_a(), Constants::WGS84_f(),
294 real(1), real(0), real(1), real(0), real(1));
295 return azimuthalequalareanorth;
296 }
297
298 const AlbersEqualArea& AlbersEqualArea::AzimuthalEqualAreaSouth() {
299 static const AlbersEqualArea
300 azimuthalequalareasouth(Constants::WGS84_a(), Constants::WGS84_f(),
301 real(-1), real(0), real(-1), real(0), real(1));
302 return azimuthalequalareasouth;
303 }
304
305 Math::real AlbersEqualArea::txif(real tphi) const {
306 // sxi = ( sphi/(1-e2*sphi^2) + atanhee(sphi) ) /
307 // ( 1/(1-e2) + atanhee(1) )
308 //
309 // txi = ( sphi/(1-e2*sphi^2) + atanhee(sphi) ) /
310 // sqrt( ( (1+e2*sphi)*(1-sphi)/( (1-e2*sphi^2) * (1-e2) ) +
311 // atanhee((1-sphi)/(1-e2*sphi)) ) *
312 // ( (1-e2*sphi)*(1+sphi)/( (1-e2*sphi^2) * (1-e2) ) +
313 // atanhee((1+sphi)/(1+e2*sphi)) ) )
314 // = ( tphi/(1-e2*sphi^2) + atanhee(sphi, e2)/cphi ) /
315 // sqrt(
316 // ( (1+e2*sphi)/( (1-e2*sphi^2) * (1-e2) ) + Datanhee(1, sphi) ) *
317 // ( (1-e2*sphi)/( (1-e2*sphi^2) * (1-e2) ) + Datanhee(1, -sphi) ) )
318 //
319 // This function maintains odd parity
320 real
321 cphi = 1 / sqrt(1 + Math::sq(tphi)),
322 sphi = tphi * cphi,
323 es1 = _e2 * sphi,
324 es2m1 = 1 - es1 * sphi, // 1 - e2 * sphi^2
325 es2m1a = _e2m * es2m1; // (1 - e2 * sphi^2) * (1 - e2)
326 return ( tphi / es2m1 + atanhee(sphi) / cphi ) /
327 sqrt( ( (1 + es1) / es2m1a + Datanhee(1, sphi) ) *
328 ( (1 - es1) / es2m1a + Datanhee(1, -sphi) ) );
329 }
330
331 Math::real AlbersEqualArea::tphif(real txi) const {
332 real
333 tphi = txi,
334 stol = tol_ * fmax(real(1), fabs(txi));
335 // CHECK: min iterations = 1, max iterations = 2; mean = 1.99
336 for (int i = 0;
337 i < numit_ ||
338 GEOGRAPHICLIB_PANIC("Convergence failure in AlbersEqualArea");
339 ++i) {
340 // dtxi/dtphi = (scxi/scphi)^3 * 2*(1-e^2)/(qZ*(1-e^2*sphi^2)^2)
341 real
342 txia = txif(tphi),
343 tphi2 = Math::sq(tphi),
344 scphi2 = 1 + tphi2,
345 scterm = scphi2/(1 + Math::sq(txia)),
346 dtphi = (txi - txia) * scterm * sqrt(scterm) *
347 _qx * Math::sq(1 - _e2 * tphi2 / scphi2);
348 tphi += dtphi;
349 if (!(fabs(dtphi) >= stol))
350 break;
351 }
352 return tphi;
353 }
354
355 // return atanh(sqrt(x))/sqrt(x) - 1 = x/3 + x^2/5 + x^3/7 + ...
356 // typical x < e^2 = 2*f
357 Math::real AlbersEqualArea::atanhxm1(real x) {
358 real s = 0;
359 if (fabs(x) < real(0.5)) {
360 static const real lg2eps_ = -log2(numeric_limits<real>::epsilon() / 2);
361 int e;
362 (void) frexp(x, &e);
363 e = -e;
364 // x = [0.5,1) * 2^(-e)
365 // estimate n s.t. x^n/(2*n+1) < x/3 * epsilon/2
366 // a stronger condition is x^(n-1) < epsilon/2
367 // taking log2 of both sides, a stronger condition is
368 // (n-1)*(-e) < -lg2eps or (n-1)*e > lg2eps or n > ceiling(lg2eps/e)+1
369 int n = x == 0 ? 1 : int(ceil(lg2eps_ / e)) + 1;
370 while (n--) // iterating from n-1 down to 0
371 s = x * s + (n ? 1 : 0)/Math::real(2*n + 1);
372 } else {
373 real xs = sqrt(fabs(x));
374 s = (x > 0 ? atanh(xs) : atan(xs)) / xs - 1;
375 }
376 return s;
377 }
378
379 // return (Datanhee(1,y) - Datanhee(1,x))/(y-x)
380 Math::real AlbersEqualArea::DDatanhee(real x, real y) const {
381 // This function is called with x = sphi1, y = sphi2, phi1 <= phi2, sphi2
382 // >= 0, abs(sphi1) <= phi2. However for safety's sake we enforce x <= y.
383 if (y < x) swap(x, y); // ensure that x <= y
384 real q1 = fabs(_e2),
385 q2 = fabs(2 * _e / _e2m * (1 - x));
386 return
387 x <= 0 || !(fmin(q1, q2) < real(0.75)) ? DDatanhee0(x, y) :
388 (q1 < q2 ? DDatanhee1(x, y) : DDatanhee2(x, y));
389 }
390
391 // Rearrange difference so that 1 - x is in the denominator, then do a
392 // straight divided difference.
393 Math::real AlbersEqualArea::DDatanhee0(real x, real y) const {
394 return (Datanhee(1, y) - Datanhee(x, y))/(1 - x);
395 }
396
397 // The expansion for e2 small
398 Math::real AlbersEqualArea::DDatanhee1(real x, real y) const {
399 // The series in e2 is
400 // sum( c[l] * e2^l, l, 1, N)
401 // where
402 // c[l] = sum( x^i * y^j; i >= 0, j >= 0, i+j < 2*l) / (2*l + 1)
403 // = ( (x-y) - (1-y) * x^(2*l+1) + (1-x) * y^(2*l+1) ) /
404 // ( (2*l+1) * (x-y) * (1-y) * (1-x) )
405 // For x = y = 1, c[l] = l
406 //
407 // In the limit x,y -> 1,
408 //
409 // DDatanhee -> e2/(1-e2)^2 = sum(l * e2^l, l, 1, inf)
410 //
411 // Use if e2 is sufficiently small.
412 real s = 0;
413 real z = 1, k = 1, t = 0, c = 0, en = 1;
414 while (true) {
415 t = y * t + z; c += t; z *= x;
416 t = y * t + z; c += t; z *= x;
417 k += 2; en *= _e2;
418 // Here en[l] = e2^l, k[l] = 2*l + 1,
419 // c[l] = sum( x^i * y^j; i >= 0, j >= 0, i+j < 2*l) / (2*l + 1)
420 // Taylor expansion is
421 // s = sum( c[l] * e2^l, l, 1, N)
422 real ds = en * c / k;
423 s += ds;
424 if (!(fabs(ds) > fabs(s) * eps_/2))
425 break; // Iterate until the added term is sufficiently small
426 }
427 return s;
428 }
429
430 // The expansion for x (and y) close to 1
431 Math::real AlbersEqualArea::DDatanhee2(real x, real y) const {
432 // If x and y are both close to 1, expand in Taylor series in dx = 1-x and
433 // dy = 1-y:
434 //
435 // DDatanhee = sum(C_m * (dx^(m+1) - dy^(m+1)) / (dx - dy), m, 0, inf)
436 //
437 // where
438 //
439 // C_m = sum( (m+2)!! / (m+2-2*k)!! *
440 // ((m+1)/2)! / ((m+1)/2-k)! /
441 // (k! * (2*k-1)!!) *
442 // e2^((m+1)/2+k),
443 // k, 0, (m+1)/2) * (-1)^m / ((m+2) * (1-e2)^(m+2))
444 // for m odd, and
445 //
446 // C_m = sum( 2 * (m+1)!! / (m+1-2*k)!! *
447 // (m/2+1)! / (m/2-k)! /
448 // (k! * (2*k+1)!!) *
449 // e2^(m/2+1+k),
450 // k, 0, m/2)) * (-1)^m / ((m+2) * (1-e2)^(m+2))
451 // for m even.
452 //
453 // Here i!! is the double factorial extended to negative i with
454 // i!! = (i+2)!!/(i+2).
455 //
456 // Note that
457 // (dx^(m+1) - dy^(m+1)) / (dx - dy) =
458 // dx^m + dx^(m-1)*dy ... + dx*dy^(m-1) + dy^m
459 //
460 // Leading (m = 0) term is e2 / (1 - e2)^2
461 //
462 // Magnitude of mth term relative to the leading term scales as
463 //
464 // 2*(2*e/(1-e2)*dx)^m
465 //
466 // So use series if (2*e/(1-e2)*dx) is sufficiently small
467 real s, dx = 1 - x, dy = 1 - y, xy = 1, yy = 1, ee = _e2 / Math::sq(_e2m);
468 s = ee;
469 for (int m = 1; ; ++m) {
470 real c = m + 2, t = c;
471 yy *= dy; // yy = dy^m
472 xy = dx * xy + yy;
473 // Now xy = dx^m + dx^(m-1)*dy ... + dx*dy^(m-1) + dy^m
474 // = (dx^(m+1) - dy^(m+1)) / (dx - dy)
475 // max value = (m+1) * max(dx,dy)^m
476 ee /= -_e2m;
477 if (m % 2 == 0) ee *= _e2;
478 // Now ee = (-1)^m * e2^(floor(m/2)+1) / (1-e2)^(m+2)
479 int kmax = (m+1)/2;
480 for (int k = kmax - 1; k >= 0; --k) {
481 // max coeff is less than 2^(m+1)
482 c *= (k + 1) * (2 * (k + m - 2*kmax) + 3);
483 c /= (kmax - k) * (2 * (kmax - k) + 1);
484 // Horner sum for inner _e2 series
485 t = _e2 * t + c;
486 }
487 // Straight sum for outer m series
488 real ds = t * ee * xy / (m + 2);
489 s = s + ds;
490 if (!(fabs(ds) > fabs(s) * eps_/2))
491 break; // Iterate until the added term is sufficiently small
492 }
493 return s;
494 }
495
496 void AlbersEqualArea::Forward(real lon0, real lat, real lon,
497 real& x, real& y, real& gamma, real& k) const {
498 lon = Math::AngDiff(lon0, lon);
499 lat *= _sign;
500 real sphi, cphi;
501 Math::sincosd(Math::LatFix(lat) * _sign, sphi, cphi);
502 cphi = fmax(epsx_, cphi);
503 real
504 lam = lon * Math::degree(),
505 tphi = sphi/cphi, txi = txif(tphi), sxi = txi/hyp(txi),
506 dq = _qZ * Dsn(txi, _txi0, sxi, _sxi0) * (txi - _txi0),
507 drho = - _a * dq / (sqrt(_m02 - _n0 * dq) + _nrho0 / _a),
508 theta = _k2 * _n0 * lam, stheta = sin(theta), ctheta = cos(theta),
509 t = _nrho0 + _n0 * drho;
510 x = t * (_n0 != 0 ? stheta / _n0 : _k2 * lam) / _k0;
511 y = (_nrho0 *
512 (_n0 != 0 ?
513 (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n0 :
514 0)
515 - drho * ctheta) / _k0;
516 k = _k0 * (t != 0 ? t * hyp(_fm * tphi) / _a : 1);
517 y *= _sign;
518 gamma = _sign * theta / Math::degree();
519 }
520
521 void AlbersEqualArea::Reverse(real lon0, real x, real y,
522 real& lat, real& lon,
523 real& gamma, real& k) const {
524 y *= _sign;
525 real
526 nx = _k0 * _n0 * x, ny = _k0 * _n0 * y, y1 = _nrho0 - ny,
527 den = hypot(nx, y1) + _nrho0, // 0 implies origin with polar aspect
528 drho = den != 0 ? (_k0*x*nx - 2*_k0*y*_nrho0 + _k0*y*ny) / den : 0,
529 // dsxia = scxi0 * dsxi
530 dsxia = - _scxi0 * (2 * _nrho0 + _n0 * drho) * drho /
531 (Math::sq(_a) * _qZ),
532 txi = (_txi0 + dsxia) / sqrt(fmax(1 - dsxia * (2*_txi0 + dsxia), epsx2_)),
533 tphi = tphif(txi),
534 theta = atan2(nx, y1),
535 lam = _n0 != 0 ? theta / (_k2 * _n0) : x / (y1 * _k0);
536 gamma = _sign * theta / Math::degree();
537 lat = Math::atand(_sign * tphi);
538 lon = lam / Math::degree();
539 lon = Math::AngNormalize(lon + Math::AngNormalize(lon0));
540 k = _k0 * (den != 0 ? (_nrho0 + _n0 * drho) * hyp(_fm * tphi) / _a : 1);
541 }
542
543 void AlbersEqualArea::SetScale(real lat, real k) {
544 if (!(isfinite(k) && k > 0))
545 throw GeographicErr("Scale is not positive");
546 if (!(fabs(lat) < Math::qd))
547 throw GeographicErr("Latitude for SetScale not in (-"
548 + to_string(Math::qd) + "d, "
549 + to_string(Math::qd) + "d)");
550 real x, y, gamma, kold;
551 Forward(0, lat, 0, x, y, gamma, kold);
552 k /= kold;
553 _k0 *= k;
554 _k2 = Math::sq(_k0);
555 }
556
557} // namespace GeographicLib
AlbersEqualArea.hpp
Header for GeographicLib::AlbersEqualArea class.
real
GeographicLib::Math::real real
Definition GeodSolve.cpp:28
GEOGRAPHICLIB_PANIC
#define GEOGRAPHICLIB_PANIC(msg)
Definition Math.hpp:62
GeographicLib::AlbersEqualArea
Albers equal area conic projection.
Definition AlbersEqualArea.hpp:60
GeographicLib::AlbersEqualArea::Reverse
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
Definition AlbersEqualArea.cpp:521
GeographicLib::AlbersEqualArea::AlbersEqualArea
AlbersEqualArea(real a, real f, real stdlat, real k0)
Definition AlbersEqualArea.cpp:16
GeographicLib::AlbersEqualArea::SetScale
void SetScale(real lat, real k=real(1))
Definition AlbersEqualArea.cpp:543
GeographicLib::AlbersEqualArea::CylindricalEqualArea
static const AlbersEqualArea & CylindricalEqualArea()
Definition AlbersEqualArea.cpp:284
GeographicLib::AlbersEqualArea::AzimuthalEqualAreaNorth
static const AlbersEqualArea & AzimuthalEqualAreaNorth()
Definition AlbersEqualArea.cpp:291
GeographicLib::AlbersEqualArea::AzimuthalEqualAreaSouth
static const AlbersEqualArea & AzimuthalEqualAreaSouth()
Definition AlbersEqualArea.cpp:298
GeographicLib::AlbersEqualArea::Forward
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
Definition AlbersEqualArea.cpp:496
GeographicLib::GeographicErr
Exception handling for GeographicLib.
Definition Constants.hpp:316
GeographicLib::Math
Mathematical functions needed by GeographicLib.
Definition Math.hpp:77
GeographicLib::Math::degree
static T degree()
Definition Math.hpp:209
GeographicLib::Math::LatFix
static T LatFix(T x)
Definition Math.hpp:309
GeographicLib::Math::sincosd
static void sincosd(T x, T &sinx, T &cosx)
Definition Math.cpp:101
GeographicLib::Math::sq
static T sq(T x)
Definition Math.hpp:221
GeographicLib::Math::qd
static constexpr int qd
degrees per quarter turn
Definition Math.hpp:145
GeographicLib::Math::AngNormalize
static T AngNormalize(T x)
Definition Math.cpp:66
GeographicLib::Math::atand
static T atand(T x)
Definition Math.cpp:218
GeographicLib::Math::AngDiff
static T AngDiff(T x, T y, T &e)
Definition Math.cpp:77
GeographicLib
Namespace for GeographicLib.
Definition Accumulator.cpp:12