18 : eps_(numeric_limits<real>::epsilon())
19 , epsx_(
Math::sq(eps_))
20 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
25 , _es((_f < 0 ? -1 : 1) * sqrt(fabs(_e2)))
27 if (!(isfinite(_a) && _a > 0))
29 if (!(isfinite(_f) && _f < 1))
31 if (!(isfinite(k0) && k0 > 0))
35 +
"d, " + to_string(
Math::qd) +
"d]");
38 Init(sphi, cphi, sphi, cphi, k0);
42 real stdlat1, real stdlat2,
44 : eps_(numeric_limits<real>::epsilon())
45 , epsx_(
Math::sq(eps_))
46 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
51 , _es((_f < 0 ? -1 : 1) * sqrt(fabs(_e2)))
53 if (!(isfinite(_a) && _a > 0))
55 if (!(isfinite(_f) && _f < 1))
57 if (!(isfinite(k1) && k1 > 0))
67 real sphi1, cphi1, sphi2, cphi2;
70 Init(sphi1, cphi1, sphi2, cphi2, k1);
74 real sinlat1, real coslat1,
75 real sinlat2, real coslat2,
77 : eps_(numeric_limits<real>::epsilon())
78 , epsx_(
Math::sq(eps_))
79 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
84 , _es((_f < 0 ? -1 : 1) * sqrt(fabs(_e2)))
86 if (!(isfinite(_a) && _a > 0))
88 if (!(isfinite(_f) && _f < 1))
90 if (!(isfinite(k1) && k1 > 0))
100 if (!(fabs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
101 throw GeographicErr(
"Bad sine/cosine of standard latitude 1");
102 if (!(fabs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
103 throw GeographicErr(
"Bad sine/cosine of standard latitude 2");
104 if (coslat1 == 0 || coslat2 == 0)
105 if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
107 (
"Standard latitudes must be equal is either is a pole");
108 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
111 void LambertConformalConic::Init(
real sphi1,
real cphi1,
115 r = hypot(sphi1, cphi1);
116 sphi1 /= r; cphi1 /= r;
117 r = hypot(sphi2, cphi2);
118 sphi2 /= r; cphi2 /= r;
120 bool polar = (cphi1 == 0);
121 cphi1 = fmax(epsx_, cphi1);
122 cphi2 = fmax(epsx_, cphi2);
124 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
126 sphi1 *= _sign; sphi2 *= _sign;
128 swap(sphi1, sphi2); swap(cphi1, cphi2);
131 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
151 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
152 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
155 xi1 =
Math::eatanhe(sphi1, _es), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
156 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
158 xi2 =
Math::eatanhe(sphi2, _es), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
159 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
161 if (tphi2 - tphi1 != 0) {
165 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
167 real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
168 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
172 _nc = sqrt((1 - _n) * (1 + _n));
194 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
196 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
199 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
200 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
201 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
202 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
203 t = Dlog1p(a2, a1) / den;
206 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
207 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
208 (4 * scbet1 * scbet2) ) * _fm;
215 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
216 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
228 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
230 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
231 1 / (scbet1 + _fm * scphi1) );
247 shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
250 dxiZ1 = Deatanhe(
real(1), sphi1)/(scphi1*(tphi1+scphi1)),
251 dxiZ2 = Deatanhe(
real(1), sphi2)/(scphi2*(tphi2+scphi2)),
252 dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
253 dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
254 dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
255 dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
257 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
258 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
260 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
263 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
265 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
266 - ( (scphi1 + scphi2)/2
267 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
269 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
270 + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2)
271 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
272 - (dchxiZ1 + dchxiZ2)/2 ),
274 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
275 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
277 _nc = sqrt(fmax(
real(0), t) * (1 + _n));
280 real r = hypot(_n, _nc);
293 _scbet0 = hyp(_fm * tphi0);
295 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
296 _psi0 = asinh(_tchi0);
299 _t0nm1 = expm1(- _n * _psi0);
302 _scale = _a * k1 / scbet1 *
305 exp( - (
Math::sq(_nc)/(1 + _n)) * psi1 )
306 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
310 _k0 = k1 * (_scbet0/scbet1) *
312 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
313 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
315 _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
319 sphi = -1, cphi = epsx_,
322 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
324 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0);
325 _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ?
326 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
327 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
328 - (_t0nm1 + 1))/(-_n) :
329 Dexp(-_n * psi, -_n * _psi0) * dpsi);
342 real& gamma, real& k)
const {
356 cphi = fmax(epsx_, cphi);
359 tphi = sphi/cphi, scbet = hyp(_fm * tphi),
361 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
363 theta = _n * lam, stheta = sin(theta), ctheta = cos(theta),
364 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0),
365 drho = - _scale * (2 * _nc < 1 && dpsi != 0 ?
366 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
367 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
368 - (_t0nm1 + 1))/(-_n) :
369 Dexp(-_n * psi, -_n * _psi0) * dpsi);
370 x = (_nrho0 + _n * drho) * (_n != 0 ? stheta / _n : lam);
373 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) / _n : 0)
375 k = _k0 * (scbet/_scbet0) /
376 (exp( - (
Math::sq(_nc)/(1 + _n)) * dpsi )
377 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
383 real& lat, real& lon,
384 real& gamma, real& k)
const {
400 nx = _n * x, ny = _n != 0 ? _n * y : 0, y1 = _nrho0 - ny,
401 den = hypot(nx, y1) + _nrho0,
403 drho = ((den != 0 && isfinite(den))
404 ? (x*nx + y * (ny - 2*_nrho0)) / den
406 drho = fmin(drho, _drhomax);
408 drho = fmax(drho, -_drhomax);
410 tnm1 = _t0nm1 + _n * drho/_scale,
411 dpsi = (den == 0 ? 0 :
412 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
418 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
419 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
420 tchi = _tchi0 + dtchi;
429 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
430 sh = sinh( -
Math::sq(_nc)/(_n * (1 + _n)) *
431 (2 * tn > 1 ? log1p(tnm1) : log(tn)) );
432 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
436 gamma = atan2(nx, y1);
439 scbet = hyp(_fm * tphi), scchi = hyp(tchi),
440 lam = _n != 0 ? gamma / _n : x / y1;
444 k = _k0 * (scbet/_scbet0) /
445 (exp(_nc != 0 ? - (
Math::sq(_nc)/(1 + _n)) * dpsi : 0)
446 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
451 if (!(isfinite(k) && k > 0))
457 if (fabs(lat) ==
Math::qd && !(_nc == 0 && lat * _n > 0))
458 throw GeographicErr(
"Incompatible polar latitude in SetScale");
459 real x, y, gamma, kold;
460 Forward(0, lat, 0, x, y, gamma, kold);
GeographicLib::Math::real real
Exception handling for GeographicLib.
Mathematical functions needed by GeographicLib.
static void sincosd(T x, T &sinx, T &cosx)
static constexpr int qd
degrees per quarter turn
static T tauf(T taup, T es)
static T AngNormalize(T x)
static T AngDiff(T x, T y, T &e)
static T eatanhe(T x, T es)
Namespace for GeographicLib.