GeographicLib: AlbersEqualArea.cpp Source File (original) (raw)
1
2
3
4
5
6
7
8
9
11
13
14 using namespace std;
15
17 : eps_(numeric_limits::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))
37 if (!(fabs(stdlat) <= Math::qd))
39 + "d, " + to_string(Math::qd) + "d]");
40 real sphi, cphi;
42 Init(sphi, cphi, sphi, cphi, k0);
43 }
44
46 real k1)
47 : eps_(numeric_limits::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))
67 if (!(fabs(stdlat1) <= Math::qd))
68 throw GeographicErr("Standard latitude 1 not in [-"
70 + to_string(Math::qd) + "d]");
71 if (!(fabs(stdlat2) <= Math::qd))
72 throw GeographicErr("Standard latitude 2 not in [-"
74 + to_string(Math::qd) + "d]");
75 real sphi1, cphi1, sphi2, cphi2;
78 Init(sphi1, cphi1, sphi2, cphi2, k1);
79 }
80
82 real sinlat1, real coslat1,
83 real sinlat2, real coslat2,
84 real k1)
85 : eps_(numeric_limits::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))
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)
119 ("Standard latitudes cannot be opposite poles");
120 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
121 }
122
123 void AlbersEqualArea::Init(real sphi1, real cphi1,
125 {
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);
134 cphi2 = fmax(epsx_, cphi2);
135
136 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
137
138 sphi1 *= _sign; sphi2 *= _sign;
139 if (sphi1 > sphi2) {
140 swap(sphi1, sphi2); swap(cphi1, cphi2);
141 }
143 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2;
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167 real tphi0, C;
168 if (polar || tphi1 == tphi2) {
169 tphi0 = tphi2;
170 C = 1;
171 } else {
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
181
182
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
189 s = 2 * dtbet2 / den,
190
191
192
193
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
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_ ||
213 ++i) {
214
215
216
217
218
219
220
221
222
223
224
225
226
227
228
229
230
231
232
233
234
235
236
237
238
239
240
241
242
243
244
245
247 scphi02 = 1 + Math::sq(tphi0), scphi0 = sqrt(scphi02),
248
249 sphi0 = tphi0 / scphi0, sphi0m = 1/(scphi0 * (tphi0 + scphi0)),
250
251 g = (1 + Math::sq( _fm * tphi0 )) * sphi0,
252
253 dg = _e2m * scphi02 * (1 + 2 * Math::sq(tphi0)) + _e2,
254 D = sphi0m * (1 - _e2*(1 + 2*sphi0*(1+sphi0))) / (_e2m * (1+sphi0)),
255
256 dD = -2 * (1 - _e2*Math::sq(sphi0) * (2*sphi0+3)) /
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
264 dAB = (2 * _e2 * (2 - _e2 * (1 + Math::sq(sphi0))) /
266 u = sm1 * g - s/_qZ * ( D - g * (A + B) ),
267
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;
282 }
283
287 real(0), real(1), real(0), real(1), real(1));
288 return cylindricalequalarea;
289 }
290
294 real(1), real(0), real(1), real(0), real(1));
295 return azimuthalequalareanorth;
296 }
297
301 real(-1), real(0), real(-1), real(0), real(1));
302 return azimuthalequalareasouth;
303 }
304
306
307
308
309
310
311
312
313
314
315
316
317
318
319
321 cphi = 1 / sqrt(1 + Math::sq(tphi)),
322 sphi = tphi * cphi,
323 es1 = _e2 * sphi,
324 es2m1 = 1 - es1 * sphi,
325 es2m1a = _e2m * es2m1;
326 return ( tphi / es2m1 + atanhee(sphi) / cphi ) /
327 sqrt( ( (1 + es1) / es2m1a + Datanhee(1, sphi) ) *
328 ( (1 - es1) / es2m1a + Datanhee(1, -sphi) ) );
329 }
330
333 tphi = txi,
334 stol = tol_ * fmax(real(1), fabs(txi));
335
336 for (int i = 0;
337 i < numit_ ||
339 ++i) {
340
342 txia = txif(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
356
359 if (fabs(x) < real(0.5)) {
360 static const real lg2eps_ = -log2(numeric_limits::epsilon() / 2);
361 int e;
362 (void) frexp(x, &e);
363 e = -e;
364
365
366
367
368
369 int n = x == 0 ? 1 : int(ceil(lg2eps_ / e)) + 1;
370 while (n--)
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
381
382
383 if (y < x) swap(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
392
394 return (Datanhee(1, y) - Datanhee(x, y))/(1 - x);
395 }
396
397
399
400
401
402
403
404
405
406
407
408
409
410
411
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
419
420
421
422 real ds = en * c / k;
423 s += ds;
424 if (!(fabs(ds) > fabs(s) * eps_/2))
425 break;
426 }
427 return s;
428 }
429
430
432
433
434
435
436
437
438
439
440
441
442
443
444
445
446
447
448
449
450
451
452
453
454
455
456
457
458
459
460
461
462
463
464
465
466
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;
472 xy = dx * xy + yy;
473
474
475
476 ee /= -_e2m;
477 if (m % 2 == 0) ee *= _e2;
478
479 int kmax = (m+1)/2;
480 for (int k = kmax - 1; k >= 0; --k) {
481
482 c *= (k + 1) * (2 * (k + m - 2*kmax) + 3);
483 c /= (kmax - k) * (2 * (kmax - k) + 1);
484
485 t = _e2 * t + c;
486 }
487
488 real ds = t * ee * xy / (m + 2);
489 s = s + ds;
490 if (!(fabs(ds) > fabs(s) * eps_/2))
491 break;
492 }
493 return s;
494 }
495
497 real& x, real& y, real& gamma, real& k) const {
499 lat *= _sign;
500 real sphi, cphi;
502 cphi = fmax(epsx_, cphi);
503 real
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;
519 }
520
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,
528 drho = den != 0 ? (_k0*x*nx - 2*_k0*y*_nrho0 + _k0*y*ny) / den : 0,
529
530 dsxia = - _scxi0 * (2 * _nrho0 + _n0 * drho) * drho /
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);
540 k = _k0 * (den != 0 ? (_nrho0 + _n0 * drho) * hyp(_fm * tphi) / _a : 1);
541 }
542
544 if (!(isfinite(k) && k > 0))
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;
555 }
556
557}
Header for GeographicLib::AlbersEqualArea class.
GeographicLib::Math::real real
#define GEOGRAPHICLIB_PANIC(msg)
Albers equal area conic projection.
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
Definition AlbersEqualArea.cpp:521
AlbersEqualArea(real a, real f, real stdlat, real k0)
Definition AlbersEqualArea.cpp:16
void SetScale(real lat, real k=real(1))
Definition AlbersEqualArea.cpp:543
static const AlbersEqualArea & CylindricalEqualArea()
Definition AlbersEqualArea.cpp:284
static const AlbersEqualArea & AzimuthalEqualAreaNorth()
Definition AlbersEqualArea.cpp:291
static const AlbersEqualArea & AzimuthalEqualAreaSouth()
Definition AlbersEqualArea.cpp:298
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
Definition AlbersEqualArea.cpp:496
Exception handling for GeographicLib.
Mathematical functions needed by GeographicLib.
static void sincosd(T x, T &sinx, T &cosx)
static constexpr int qd
degrees per quarter turn
static T AngNormalize(T x)
static T AngDiff(T x, T y, T &e)
Namespace for GeographicLib.
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)