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