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