12 namespace GeographicLib {
17 numeric_limits<real>::epsilon();
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))
27 , _f(f <= 1 ? f : 1/f)
39 if (!(abs(stdlat) <= 90))
42 phi = stdlat * Math::degree<real>(),
44 cphi = abs(stdlat) != 90 ? cos(phi) : 0;
45 Init(sphi, cphi, sphi, cphi, k0);
49 real stdlat1, real stdlat2,
52 , _f(f <= 1 ? f : 1/f)
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]");
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);
76 real sinlat1, real coslat1,
77 real sinlat2, real coslat2,
80 , _f(f <= 1 ? f : 1/f)
93 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
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))
103 (
"Standard latitudes must be equal is either is a pole");
104 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
107 void LambertConformalConic::Init(
real sphi1,
real cphi1,
112 sphi1 /= r; cphi1 /= r;
114 sphi2 /= r; cphi2 /= r;
116 bool polar = (cphi1 == 0);
117 cphi1 = max(epsx_, cphi1);
118 cphi2 = max(epsx_, cphi2);
120 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
122 sphi1 *= _sign; sphi2 *= _sign;
124 swap(sphi1, sphi2); swap(cphi1, cphi2);
127 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
147 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
148 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
151 xi1 = eatanhe(sphi1), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
152 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
154 xi2 = eatanhe(sphi2), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
155 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
157 if (tphi2 - tphi1 != 0) {
161 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
163 real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
164 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
168 _nc = sqrt((1 - _n) * (1 + _n));
190 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
192 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
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;
202 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
203 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
204 (4 * scbet1 * scbet2) ) * _fm;
211 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
212 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
224 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
226 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
227 1 / (scbet1 + _fm * scphi1) );
242 xiZ = eatanhe(
real(1)), shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
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,
252 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
253 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
255 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
258 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
260 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
261 - ( (scphi1 + scphi2)/2
262 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
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 ),
269 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
270 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
272 _nc = sqrt(max(
real(0), t) * (1 + _n));
288 _scbet0 = hyp(_fm * tphi0);
289 real shxi0 = sinh(eatanhe(_n));
290 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
293 _lat0 = atan(_sign * tphi0) / Math::degree<real>();
297 _scale = _a * k1 / scbet1 *
300 exp( - (
Math::sq(_nc)/(1 + _n)) * psi1 )
301 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
305 _k0 = k1 * (_scbet0/scbet1) *
307 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
308 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
310 _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
314 sphi = -1, cphi = epsx_,
316 scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
317 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(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);
328 const LambertConformalConic
330 Constants::WGS84_f<real>(),
334 real& x, real& y, real& gamma, real& k)
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),
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);
365 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) / _n : 0)
367 k = _k0 * (scbet/_scbet0) /
368 (exp( - (
Math::sq(_nc)/(1 + _n)) * dpsi )
369 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
371 gamma = _sign * theta / Math::degree<real>();
375 real& lat, real& lon,
376 real& gamma, real& k)
393 nx = _n * x, ny = _n ? _n * y : 0, y1 = _nrho0 - ny,
397 ? (x*nx + y * (ny - 2*_nrho0)) / den
399 drho = min(drho, _drhomax);
401 drho = max(drho, -_drhomax);
403 tnm1 = _t0nm1 + _n * drho/_scale,
404 dpsi = (den == 0 ? 0 :
405 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
411 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
412 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
413 tchi = _tchi0 + dtchi;
422 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
423 sh = sinh( -
Math::sq(_nc)/(_n * (1 + _n)) *
425 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
432 stol = tol_ * max(
real(1), abs(tchi));
434 for (
int i = 0; i < numit_; ++i) {
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) );
442 if (!(abs(dtphi) >= stol))
446 gamma = atan2(nx, y1);
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>();
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>();
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);
static T AngNormalize(T x)
GeographicLib::Math::real real
static bool isfinite(T x)
static T AngDiff(T x, T y)
Exception handling for GeographicLib.