18 : eps_(numeric_limits<real>::epsilon())
19 , epsx_(
Math::sq(eps_))
20 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
22 , _f(f <= 1 ? f : 1/f)
25 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
33 if (!(abs(stdlat) <= 90))
38 cphi = abs(stdlat) != 90 ? cos(phi) : 0;
39 Init(sphi, cphi, sphi, cphi, k0);
43 real stdlat1, real stdlat2,
45 : eps_(numeric_limits<real>::epsilon())
46 , epsx_(
Math::sq(eps_))
47 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
49 , _f(f <= 1 ? f : 1/f)
52 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
60 if (!(abs(stdlat1) <= 90))
61 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
62 if (!(abs(stdlat2) <= 90))
63 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
67 Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
68 sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
72 real sinlat1, real coslat1,
73 real sinlat2, real coslat2,
75 : eps_(numeric_limits<real>::epsilon())
76 , epsx_(
Math::sq(eps_))
77 , ahypover_(
Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
79 , _f(f <= 1 ? f : 1/f)
82 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2)))
91 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
93 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
94 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
95 throw GeographicErr(
"Bad sine/cosine of standard latitude 1");
96 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
97 throw GeographicErr(
"Bad sine/cosine of standard latitude 2");
98 if (coslat1 == 0 || coslat2 == 0)
99 if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
101 (
"Standard latitudes must be equal is either is a pole");
102 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
105 void LambertConformalConic::Init(
real sphi1,
real cphi1,
110 sphi1 /= r; cphi1 /= r;
112 sphi2 /= r; cphi2 /= r;
114 bool polar = (cphi1 == 0);
115 cphi1 = max(epsx_, cphi1);
116 cphi2 = max(epsx_, cphi2);
118 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
120 sphi1 *= _sign; sphi2 *= _sign;
122 swap(sphi1, sphi2); swap(cphi1, cphi2);
125 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
145 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
146 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
149 xi1 =
Math::eatanhe(sphi1, _es), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
150 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
152 xi2 =
Math::eatanhe(sphi2, _es), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
153 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
155 if (tphi2 - tphi1 != 0) {
159 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
161 real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
162 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
166 _nc = sqrt((1 - _n) * (1 + _n));
188 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
190 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
193 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
194 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
195 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
196 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
197 t = Dlog1p(a2, a1) / den;
200 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
201 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
202 (4 * scbet1 * scbet2) ) * _fm;
209 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
210 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
222 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
224 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
225 1 / (scbet1 + _fm * scphi1) );
241 shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
244 dxiZ1 = Deatanhe(
real(1), sphi1)/(scphi1*(tphi1+scphi1)),
245 dxiZ2 = Deatanhe(
real(1), sphi2)/(scphi2*(tphi2+scphi2)),
246 dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
247 dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
248 dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
249 dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
251 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
252 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
254 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
257 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
259 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
260 - ( (scphi1 + scphi2)/2
261 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
263 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
264 + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2)
265 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
266 - (dchxiZ1 + dchxiZ2)/2 ),
268 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
269 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
271 _nc = sqrt(max(
real(0), t) * (1 + _n));
287 _scbet0 = hyp(_fm * tphi0);
289 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
296 _scale = _a * k1 / scbet1 *
299 exp( - (
Math::sq(_nc)/(1 + _n)) * psi1 )
300 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
304 _k0 = k1 * (_scbet0/scbet1) *
306 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
307 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
309 _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
313 sphi = -1, cphi = epsx_,
316 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
318 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0);
319 _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ?
320 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
321 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
322 - (_t0nm1 + 1))/(-_n) :
323 Dexp(-_n * psi, -_n * _psi0) * dpsi);
335 real& x, real& y, real& gamma, real& k)
351 sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx_,
352 tphi = sphi/cphi, scbet = hyp(_fm * tphi),
354 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
356 theta = _n * lam, stheta = sin(theta), ctheta = cos(theta),
357 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0),
358 drho = - _scale * (2 * _nc < 1 && dpsi != 0 ?
359 (exp(
Math::sq(_nc)/(1 + _n) * psi ) *
360 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
361 - (_t0nm1 + 1))/(-_n) :
362 Dexp(-_n * psi, -_n * _psi0) * dpsi);
363 x = (_nrho0 + _n * drho) * (_n ? stheta / _n : lam);
366 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) / _n : 0)
368 k = _k0 * (scbet/_scbet0) /
369 (exp( - (
Math::sq(_nc)/(1 + _n)) * dpsi )
370 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
376 real& lat, real& lon,
377 real& gamma, real& k)
394 nx = _n * x, ny = _n ? _n * y : 0, y1 = _nrho0 - ny,
398 ? (x*nx + y * (ny - 2*_nrho0)) / den
400 drho = min(drho, _drhomax);
402 drho = max(drho, -_drhomax);
404 tnm1 = _t0nm1 + _n * drho/_scale,
405 dpsi = (den == 0 ? 0 :
406 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
412 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
413 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
414 tchi = _tchi0 + dtchi;
423 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
424 sh = sinh( -
Math::sq(_nc)/(_n * (1 + _n)) *
426 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
430 gamma = atan2(nx, y1);
433 phi = _sign * atan(tphi),
434 scbet = hyp(_fm * tphi), scchi = hyp(tchi),
435 lam = _n ? gamma / _n : x / y1;
439 k = _k0 * (scbet/_scbet0) /
440 (exp(_nc ? - (
Math::sq(_nc)/(1 + _n)) * dpsi : 0)
441 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
448 if (!(abs(lat) <= 90))
449 throw GeographicErr(
"Latitude for SetScale not in [-90d, 90d]");
450 if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0))
451 throw GeographicErr(
"Incompatible polar latitude in SetScale");
452 real x, y, gamma, kold;
453 Forward(0, lat, 0, x, y, gamma, kold);
static T AngNormalize(T x)
GeographicLib::Math::real real
static T eatanhe(T x, T es)
static bool isfinite(T x)
Mathematical functions needed by GeographicLib.
Namespace for GeographicLib.
static T AngDiff(T x, T y)
static T tauf(T taup, T es)
Exception handling for GeographicLib.