GeographicLib  1.35
LambertConformalConic.cpp
Go to the documentation of this file.
1 /**
2  * \file LambertConformalConic.cpp
3  * \brief Implementation for GeographicLib::LambertConformalConic 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 LambertConformalConic::eps_ =
17  numeric_limits<real>::epsilon();
18  const Math::real LambertConformalConic::epsx_ = Math::sq(eps_);
19  const Math::real LambertConformalConic::tol_ = real(0.1) * sqrt(eps_);
20  const Math::real LambertConformalConic::ahypover_ =
21  real(numeric_limits<real>::digits) * log(real(numeric_limits<real>::radix))
22  + 2;
23 
25  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  {
33  if (!(Math::isfinite(_a) && _a > 0))
34  throw GeographicErr("Major radius is not positive");
35  if (!(Math::isfinite(_f) && _f < 1))
36  throw GeographicErr("Minor radius is not positive");
37  if (!(Math::isfinite(k0) && k0 > 0))
38  throw GeographicErr("Scale is not positive");
39  if (!(abs(stdlat) <= 90))
40  throw GeographicErr("Standard latitude not in [-90d, 90d]");
41  real
42  phi = stdlat * Math::degree<real>(),
43  sphi = sin(phi),
44  cphi = abs(stdlat) != 90 ? cos(phi) : 0;
45  Init(sphi, cphi, sphi, cphi, k0);
46  }
47 
49  real stdlat1, real stdlat2,
50  real k1)
51  : _a(a)
52  , _f(f <= 1 ? f : 1/f)
53  , _fm(1 - _f)
54  , _e2(_f * (2 - _f))
55  , _e(sqrt(abs(_e2)))
56  , _e2m(1 - _e2)
57  {
58  if (!(Math::isfinite(_a) && _a > 0))
59  throw GeographicErr("Major radius is not positive");
60  if (!(Math::isfinite(_f) && _f < 1))
61  throw GeographicErr("Minor radius is not positive");
62  if (!(Math::isfinite(k1) && k1 > 0))
63  throw GeographicErr("Scale is not positive");
64  if (!(abs(stdlat1) <= 90))
65  throw GeographicErr("Standard latitude 1 not in [-90d, 90d]");
66  if (!(abs(stdlat2) <= 90))
67  throw GeographicErr("Standard latitude 2 not in [-90d, 90d]");
68  real
69  phi1 = stdlat1 * Math::degree<real>(),
70  phi2 = stdlat2 * Math::degree<real>();
71  Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
72  sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
73  }
74 
76  real sinlat1, real coslat1,
77  real sinlat2, real coslat2,
78  real k1)
79  : _a(a)
80  , _f(f <= 1 ? f : 1/f)
81  , _fm(1 - _f)
82  , _e2(_f * (2 - _f))
83  , _e(sqrt(abs(_e2)))
84  , _e2m(1 - _e2)
85  {
86  if (!(Math::isfinite(_a) && _a > 0))
87  throw GeographicErr("Major radius is not positive");
88  if (!(Math::isfinite(_f) && _f < 1))
89  throw GeographicErr("Minor radius is not positive");
90  if (!(Math::isfinite(k1) && k1 > 0))
91  throw GeographicErr("Scale is not positive");
92  if (!(coslat1 >= 0))
93  throw GeographicErr("Standard latitude 1 not in [-90d, 90d]");
94  if (!(coslat2 >= 0))
95  throw GeographicErr("Standard latitude 2 not in [-90d, 90d]");
96  if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
97  throw GeographicErr("Bad sine/cosine of standard latitude 1");
98  if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
99  throw GeographicErr("Bad sine/cosine of standard latitude 2");
100  if (coslat1 == 0 || coslat2 == 0)
101  if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
102  throw GeographicErr
103  ("Standard latitudes must be equal is either is a pole");
104  Init(sinlat1, coslat1, sinlat2, coslat2, k1);
105  }
106 
107  void LambertConformalConic::Init(real sphi1, real cphi1,
108  real sphi2, real cphi2, real k1) throw() {
109  {
110  real r;
111  r = Math::hypot(sphi1, cphi1);
112  sphi1 /= r; cphi1 /= r;
113  r = Math::hypot(sphi2, cphi2);
114  sphi2 /= r; cphi2 /= r;
115  }
116  bool polar = (cphi1 == 0);
117  cphi1 = max(epsx_, cphi1); // Avoid singularities at poles
118  cphi2 = max(epsx_, cphi2);
119  // Determine hemisphere of tangent latitude
120  _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
121  // Internally work with tangent latitude positive
122  sphi1 *= _sign; sphi2 *= _sign;
123  if (sphi1 > sphi2) {
124  swap(sphi1, sphi2); swap(cphi1, cphi2); // Make phi1 < phi2
125  }
126  real
127  tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
128  //
129  // Snyder: 15-8: n = (log(m1) - log(m2))/(log(t1)-log(t2))
130  //
131  // m = cos(bet) = 1/sec(bet) = 1/sqrt(1+tan(bet)^2)
132  // bet = parametric lat, tan(bet) = (1-f)*tan(phi)
133  //
134  // t = tan(pi/4-chi/2) = 1/(sec(chi) + tan(chi)) = sec(chi) - tan(chi)
135  // log(t) = -asinh(tan(chi)) = -psi
136  // chi = conformal lat
137  // tan(chi) = tan(phi)*cosh(xi) - sinh(xi)*sec(phi)
138  // xi = eatanhe(sin(phi)), eatanhe(x) = e * atanh(e*x)
139  //
140  // n = (log(sec(bet2))-log(sec(bet1)))/(asinh(tan(chi2))-asinh(tan(chi1)))
141  //
142  // Let log(sec(bet)) = b(tphi), asinh(tan(chi)) = c(tphi)
143  // Then n = Db(tphi2, tphi1)/Dc(tphi2, tphi1)
144  // In limit tphi2 -> tphi1, n -> sphi1
145  //
146  real
147  tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
148  tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
149  real
150  scphi1 = 1/cphi1,
151  xi1 = eatanhe(sphi1), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
152  tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
153  scphi2 = 1/cphi2,
154  xi2 = eatanhe(sphi2), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
155  tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
156  psi1 = Math::asinh(tchi1);
157  if (tphi2 - tphi1 != 0) {
158  // Db(tphi2, tphi1)
159  real num = Dlog1p(Math::sq(tbet2)/(1 + scbet2),
160  Math::sq(tbet1)/(1 + scbet1))
161  * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
162  // Dc(tphi2, tphi1)
163  real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
164  - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
165  _n = num/den;
166 
167  if (_n < 0.25)
168  _nc = sqrt((1 - _n) * (1 + _n));
169  else {
170  // Compute nc = cos(phi0) = sqrt((1 - n) * (1 + n)), evaluating 1 - n
171  // carefully. First write
172  //
173  // Dc(tphi2, tphi1) * (tphi2 - tphi1)
174  // = log(tchi2 + scchi2) - log(tchi1 + scchi1)
175  //
176  // then den * (1 - n) =
177  // (log((tchi2 + scchi2)/(2*scbet2)) - log((tchi1 + scchi1)/(2*scbet1)))
178  // / (tphi2 - tphi1)
179  // = Dlog1p(a2, a1) * (tchi2+scchi2 + tchi1+scchi1)/(4*scbet1*scbet2)
180  // * fm * Q
181  //
182  // where
183  // a1 = ( (tchi1 - scbet1) + (scchi1 - scbet1) ) / (2 * scbet1)
184  // Q = ((scbet2 + scbet1)/fm)/((scchi2 + scchi1)/D(tchi2, tchi1))
185  // - (tbet2 + tbet1)/(scbet2 + scbet1)
186  real t;
187  {
188  real
189  // s1 = (scbet1 - scchi1) * (scbet1 + scchi1)
190  s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
191  Math::sq(shxi1) * (1 + 2 * Math::sq(tphi1))),
192  s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
193  Math::sq(shxi2) * (1 + 2 * Math::sq(tphi2))),
194  // t1 = scbet1 - tchi1
195  t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
196  t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
197  a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
198  a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
199  t = Dlog1p(a2, a1) / den;
200  }
201  // multiply by (tchi2 + scchi2 + tchi1 + scchi1)/(4*scbet1*scbet2) * fm
202  t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
203  (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
204  (4 * scbet1 * scbet2) ) * _fm;
205 
206  // Rewrite
207  // Q = (1 - (tbet2 + tbet1)/(scbet2 + scbet1)) -
208  // (1 - ((scbet2 + scbet1)/fm)/((scchi2 + scchi1)/D(tchi2, tchi1)))
209  // = tbm - tam
210  // where
211  real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
212  (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
213  (scbet1+scbet2) );
214 
215  // tam = (1 - ((scbet2+scbet1)/fm)/((scchi2+scchi1)/D(tchi2, tchi1)))
216  //
217  // Let
218  // (scbet2 + scbet1)/fm = scphi2 + scphi1 + dbet
219  // (scchi2 + scchi1)/D(tchi2, tchi1) = scphi2 + scphi1 + dchi
220  // then
221  // tam = D(tchi2, tchi1) * (dchi - dbet) / (scchi1 + scchi2)
222  real
223  // D(tchi2, tchi1)
224  dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
225  // (scbet2 + scbet1)/fm - (scphi2 + scphi1)
226  dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
227  1 / (scbet1 + _fm * scphi1) );
228 
229  // dchi = (scchi2 + scchi1)/D(tchi2, tchi1) - (scphi2 + scphi1)
230  // Let
231  // tzet = chxiZ * tphi - shxiZ * scphi
232  // tchi = tzet + nu
233  // scchi = sczet + mu
234  // where
235  // xiZ = eatanhe(1), shxiZ = sinh(xiZ), chxiZ = cosh(xiZ)
236  // nu = scphi * (shxiZ - shxi) - tphi * (chxiZ - chxi)
237  // mu = - scphi * (chxiZ - chxi) + tphi * (shxiZ - shxi)
238  // then
239  // dchi = ((mu2 + mu1) - D(nu2, nu1) * (scphi2 + scphi1)) /
240  // D(tchi2, tchi1)
241  real
242  xiZ = eatanhe(real(1)), shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
243  // These are differences not divided differences
244  // dxiZ1 = xiZ - xi1; dshxiZ1 = shxiZ - shxi; dchxiZ1 = chxiZ - chxi
245  dxiZ1 = Deatanhe(real(1), sphi1)/(scphi1*(tphi1+scphi1)),
246  dxiZ2 = Deatanhe(real(1), sphi2)/(scphi2*(tphi2+scphi2)),
247  dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
248  dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
249  dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
250  dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
251  // mu1 + mu2
252  amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
253  - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
254  // D(xi2, xi1)
255  dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
256  // D(nu2, nu1)
257  dnu12 =
258  ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
259  // Use divided differences
260  (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
261  - ( (scphi1 + scphi2)/2
262  * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
263  // Use ratio of differences
264  (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
265  + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2)
266  * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
267  - (dchxiZ1 + dchxiZ2)/2 ),
268  // dtchi * dchi
269  dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
270  tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
271  t *= tbm - tam;
272  _nc = sqrt(max(real(0), t) * (1 + _n));
273  }
274  {
275  real r = Math::hypot(_n, _nc);
276  _n /= r;
277  _nc /= r;
278  }
279  tphi0 = _n / _nc;
280  } else {
281  tphi0 = tphi1;
282  _nc = 1/hyp(tphi0);
283  _n = tphi0 * _nc;
284  if (polar)
285  _nc = 0;
286  }
287 
288  _scbet0 = hyp(_fm * tphi0);
289  real shxi0 = sinh(eatanhe(_n));
290  _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
291  _psi0 = Math::asinh(_tchi0);
292 
293  _lat0 = atan(_sign * tphi0) / Math::degree<real>();
294  _t0nm1 = Math::expm1(- _n * _psi0); // Snyder's t0^n - 1
295  // a * k1 * m1/t1^n = a * k1 * m2/t2^n = a * k1 * n * (Snyder's F)
296  // = a * k1 / (scbet1 * exp(-n * psi1))
297  _scale = _a * k1 / scbet1 *
298  // exp(n * psi1) = exp(- (1 - n) * psi1) * exp(psi1)
299  // with (1-n) = nc^2/(1+n) and exp(-psi1) = scchi1 + tchi1
300  exp( - (Math::sq(_nc)/(1 + _n)) * psi1 )
301  * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
302  // Scale at phi0 = k0 = k1 * (scbet0*exp(-n*psi0))/(scbet1*exp(-n*psi1))
303  // = k1 * scbet0/scbet1 * exp(n * (psi1 - psi0))
304  // psi1 - psi0 = Dasinh(tchi1, tchi0) * (tchi1 - tchi0)
305  _k0 = k1 * (_scbet0/scbet1) *
306  exp( - (Math::sq(_nc)/(1 + _n)) *
307  Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
308  * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
309  (_scchi0 + _tchi0);
310  _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
311  {
312  // Figure _drhomax using code at beginning of Forward with lat = -90
313  real
314  sphi = -1, cphi = epsx_,
315  tphi = sphi/cphi,
316  scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
317  tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
318  psi = Math::asinh(tchi),
319  dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0);
320  _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ?
321  (exp(Math::sq(_nc)/(1 + _n) * psi ) *
322  (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
323  - (_t0nm1 + 1))/(-_n) :
324  Dexp(-_n * psi, -_n * _psi0) * dpsi);
325  }
326  }
327 
328  const LambertConformalConic
329  LambertConformalConic::Mercator(Constants::WGS84_a<real>(),
330  Constants::WGS84_f<real>(),
331  real(0), real(1));
332 
333  void LambertConformalConic::Forward(real lon0, real lat, real lon,
334  real& x, real& y, real& gamma, real& k)
335  const throw() {
337  // From Snyder, we have
338  //
339  // theta = n * lambda
340  // x = rho * sin(theta)
341  // = (nrho0 + n * drho) * sin(theta)/n
342  // y = rho0 - rho * cos(theta)
343  // = nrho0 * (1-cos(theta))/n - drho * cos(theta)
344  //
345  // where nrho0 = n * rho0, drho = rho - rho0
346  // and drho is evaluated with divided differences
347  real
348  lam = lon * Math::degree<real>(),
349  phi = _sign * lat * Math::degree<real>(),
350  sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx_,
351  tphi = sphi/cphi, scbet = hyp(_fm * tphi),
352  scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
353  tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
354  psi = Math::asinh(tchi),
355  theta = _n * lam, stheta = sin(theta), ctheta = cos(theta),
356  dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0),
357  drho = - _scale * (2 * _nc < 1 && dpsi != 0 ?
358  (exp(Math::sq(_nc)/(1 + _n) * psi ) *
359  (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
360  - (_t0nm1 + 1))/(-_n) :
361  Dexp(-_n * psi, -_n * _psi0) * dpsi);
362  x = (_nrho0 + _n * drho) * (_n != 0 ? stheta / _n : lam);
363  y = _nrho0 *
364  (_n != 0 ?
365  (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n : 0)
366  - drho * ctheta;
367  k = _k0 * (scbet/_scbet0) /
368  (exp( - (Math::sq(_nc)/(1 + _n)) * dpsi )
369  * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
370  y *= _sign;
371  gamma = _sign * theta / Math::degree<real>();
372  }
373 
374  void LambertConformalConic::Reverse(real lon0, real x, real y,
375  real& lat, real& lon,
376  real& gamma, real& k)
377  const throw() {
378  // From Snyder, we have
379  //
380  // x = rho * sin(theta)
381  // rho0 - y = rho * cos(theta)
382  //
383  // rho = hypot(x, rho0 - y)
384  // drho = (n*x^2 - 2*y*nrho0 + n*y^2)/(hypot(n*x, nrho0-n*y) + nrho0)
385  // theta = atan2(n*x, nrho0-n*y)
386  //
387  // From drho, obtain t^n-1
388  // psi = -log(t), so
389  // dpsi = - Dlog1p(t^n-1, t0^n-1) * drho / scale
390  y *= _sign;
391  real
392  // Guard against 0 * inf in computation of ny
393  nx = _n * x, ny = _n ? _n * y : 0, y1 = _nrho0 - ny,
394  den = Math::hypot(nx, y1) + _nrho0, // 0 implies origin with polar aspect
395  // isfinite test is to avoid inf/inf
396  drho = ((den != 0 && Math::isfinite(den))
397  ? (x*nx + y * (ny - 2*_nrho0)) / den
398  : den);
399  drho = min(drho, _drhomax);
400  if (_n == 0)
401  drho = max(drho, -_drhomax);
402  real
403  tnm1 = _t0nm1 + _n * drho/_scale,
404  dpsi = (den == 0 ? 0 :
405  (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
406  ahypover_));
407  real tchi;
408  if (2 * _n <= 1) {
409  // tchi = sinh(psi)
410  real
411  psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
412  dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
413  tchi = _tchi0 + dtchi; // Update tchi using divided difference
414  } else {
415  // tchi = sinh(-1/n * log(tn))
416  // = sinh((1-1/n) * log(tn) - log(tn))
417  // = + sinh((1-1/n) * log(tn)) * cosh(log(tn))
418  // - cosh((1-1/n) * log(tn)) * sinh(log(tn))
419  // (1-1/n) = - nc^2/(n*(1+n))
420  // cosh(log(tn)) = (tn + 1/tn)/2; sinh(log(tn)) = (tn - 1/tn)/2
421  real
422  tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
423  sh = sinh( -Math::sq(_nc)/(_n * (1 + _n)) *
424  (2 * tn > 1 ? Math::log1p(tnm1) : log(tn)) );
425  tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
426  }
427 
428  // Use Newton's method to solve for tphi
429  real
430  // See comment in TransverseMercator.cpp about the initial guess
431  tphi = tchi/_e2m,
432  stol = tol_ * max(real(1), abs(tchi));
433  // min iterations = 1, max iterations = 2; mean = 1.94
434  for (int i = 0; i < numit_; ++i) {
435  real
436  scphi = hyp(tphi),
437  shxi = sinh( eatanhe( tphi / scphi ) ),
438  tchia = hyp(shxi) * tphi - shxi * scphi,
439  dtphi = (tchi - tchia) * (1 + _e2m * Math::sq(tphi)) /
440  ( _e2m * scphi * hyp(tchia) );
441  tphi += dtphi;
442  if (!(abs(dtphi) >= stol))
443  break;
444  }
445  // log(t) = -asinh(tan(chi)) = -psi
446  gamma = atan2(nx, y1);
447  real
448  phi = _sign * atan(tphi),
449  scbet = hyp(_fm * tphi), scchi = hyp(tchi),
450  lam = _n != 0 ? gamma / _n : x / y1;
451  lat = phi / Math::degree<real>();
452  lon = lam / Math::degree<real>();
453  lon = Math::AngNormalize(lon + Math::AngNormalize(lon0));
454  k = _k0 * (scbet/_scbet0) /
455  (exp(_nc != 0 ? - (Math::sq(_nc)/(1 + _n)) * dpsi : 0)
456  * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
457  gamma /= _sign * Math::degree<real>();
458  }
459 
460  void LambertConformalConic::SetScale(real lat, real k) {
461  if (!(Math::isfinite(k) && k > 0))
462  throw GeographicErr("Scale is not positive");
463  if (!(abs(lat) <= 90))
464  throw GeographicErr("Latitude for SetScale not in [-90d, 90d]");
465  if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0))
466  throw GeographicErr("Incompatible polar latitude in SetScale");
467  real x, y, gamma, kold;
468  Forward(0, lat, 0, x, y, gamma, kold);
469  k /= kold;
470  _scale *= k;
471  _k0 *= k;
472  }
473 
474 } // namespace GeographicLib