14 # pragma warning (disable: 4127)
23 : eps_(numeric_limits<real>::epsilon())
24 , epsx_(
Math::sq(eps_))
25 , tol_(real(0.1) * sqrt(eps_))
26 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
28 , _f(f <= 1 ? f : 1/f)
40 if (!(abs(stdlat) <= 90))
45 cphi = abs(stdlat) != 90 ? cos(phi) : 0;
46 Init(sphi, cphi, sphi, cphi, k0);
50 real stdlat1, real stdlat2,
52 : eps_(numeric_limits<real>::epsilon())
53 , epsx_(
Math::sq(eps_))
54 , tol_(real(0.1) * sqrt(eps_))
55 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
57 , _f(f <= 1 ? f : 1/f)
69 if (!(abs(stdlat1) <= 90))
70 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
71 if (!(abs(stdlat2) <= 90))
72 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
76 Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
77 sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
81 real sinlat1, real coslat1,
82 real sinlat2, real coslat2,
84 : eps_(numeric_limits<real>::epsilon())
85 , epsx_(
Math::sq(eps_))
86 , tol_(real(0.1) * sqrt(eps_))
87 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
89 , _f(f <= 1 ? f : 1/f)
102 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
104 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
105 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
106 throw GeographicErr(
"Bad sine/cosine of standard latitude 1");
107 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
108 throw GeographicErr(
"Bad sine/cosine of standard latitude 2");
109 if (coslat1 == 0 || coslat2 == 0)
110 if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
112 (
"Standard latitudes must be equal is either is a pole");
113 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
116 void LambertConformalConic::Init(
real sphi1,
real cphi1,
121 sphi1 /= r; cphi1 /= r;
123 sphi2 /= r; cphi2 /= r;
125 bool polar = (cphi1 == 0);
126 cphi1 = max(epsx_, cphi1);
127 cphi2 = max(epsx_, cphi2);
129 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
131 sphi1 *= _sign; sphi2 *= _sign;
133 swap(sphi1, sphi2); swap(cphi1, cphi2);
136 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
156 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
157 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
160 xi1 = eatanhe(sphi1), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
161 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
163 xi2 = eatanhe(sphi2), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
164 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
166 if (tphi2 - tphi1 != 0) {
170 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
172 real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
173 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
177 _nc = sqrt((1 - _n) * (1 + _n));
199 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
201 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
204 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
205 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
206 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
207 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
208 t = Dlog1p(a2, a1) / den;
211 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
212 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
213 (4 * scbet1 * scbet2) ) * _fm;
220 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
221 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
233 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
235 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
236 1 / (scbet1 + _fm * scphi1) );
251 xiZ = eatanhe(
real(1)), shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
254 dxiZ1 = Deatanhe(
real(1), sphi1)/(scphi1*(tphi1+scphi1)),
255 dxiZ2 = Deatanhe(
real(1), sphi2)/(scphi2*(tphi2+scphi2)),
256 dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
257 dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
258 dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
259 dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
261 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
262 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
264 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
267 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
269 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
270 - ( (scphi1 + scphi2)/2
271 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
273 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
274 + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2)
275 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
276 - (dchxiZ1 + dchxiZ2)/2 ),
278 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
279 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
281 _nc = sqrt(max(
real(0), t) * (1 + _n));
297 _scbet0 = hyp(_fm * tphi0);
298 real shxi0 = sinh(eatanhe(_n));
299 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
306 _scale = _a * k1 / scbet1 *
309 exp( - (
Math::sq(_nc)/(1 + _n)) * psi1 )
310 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
314 _k0 = k1 * (_scbet0/scbet1) *
316 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
317 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
319 _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
323 sphi = -1, cphi = epsx_,
325 scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
326 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
328 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0);
329 _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ?
330 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
331 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
332 - (_t0nm1 + 1))/(-_n) :
333 Dexp(-_n * psi, -_n * _psi0) * dpsi);
345 real& x, real& y, real& gamma, real& k)
361 sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx_,
362 tphi = sphi/cphi, scbet = hyp(_fm * tphi),
363 scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
364 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
366 theta = _n * lam, stheta = sin(theta), ctheta = cos(theta),
367 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0),
368 drho = - _scale * (2 * _nc < 1 && dpsi != 0 ?
369 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
370 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
371 - (_t0nm1 + 1))/(-_n) :
372 Dexp(-_n * psi, -_n * _psi0) * dpsi);
373 x = (_nrho0 + _n * drho) * (_n ? stheta / _n : lam);
376 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) / _n : 0)
378 k = _k0 * (scbet/_scbet0) /
379 (exp( - (
Math::sq(_nc)/(1 + _n)) * dpsi )
380 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
386 real& lat, real& lon,
387 real& gamma, real& k)
404 nx = _n * x, ny = _n ? _n * y : 0, y1 = _nrho0 - ny,
408 ? (x*nx + y * (ny - 2*_nrho0)) / den
410 drho = min(drho, _drhomax);
412 drho = max(drho, -_drhomax);
414 tnm1 = _t0nm1 + _n * drho/_scale,
415 dpsi = (den == 0 ? 0 :
416 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
422 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
423 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
424 tchi = _tchi0 + dtchi;
433 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
434 sh = sinh( -
Math::sq(_nc)/(_n * (1 + _n)) *
436 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
443 stol = tol_ * max(real(1), abs(tchi));
448 shxi = sinh( eatanhe( tphi / scphi ) ),
449 tchia = hyp(shxi) * tphi - shxi * scphi,
450 dtphi = (tchi - tchia) * (1 + _e2m *
Math::sq(tphi)) /
451 ( _e2m * scphi * hyp(tchia) );
453 if (!(abs(dtphi) >= stol))
457 gamma = atan2(nx, y1);
459 phi = _sign * atan(tphi),
460 scbet = hyp(_fm * tphi), scchi = hyp(tchi),
461 lam = _n ? gamma / _n : x / y1;
465 k = _k0 * (scbet/_scbet0) /
466 (exp(_nc ? - (
Math::sq(_nc)/(1 + _n)) * dpsi : 0)
467 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
474 if (!(abs(lat) <= 90))
475 throw GeographicErr(
"Latitude for SetScale not in [-90d, 90d]");
476 if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0))
477 throw GeographicErr(
"Incompatible polar latitude in SetScale");
478 real x, y, gamma, kold;
479 Forward(0, lat, 0, x, y, gamma, kold);
static T AngNormalize(T x)
GeographicLib::Math::real real
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
Namespace for GeographicLib.
static T AngDiff(T x, T y)
Exception handling for GeographicLib.
#define GEOGRAPHICLIB_PANIC