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(abs(_e2)))
33 if (!(abs(stdlat) <= 90))
37 Init(sphi, cphi, sphi, cphi, k0);
41 real stdlat1, real stdlat2,
43 : eps_(numeric_limits<real>::epsilon())
44 , epsx_(
Math::sq(eps_))
45 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
50 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
58 if (!(abs(stdlat1) <= 90))
59 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
60 if (!(abs(stdlat2) <= 90))
61 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
62 real sphi1, cphi1, sphi2, cphi2;
65 Init(sphi1, cphi1, sphi2, cphi2, k1);
69 real sinlat1, real coslat1,
70 real sinlat2, real coslat2,
72 : eps_(numeric_limits<real>::epsilon())
73 , epsx_(
Math::sq(eps_))
74 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
79 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
88 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
90 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
91 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
92 throw GeographicErr(
"Bad sine/cosine of standard latitude 1");
93 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
94 throw GeographicErr(
"Bad sine/cosine of standard latitude 2");
95 if (coslat1 == 0 || coslat2 == 0)
96 if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
98 (
"Standard latitudes must be equal is either is a pole");
99 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
102 void LambertConformalConic::Init(
real sphi1,
real cphi1,
107 sphi1 /= r; cphi1 /= r;
109 sphi2 /= r; cphi2 /= r;
111 bool polar = (cphi1 == 0);
112 cphi1 = max(epsx_, cphi1);
113 cphi2 = max(epsx_, cphi2);
115 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
117 sphi1 *= _sign; sphi2 *= _sign;
119 swap(sphi1, sphi2);
swap(cphi1, cphi2);
122 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
142 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
143 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
146 xi1 =
Math::eatanhe(sphi1, _es), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
147 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
149 xi2 =
Math::eatanhe(sphi2, _es), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
150 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
152 if (tphi2 - tphi1 != 0) {
156 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
158 real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
159 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
163 _nc = sqrt((1 - _n) * (1 + _n));
185 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
187 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
190 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
191 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
192 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
193 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
194 t = Dlog1p(a2, a1) / den;
197 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
198 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
199 (4 * scbet1 * scbet2) ) * _fm;
206 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
207 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
219 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
221 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
222 1 / (scbet1 + _fm * scphi1) );
238 shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
241 dxiZ1 = Deatanhe(
real(1), sphi1)/(scphi1*(tphi1+scphi1)),
242 dxiZ2 = Deatanhe(
real(1), sphi2)/(scphi2*(tphi2+scphi2)),
243 dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
244 dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
245 dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
246 dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
248 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
249 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
251 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
254 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
256 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
257 - ( (scphi1 + scphi2)/2
258 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
260 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
261 + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2)
262 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
263 - (dchxiZ1 + dchxiZ2)/2 ),
265 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
266 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
268 _nc = sqrt(max(
real(0), t) * (1 + _n));
284 _scbet0 = hyp(_fm * tphi0);
286 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
293 _scale = _a * k1 / scbet1 *
296 exp( - (
Math::sq(_nc)/(1 + _n)) * psi1 )
297 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
301 _k0 = k1 * (_scbet0/scbet1) *
303 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
304 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
306 _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
310 sphi = -1, cphi = epsx_,
313 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
315 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0);
316 _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ?
317 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
318 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
319 - (_t0nm1 + 1))/(-_n) :
320 Dexp(-_n * psi, -_n * _psi0) * dpsi);
333 real& gamma, real& k)
const {
347 cphi = max(epsx_, cphi);
350 tphi = sphi/cphi, scbet = hyp(_fm * tphi),
352 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
354 theta = _n * lam, stheta = sin(theta), ctheta = cos(theta),
355 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0),
356 drho = - _scale * (2 * _nc < 1 && dpsi != 0 ?
357 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
358 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
359 - (_t0nm1 + 1))/(-_n) :
360 Dexp(-_n * psi, -_n * _psi0) * dpsi);
361 x = (_nrho0 + _n * drho) * (_n != 0 ? stheta / _n : lam);
364 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) / _n : 0)
366 k = _k0 * (scbet/_scbet0) /
367 (exp( - (
Math::sq(_nc)/(1 + _n)) * dpsi )
368 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
374 real& lat, real& lon,
375 real& gamma, real& k)
const {
391 nx = _n * x, ny = _n != 0 ? _n * y : 0, y1 = _nrho0 - ny,
395 ? (x*nx + y * (ny - 2*_nrho0)) / den
397 drho = min(drho, _drhomax);
399 drho = max(drho, -_drhomax);
401 tnm1 = _t0nm1 + _n * drho/_scale,
402 dpsi = (den == 0 ? 0 :
403 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
409 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
410 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
411 tchi = _tchi0 + dtchi;
420 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
421 sh = sinh( -
Math::sq(_nc)/(_n * (1 + _n)) *
423 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
427 gamma = atan2(nx, y1);
430 scbet = hyp(_fm * tphi), scchi = hyp(tchi),
431 lam = _n != 0 ? gamma / _n : x / y1;
435 k = _k0 * (scbet/_scbet0) /
436 (exp(_nc != 0 ? - (
Math::sq(_nc)/(1 + _n)) * dpsi : 0)
437 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
444 if (!(abs(lat) <= 90))
445 throw GeographicErr(
"Latitude for SetScale not in [-90d, 90d]");
446 if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0))
447 throw GeographicErr(
"Incompatible polar latitude in SetScale");
448 real x, y, gamma, kold;
449 Forward(0, lat, 0, x, y, gamma, kold);
static T AngNormalize(T x)
static bool isfinite(T x)
GeographicLib::Math::real real
Mathematical functions needed by GeographicLib.
static T AngDiff(T x, T y, T &e)
Namespace for GeographicLib.
void swap(GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &a, GeographicLib::NearestNeighbor< dist_t, pos_t, distfun_t > &b)
Exception handling for GeographicLib.
static T tauf(T taup, T es)
static void sincosd(T x, T &sinx, T &cosx)
static T eatanhe(T x, T es)