GeographicLib  1.35
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-2012) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  **********************************************************************/
9 
11 
12 namespace GeographicLib {
13 
14  using namespace std;
15 
16  const Math::real AlbersEqualArea::eps_ = numeric_limits<real>::epsilon();
17  const Math::real AlbersEqualArea::epsx_ = Math::sq(eps_);
18  const Math::real AlbersEqualArea::epsx2_ = Math::sq(epsx_);
19  const Math::real AlbersEqualArea::tol_ = sqrt(eps_);
20  const Math::real AlbersEqualArea::tol0_ = tol_ * sqrt(sqrt(eps_));
21  const Math::real AlbersEqualArea::ahypover_ =
22  real(numeric_limits<real>::digits) * log(real(numeric_limits<real>::radix))
23  + 2;
24 
25  AlbersEqualArea::AlbersEqualArea(real a, real f, real stdlat, real k0)
26  : _a(a)
27  , _f(f <= 1 ? f : 1/f)
28  , _fm(1 - _f)
29  , _e2(_f * (2 - _f))
30  , _e(sqrt(abs(_e2)))
31  , _e2m(1 - _e2)
32  , _qZ(1 + _e2m * atanhee(real(1)))
33  , _qx(_qZ / ( 2 * _e2m ))
34  {
35  if (!(Math::isfinite(_a) && _a > 0))
36  throw GeographicErr("Major radius is not positive");
37  if (!(Math::isfinite(_f) && _f < 1))
38  throw GeographicErr("Minor radius is not positive");
39  if (!(Math::isfinite(k0) && k0 > 0))
40  throw GeographicErr("Scale is not positive");
41  if (!(abs(stdlat) <= 90))
42  throw GeographicErr("Standard latitude not in [-90d, 90d]");
43  real
44  phi = stdlat * Math::degree<real>(),
45  sphi = sin(phi),
46  cphi = abs(stdlat) != 90 ? cos(phi) : 0;
47  Init(sphi, cphi, sphi, cphi, k0);
48  }
49 
50  AlbersEqualArea::AlbersEqualArea(real a, real f, real stdlat1, real stdlat2,
51  real k1)
52  : _a(a)
53  , _f(f <= 1 ? f : 1/f)
54  , _fm(1 - _f)
55  , _e2(_f * (2 - _f))
56  , _e(sqrt(abs(_e2)))
57  , _e2m(1 - _e2)
58  , _qZ(1 + _e2m * atanhee(real(1)))
59  , _qx(_qZ / ( 2 * _e2m ))
60  {
61  if (!(Math::isfinite(_a) && _a > 0))
62  throw GeographicErr("Major radius is not positive");
63  if (!(Math::isfinite(_f) && _f < 1))
64  throw GeographicErr("Minor radius is not positive");
65  if (!(Math::isfinite(k1) && k1 > 0))
66  throw GeographicErr("Scale is not positive");
67  if (!(abs(stdlat1) <= 90))
68  throw GeographicErr("Standard latitude 1 not in [-90d, 90d]");
69  if (!(abs(stdlat2) <= 90))
70  throw GeographicErr("Standard latitude 2 not in [-90d, 90d]");
71  real
72  phi1 = stdlat1 * Math::degree<real>(),
73  phi2 = stdlat2 * Math::degree<real>();
74  Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
75  sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
76  }
77 
79  real sinlat1, real coslat1,
80  real sinlat2, real coslat2,
81  real k1)
82  : _a(a)
83  , _f(f <= 1 ? f : 1/f)
84  , _fm(1 - _f)
85  , _e2(_f * (2 - _f))
86  , _e(sqrt(abs(_e2)))
87  , _e2m(1 - _e2)
88  , _qZ(1 + _e2m * atanhee(real(1)))
89  , _qx(_qZ / ( 2 * _e2m ))
90  {
91  if (!(Math::isfinite(_a) && _a > 0))
92  throw GeographicErr("Major radius is not positive");
93  if (!(Math::isfinite(_f) && _f < 1))
94  throw GeographicErr("Minor radius is not positive");
95  if (!(Math::isfinite(k1) && k1 > 0))
96  throw GeographicErr("Scale is not positive");
97  if (!(coslat1 >= 0))
98  throw GeographicErr("Standard latitude 1 not in [-90d, 90d]");
99  if (!(coslat2 >= 0))
100  throw GeographicErr("Standard latitude 2 not in [-90d, 90d]");
101  if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
102  throw GeographicErr("Bad sine/cosine of standard latitude 1");
103  if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
104  throw GeographicErr("Bad sine/cosine of standard latitude 2");
105  if (coslat1 == 0 && coslat2 == 0 && sinlat1 * sinlat2 <= 0)
106  throw GeographicErr
107  ("Standard latitudes cannot be opposite poles");
108  Init(sinlat1, coslat1, sinlat2, coslat2, k1);
109  }
110 
111  void AlbersEqualArea::Init(real sphi1, real cphi1,
112  real sphi2, real cphi2, real k1) throw() {
113  {
114  real r;
115  r = Math::hypot(sphi1, cphi1);
116  sphi1 /= r; cphi1 /= r;
117  r = Math::hypot(sphi2, cphi2);
118  sphi2 /= r; cphi2 /= r;
119  }
120  bool polar = (cphi1 == 0);
121  cphi1 = max(epsx_, cphi1); // Avoid singularities at poles
122  cphi2 = max(epsx_, cphi2);
123  // Determine hemisphere of tangent latitude
124  _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
125  // Internally work with tangent latitude positive
126  sphi1 *= _sign; sphi2 *= _sign;
127  if (sphi1 > sphi2) {
128  swap(sphi1, sphi2); swap(cphi1, cphi2); // Make phi1 < phi2
129  }
130  real
131  tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2;
132 
133  // q = (1-e^2)*(sphi/(1-e^2*sphi^2) - atanhee(sphi))
134  // qZ = q(pi/2) = (1 + (1-e^2)*atanhee(1))
135  // atanhee(x) = atanh(e*x)/e
136  // q = sxi * qZ
137  // dq/dphi = 2*(1-e^2)*cphi/(1-e^2*sphi^2)^2
138  //
139  // n = (m1^2-m2^2)/(q2-q1) -> sin(phi0) for phi1, phi2 -> phi0
140  // C = m1^2 + n*q1 = (m1^2*q2-m2^2*q1)/(q2-q1)
141  // let
142  // rho(pi/2)/rho(-pi/2) = (1-s)/(1+s)
143  // s = n*qZ/C
144  // = qZ * (m1^2-m2^2)/(m1^2*q2-m2^2*q1)
145  // = qZ * (scbet2^2 - scbet1^2)/(scbet2^2*q2 - scbet1^2*q1)
146  // = (scbet2^2 - scbet1^2)/(scbet2^2*sxi2 - scbet1^2*sxi1)
147  // = (tbet2^2 - tbet1^2)/(scbet2^2*sxi2 - scbet1^2*sxi1)
148  // 1-s = -((1-sxi2)*scbet2^2 - (1-sxi1)*scbet1^2)/
149  // (scbet2^2*sxi2 - scbet1^2*sxi1)
150  //
151  // Define phi0 to give same value of s, i.e.,
152  // s = sphi0 * qZ / (m0^2 + sphi0*q0)
153  // = sphi0 * scbet0^2 / (1/qZ + sphi0 * scbet0^2 * sxi0)
154 
155  real tphi0, C;
156  if (polar || tphi1 == tphi2) {
157  tphi0 = tphi2;
158  C = 1; // ignored
159  } else {
160  real
161  tbet1 = _fm * tphi1, scbet12 = 1 + Math::sq(tbet1),
162  tbet2 = _fm * tphi2, scbet22 = 1 + Math::sq(tbet2),
163  txi1 = txif(tphi1), cxi1 = 1/hyp(txi1), sxi1 = txi1 * cxi1,
164  txi2 = txif(tphi2), cxi2 = 1/hyp(txi2), sxi2 = txi2 * cxi2,
165  dtbet2 = _fm * (tbet1 + tbet2),
166  es1 = 1 - _e2 * Math::sq(sphi1), es2 = 1 - _e2 * Math::sq(sphi2),
167  /*
168  dsxi = ( (_e2 * sq(sphi2 + sphi1) + es2 + es1) / (2 * es2 * es1) +
169  Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) /
170  ( 2 * _qx ),
171  */
172  dsxi = ( (1 + _e2 * sphi1 * sphi2) / (es2 * es1) +
173  Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) /
174  ( 2 * _qx ),
175  den = (sxi2 + sxi1) * dtbet2 + (scbet22 + scbet12) * dsxi,
176  // s = (sq(tbet2) - sq(tbet1)) / (scbet22*sxi2 - scbet12*sxi1)
177  s = 2 * dtbet2 / den,
178  // 1-s = -(sq(scbet2)*(1-sxi2) - sq(scbet1)*(1-sxi1)) /
179  // (scbet22*sxi2 - scbet12*sxi1)
180  // Write
181  // sq(scbet)*(1-sxi) = sq(scbet)*(1-sphi) * (1-sxi)/(1-sphi)
182  sm1 = -Dsn(tphi2, tphi1, sphi2, sphi1) *
183  ( -( ((sphi2 <= 0 ? (1 - sxi2) / (1 - sphi2) :
184  Math::sq(cxi2/cphi2) * (1 + sphi2) / (1 + sxi2)) +
185  (sphi1 <= 0 ? (1 - sxi1) / (1 - sphi1) :
186  Math::sq(cxi1/cphi1) * (1 + sphi1) / (1 + sxi1))) ) *
187  (1 + _e2 * (sphi1 + sphi2 + sphi1 * sphi2)) /
188  (1 + (sphi1 + sphi2 + sphi1 * sphi2)) +
189  (scbet22 * (sphi2 <= 0 ? 1 - sphi2 : Math::sq(cphi2) / ( 1 + sphi2)) +
190  scbet12 * (sphi1 <= 0 ? 1 - sphi1 : Math::sq(cphi1) / ( 1 + sphi1)))
191  * (_e2 * (1 + sphi1 + sphi2 + _e2 * sphi1 * sphi2)/(es1 * es2)
192  +_e2m * DDatanhee(sphi1, sphi2) ) / _qZ ) / den;
193  // C = (scbet22*sxi2 - scbet12*sxi1) / (scbet22 * scbet12 * (sx2 - sx1))
194  C = den / (2 * scbet12 * scbet22 * dsxi);
195  tphi0 = (tphi2 + tphi1)/2;
196  real stol = tol0_ * max(real(1), abs(tphi0));
197  for (int i = 0; i < 2*numit0_; ++i) {
198  // Solve (scbet0^2 * sphi0) / (1/qZ + scbet0^2 * sphi0 * sxi0) = s
199  // for tphi0 by Newton's method on
200  // v(tphi0) = (scbet0^2 * sphi0) - s * (1/qZ + scbet0^2 * sphi0 * sxi0)
201  // = 0
202  // Alt:
203  // (scbet0^2 * sphi0) / (1/qZ - scbet0^2 * sphi0 * (1-sxi0)) = s / (1-s)
204  // w(tphi0) = (1-s) * (scbet0^2 * sphi0)
205  // - s * (1/qZ - scbet0^2 * sphi0 * (1-sxi0))
206  // = (1-s) * (scbet0^2 * sphi0)
207  // - S/qZ * (1 - scbet0^2 * sphi0 * (qZ-q0))
208  // Now
209  // qZ-q0 = (1+e2*sphi0)*(1-sphi0)/(1-e2*sphi0^2) +
210  // (1-e2)*atanhee((1-sphi0)/(1-e2*sphi0))
211  // In limit sphi0 -> 1, qZ-q0 -> 2*(1-sphi0)/(1-e2), so wrte
212  // qZ-q0 = 2*(1-sphi0)/(1-e2) + A + B
213  // A = (1-sphi0)*( (1+e2*sphi0)/(1-e2*sphi0^2) - (1+e2)/(1-e2) )
214  // = -e2 *(1-sphi0)^2 * (2+(1+e2)*sphi0) / ((1-e2)*(1-e2*sphi0^2))
215  // B = (1-e2)*atanhee((1-sphi0)/(1-e2*sphi0)) - (1-sphi0)
216  // = (1-sphi0)*(1-e2)/(1-e2*sphi0)*
217  // ((atanhee(x)/x-1) - e2*(1-sphi0)/(1-e2))
218  // x = (1-sphi0)/(1-e2*sphi0), atanhee(x)/x = atanh(e*x)/(e*x)
219  //
220  // 1 - scbet0^2 * sphi0 * (qZ-q0)
221  // = 1 - scbet0^2 * sphi0 * (2*(1-sphi0)/(1-e2) + A + B)
222  // = D - scbet0^2 * sphi0 * (A + B)
223  // D = 1 - scbet0^2 * sphi0 * 2*(1-sphi0)/(1-e2)
224  // = (1-sphi0)*(1-e2*(1+2*sphi0*(1+sphi0)))/((1-e2)*(1+sphi0))
225  // dD/dsphi0 = -2*(1-e2*sphi0^2*(2*sphi0+3))/((1-e2)*(1+sphi0)^2)
226  // d(A+B)/dsphi0 = 2*(1-sphi0^2)*e2*(2-e2*(1+sphi0^2))/
227  // ((1-e2)*(1-e2*sphi0^2)^2)
228 
229  real
230  scphi02 = 1 + Math::sq(tphi0), scphi0 = sqrt(scphi02),
231  // sphi0m = 1-sin(phi0) = 1/( sec(phi0) * (tan(phi0) + sec(phi0)) )
232  sphi0 = tphi0 / scphi0, sphi0m = 1/(scphi0 * (tphi0 + scphi0)),
233  // scbet0^2 * sphi0
234  g = (1 + Math::sq( _fm * tphi0 )) * sphi0,
235  // dg/dsphi0 = dg/dtphi0 * scphi0^3
236  dg = _e2m * scphi02 * (1 + 2 * Math::sq(tphi0)) + _e2,
237  D = sphi0m * (1 - _e2*(1 + 2*sphi0*(1+sphi0))) / (_e2m * (1+sphi0)),
238  // dD/dsphi0
239  dD = -2 * (1 - _e2*Math::sq(sphi0) * (2*sphi0+3)) /
240  (_e2m * Math::sq(1+sphi0)),
241  A = -_e2 * Math::sq(sphi0m) * (2+(1+_e2)*sphi0) /
242  (_e2m*(1-_e2*Math::sq(sphi0))),
243  B = (sphi0m * _e2m / (1 - _e2*sphi0) *
244  (atanhxm1(_e2 *
245  Math::sq(sphi0m / (1-_e2*sphi0))) - _e2*sphi0m/_e2m)),
246  // d(A+B)/dsphi0
247  dAB = (2 * _e2 * (2 - _e2 * (1 + Math::sq(sphi0))) /
248  (_e2m * Math::sq(1 - _e2*Math::sq(sphi0)) * scphi02)),
249  u = sm1 * g - s/_qZ * ( D - g * (A + B) ),
250  // du/dsphi0
251  du = sm1 * dg - s/_qZ * (dD - dg * (A + B) - g * dAB),
252  dtu = -u/du * (scphi0 * scphi02);
253  tphi0 += dtu;
254  if (!(abs(dtu) >= stol))
255  break;
256  }
257  }
258  _txi0 = txif(tphi0); _scxi0 = hyp(_txi0); _sxi0 = _txi0 / _scxi0;
259  _n0 = tphi0/hyp(tphi0);
260  _m02 = 1 / (1 + Math::sq(_fm * tphi0));
261  _nrho0 = polar ? 0 : _a * sqrt(_m02);
262  _k0 = sqrt(tphi1 == tphi2 ? 1 : C / (_m02 + _n0 * _qZ * _sxi0)) * k1;
263  _k2 = Math::sq(_k0);
264  _lat0 = _sign * atan(tphi0)/Math::degree<real>();
265  }
266 
267  const AlbersEqualArea
268  AlbersEqualArea::CylindricalEqualArea(Constants::WGS84_a<real>(),
269  Constants::WGS84_f<real>(),
270  real(0), real(1), real(0), real(1),
271  real(1));
272 
273  const AlbersEqualArea
274  AlbersEqualArea::AzimuthalEqualAreaNorth(Constants::WGS84_a<real>(),
275  Constants::WGS84_f<real>(),
276  real(1), real(0), real(1), real(0),
277  real(1));
278 
279  const AlbersEqualArea
280  AlbersEqualArea::AzimuthalEqualAreaSouth(Constants::WGS84_a<real>(),
281  Constants::WGS84_f<real>(),
282  real(-1), real(0), real(-1), real(0),
283  real(1));
284 
285  Math::real AlbersEqualArea::txif(real tphi) const throw() {
286  // sxi = ( sphi/(1-e2*sphi^2) + atanhee(sphi) ) /
287  // ( 1/(1-e2) + atanhee(1) )
288  //
289  // txi = ( sphi/(1-e2*sphi^2) + atanhee(sphi) ) /
290  // sqrt( ( (1+e2*sphi)*(1-sphi)/( (1-e2*sphi^2) * (1-e2) ) +
291  // atanhee((1-sphi)/(1-e2*sphi)) ) *
292  // ( (1-e2*sphi)*(1+sphi)/( (1-e2*sphi^2) * (1-e2) ) +
293  // atanhee((1+sphi)/(1+e2*sphi)) ) )
294  //
295  // subst 1-sphi = cphi^2/(1+sphi)
296  int s = tphi < 0 ? -1 : 1; // Enforce odd parity
297  tphi *= s;
298  real
299  cphi2 = 1 / (1 + Math::sq(tphi)),
300  sphi = tphi * sqrt(cphi2),
301  es1 = _e2 * sphi,
302  es2m1 = 1 - es1 * sphi,
303  sp1 = 1 + sphi,
304  es1m1 = (1 - es1) * sp1,
305  es2m1a = _e2m * es2m1,
306  es1p1 = sp1 / (1 + es1);
307  return s * ( sphi / es2m1 + atanhee(sphi) ) /
308  sqrt( ( cphi2 / (es1p1 * es2m1a) + atanhee(cphi2 / es1m1) ) *
309  ( es1m1 / es2m1a + atanhee(es1p1) ) );
310  }
311 
312  Math::real AlbersEqualArea::tphif(real txi) const throw() {
313  real
314  tphi = txi,
315  stol = tol_ * max(real(1), abs(txi));
316  // CHECK: min iterations = 1, max iterations = 2; mean = 1.99
317  for (int i = 0; i < numit_; ++i) {
318  // dtxi/dtphi = (scxi/scphi)^3 * 2*(1-e^2)/(qZ*(1-e^2*sphi^2)^2)
319  real
320  txia = txif(tphi),
321  tphi2 = Math::sq(tphi),
322  scphi2 = 1 + tphi2,
323  scterm = scphi2/(1 + Math::sq(txia)),
324  dtphi = (txi - txia) * scterm * sqrt(scterm) *
325  _qx * Math::sq(1 - _e2 * tphi2 / scphi2);
326  tphi += dtphi;
327  if (!(abs(dtphi) >= stol))
328  break;
329  }
330  return tphi;
331  }
332 
333  // return atanh(sqrt(x))/sqrt(x) - 1 = y/3 + y^2/5 + y^3/7 + ...
334  // typical x < e^2 = 2*f
335  Math::real AlbersEqualArea::atanhxm1(real x) throw() {
336  real s = 0;
337  if (abs(x) < real(0.5)) {
338  real os = -1, y = 1, k = 1;
339  while (os != s) {
340  os = s;
341  y *= x; // y = x^n
342  k += 2; // k = 2*n + 1
343  s += y/k; // sum( x^n/(2*n + 1) )
344  }
345  } else {
346  real xs = sqrt(abs(x));
347  s = (x > 0 ? Math::atanh(xs) : atan(xs)) / xs - 1;
348  }
349  return s;
350  }
351 
352  // return (Datanhee(1,y) - Datanhee(1,x))/(y-x)
353  Math::real AlbersEqualArea::DDatanhee(real x, real y) const throw() {
354  real s = 0;
355  if (_e2 * (abs(x) + abs(y)) < real(0.5)) {
356  real os = -1, z = 1, k = 1, t = 0, c = 0, en = 1;
357  while (os != s) {
358  os = s;
359  t = y * t + z; c += t; z *= x;
360  t = y * t + z; c += t; z *= x;
361  k += 2; en *= _e2;
362  // Here en[l] = e2^l, k[l] = 2*l + 1,
363  // c[l] = sum( x^i * y^j; i >= 0, j >= 0, i+j < 2*l)
364  s += en * c / k;
365  }
366  // Taylor expansion is
367  // s = sum( c[l] * e2^l / (2*l + 1), l, 1, N)
368  } else
369  s = (Datanhee(1, y) - Datanhee(x, y))/(1 - x);
370  return s;
371  }
372 
373  void AlbersEqualArea::Forward(real lon0, real lat, real lon,
374  real& x, real& y, real& gamma, real& k)
375  const throw() {
377  lat *= _sign;
378  real
379  lam = lon * Math::degree<real>(),
380  phi = lat * Math::degree<real>(),
381  sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx_,
382  tphi = sphi/cphi, txi = txif(tphi), sxi = txi/hyp(txi),
383  dq = _qZ * Dsn(txi, _txi0, sxi, _sxi0) * (txi - _txi0),
384  drho = - _a * dq / (sqrt(_m02 - _n0 * dq) + _nrho0 / _a),
385  theta = _k2 * _n0 * lam, stheta = sin(theta), ctheta = cos(theta),
386  t = _nrho0 + _n0 * drho;
387  x = t * (_n0 != 0 ? stheta / _n0 : _k2 * lam) / _k0;
388  y = (_nrho0 *
389  (_n0 != 0 ?
390  (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n0 :
391  0)
392  - drho * ctheta) / _k0;
393  k = _k0 * (t != 0 ? t * hyp(_fm * tphi) / _a : 1);
394  y *= _sign;
395  gamma = _sign * theta / Math::degree<real>();
396  }
397 
398  void AlbersEqualArea::Reverse(real lon0, real x, real y,
399  real& lat, real& lon,
400  real& gamma, real& k)
401  const throw() {
402  y *= _sign;
403  real
404  nx = _k0 * _n0 * x, ny = _k0 * _n0 * y, y1 = _nrho0 - ny,
405  den = Math::hypot(nx, y1) + _nrho0, // 0 implies origin with polar aspect
406  drho = den != 0 ? (_k0*x*nx - 2*_k0*y*_nrho0 + _k0*y*ny) / den : 0,
407  // dsxia = scxi0 * dsxi
408  dsxia = - _scxi0 * (2 * _nrho0 + _n0 * drho) * drho /
409  (Math::sq(_a) * _qZ),
410  txi = (_txi0 + dsxia) / sqrt(max(1 - dsxia * (2*_txi0 + dsxia), epsx2_)),
411  tphi = tphif(txi),
412  phi = _sign * atan(tphi),
413  theta = atan2(nx, y1),
414  lam = _n0 != 0 ? theta / (_k2 * _n0) : x / (y1 * _k0);
415  gamma = _sign * theta / Math::degree<real>();
416  lat = phi / Math::degree<real>();
417  lon = lam / Math::degree<real>();
418  lon = Math::AngNormalize(lon + Math::AngNormalize(lon0));
419  k = _k0 * (den != 0 ? (_nrho0 + _n0 * drho) * hyp(_fm * tphi) / _a : 1);
420  }
421 
422  void AlbersEqualArea::SetScale(real lat, real k) {
423  if (!(Math::isfinite(k) && k > 0))
424  throw GeographicErr("Scale is not positive");
425  if (!(abs(lat) < 90))
426  throw GeographicErr("Latitude for SetScale not in (-90d, 90d)");
427  real x, y, gamma, kold;
428  Forward(0, lat, 0, x, y, gamma, kold);
429  k /= kold;
430  _k0 *= k;
431  _k2 = Math::sq(_k0);
432  }
433 
434 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:388
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
static const AlbersEqualArea AzimuthalEqualAreaSouth
GeographicLib::Math::real real
Definition: GeodSolve.cpp:40
AlbersEqualArea(real a, real f, real stdlat, real k0)
static bool isfinite(T x)
Definition: Math.hpp:435
static const AlbersEqualArea CylindricalEqualArea
static T atanh(T x)
Definition: Math.hpp:315
static const AlbersEqualArea AzimuthalEqualAreaNorth
Header for GeographicLib::AlbersEqualArea class.
static T hypot(T x, T y)
Definition: Math.hpp:165
static T sq(T x)
Definition: Math.hpp:153
static T AngDiff(T x, T y)
Definition: Math.hpp:418
Exception handling for GeographicLib.
Definition: Constants.hpp:320
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
void SetScale(real lat, real k=real(1))