14 # pragma warning (disable: 4127) 22 : eps_(numeric_limits<real>::epsilon())
23 , epsx_(
Math::sq(eps_))
24 , epsx2_(
Math::sq(epsx_))
26 , tol0_(tol_ * sqrt(sqrt(eps_)))
28 , _f(f <= 1 ? f : 1/f)
33 , _qZ(1 + _e2m * atanhee(real(1)))
34 , _qx(_qZ / ( 2 * _e2m ))
42 if (!(abs(stdlat) <= 90))
47 cphi = abs(stdlat) != 90 ? cos(phi) : 0;
48 Init(sphi, cphi, sphi, cphi, k0);
53 : eps_(numeric_limits<real>::epsilon())
54 , epsx_(
Math::sq(eps_))
55 , epsx2_(
Math::sq(epsx_))
57 , tol0_(tol_ * sqrt(sqrt(eps_)))
59 , _f(f <= 1 ? f : 1/f)
64 , _qZ(1 + _e2m * atanhee(real(1)))
65 , _qx(_qZ / ( 2 * _e2m ))
73 if (!(abs(stdlat1) <= 90))
74 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
75 if (!(abs(stdlat2) <= 90))
76 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
80 Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
81 sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
85 real sinlat1, real coslat1,
86 real sinlat2, real coslat2,
88 : eps_(numeric_limits<real>::epsilon())
89 , epsx_(
Math::sq(eps_))
90 , epsx2_(
Math::sq(epsx_))
92 , tol0_(tol_ * sqrt(sqrt(eps_)))
94 , _f(f <= 1 ? f : 1/f)
99 , _qZ(1 + _e2m * atanhee(real(1)))
100 , _qx(_qZ / ( 2 * _e2m ))
109 throw GeographicErr(
"Standard latitude 1 not in [-90d, 90d]");
111 throw GeographicErr(
"Standard latitude 2 not in [-90d, 90d]");
112 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
113 throw GeographicErr(
"Bad sine/cosine of standard latitude 1");
114 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
115 throw GeographicErr(
"Bad sine/cosine of standard latitude 2");
116 if (coslat1 == 0 && coslat2 == 0 && sinlat1 * sinlat2 <= 0)
118 (
"Standard latitudes cannot be opposite poles");
119 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
122 void AlbersEqualArea::Init(real sphi1, real cphi1,
123 real sphi2, real cphi2, real k1) {
127 sphi1 /= r; cphi1 /= r;
129 sphi2 /= r; cphi2 /= r;
131 bool polar = (cphi1 == 0);
132 cphi1 = max(epsx_, cphi1);
133 cphi2 = max(epsx_, cphi2);
135 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
137 sphi1 *= _sign; sphi2 *= _sign;
139 swap(sphi1, sphi2); swap(cphi1, cphi2);
142 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2;
167 if (polar || tphi1 == tphi2) {
172 tbet1 = _fm * tphi1, scbet12 = 1 +
Math::sq(tbet1),
173 tbet2 = _fm * tphi2, scbet22 = 1 +
Math::sq(tbet2),
174 txi1 = txif(tphi1), cxi1 = 1/hyp(txi1), sxi1 = txi1 * cxi1,
175 txi2 = txif(tphi2), cxi2 = 1/hyp(txi2), sxi2 = txi2 * cxi2,
176 dtbet2 = _fm * (tbet1 + tbet2),
183 dsxi = ( (1 + _e2 * sphi1 * sphi2) / (es2 * es1) +
184 Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) /
186 den = (sxi2 + sxi1) * dtbet2 + (scbet22 + scbet12) * dsxi,
188 s = 2 * dtbet2 / den,
193 sm1 = -Dsn(tphi2, tphi1, sphi2, sphi1) *
194 ( -( ((sphi2 <= 0 ? (1 - sxi2) / (1 - sphi2) :
195 Math::sq(cxi2/cphi2) * (1 + sphi2) / (1 + sxi2)) +
196 (sphi1 <= 0 ? (1 - sxi1) / (1 - sphi1) :
197 Math::sq(cxi1/cphi1) * (1 + sphi1) / (1 + sxi1))) ) *
198 (1 + _e2 * (sphi1 + sphi2 + sphi1 * sphi2)) /
199 (1 + (sphi1 + sphi2 + sphi1 * sphi2)) +
200 (scbet22 * (sphi2 <= 0 ? 1 - sphi2 :
Math::sq(cphi2) / ( 1 + sphi2)) +
201 scbet12 * (sphi1 <= 0 ? 1 - sphi1 :
Math::sq(cphi1) / ( 1 + sphi1)))
202 * (_e2 * (1 + sphi1 + sphi2 + _e2 * sphi1 * sphi2)/(es1 * es2)
203 +_e2m * DDatanhee(sphi1, sphi2) ) / _qZ ) / den;
205 C = den / (2 * scbet12 * scbet22 * dsxi);
206 tphi0 = (tphi2 + tphi1)/2;
207 real stol = tol0_ * max(real(1), abs(tphi0));
241 scphi02 = 1 +
Math::sq(tphi0), scphi0 = sqrt(scphi02),
243 sphi0 = tphi0 / scphi0, sphi0m = 1/(scphi0 * (tphi0 + scphi0)),
245 g = (1 +
Math::sq( _fm * tphi0 )) * sphi0,
247 dg = _e2m * scphi02 * (1 + 2 *
Math::sq(tphi0)) + _e2,
248 D = sphi0m * (1 - _e2*(1 + 2*sphi0*(1+sphi0))) / (_e2m * (1+sphi0)),
250 dD = -2 * (1 - _e2*
Math::sq(sphi0) * (2*sphi0+3)) /
252 A = -_e2 *
Math::sq(sphi0m) * (2+(1+_e2)*sphi0) /
254 B = (sphi0m * _e2m / (1 - _e2*sphi0) *
256 Math::sq(sphi0m / (1-_e2*sphi0))) - _e2*sphi0m/_e2m)),
258 dAB = (2 * _e2 * (2 - _e2 * (1 +
Math::sq(sphi0))) /
260 u = sm1 * g - s/_qZ * ( D - g * (A + B) ),
262 du = sm1 * dg - s/_qZ * (dD - dg * (A + B) - g * dAB),
263 dtu = -u/du * (scphi0 * scphi02);
265 if (!(abs(dtu) >= stol))
269 _txi0 = txif(tphi0); _scxi0 = hyp(_txi0); _sxi0 = _txi0 / _scxi0;
270 _n0 = tphi0/hyp(tphi0);
271 _m02 = 1 / (1 +
Math::sq(_fm * tphi0));
272 _nrho0 = polar ? 0 : _a * sqrt(_m02);
273 _k0 = sqrt(tphi1 == tphi2 ? 1 : C / (_m02 + _n0 * _qZ * _sxi0)) * k1;
281 real(0), real(1), real(0), real(1), real(1));
282 return cylindricalequalarea;
288 real(1), real(0), real(1), real(0), real(1));
289 return azimuthalequalareanorth;
295 real(-1), real(0), real(-1), real(0), real(1));
296 return azimuthalequalareasouth;
299 Math::real AlbersEqualArea::txif(real tphi)
const {
310 int s = tphi < 0 ? -1 : 1;
314 sphi = tphi * sqrt(cphi2),
316 es2m1 = 1 - es1 * sphi,
318 es1m1 = (1 - es1) * sp1,
319 es2m1a = _e2m * es2m1,
320 es1p1 = sp1 / (1 + es1);
321 return s * ( sphi / es2m1 + atanhee(sphi) ) /
322 sqrt( ( cphi2 / (es1p1 * es2m1a) + atanhee(cphi2 / es1m1) ) *
323 ( es1m1 / es2m1a + atanhee(es1p1) ) );
326 Math::real AlbersEqualArea::tphif(real txi)
const {
329 stol = tol_ * max(real(1), abs(txi));
337 scterm = scphi2/(1 +
Math::sq(txia)),
338 dtphi = (txi - txia) * scterm * sqrt(scterm) *
339 _qx *
Math::sq(1 - _e2 * tphi2 / scphi2);
341 if (!(abs(dtphi) >= stol))
349 Math::real AlbersEqualArea::atanhxm1(real x) {
351 if (abs(x) < real(0.5)) {
352 real os = -1, y = 1, k = 1;
360 real xs = sqrt(abs(x));
367 Math::real AlbersEqualArea::DDatanhee(real x, real y)
const {
369 if (_e2 * (abs(x) + abs(y)) < real(0.5)) {
370 real os = -1, z = 1, k = 1, t = 0, c = 0, en = 1;
373 t = y * t + z; c += t; z *= x;
374 t = y * t + z; c += t; z *= x;
383 s = (Datanhee(1, y) - Datanhee(x, y))/(1 - x);
388 real& x, real& y, real& gamma, real& k)
395 sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx_,
396 tphi = sphi/cphi, txi = txif(tphi), sxi = txi/hyp(txi),
397 dq = _qZ * Dsn(txi, _txi0, sxi, _sxi0) * (txi - _txi0),
398 drho = - _a * dq / (sqrt(_m02 - _n0 * dq) + _nrho0 / _a),
399 theta = _k2 * _n0 * lam, stheta = sin(theta), ctheta = cos(theta),
400 t = _nrho0 + _n0 * drho;
401 x = t * (_n0 ? stheta / _n0 : _k2 * lam) / _k0;
404 (ctheta < 0 ? 1 - ctheta :
Math::sq(stheta)/(1 + ctheta)) / _n0 :
406 - drho * ctheta) / _k0;
407 k = _k0 * (t ? t * hyp(_fm * tphi) / _a : 1);
413 real& lat, real& lon,
414 real& gamma, real& k)
418 nx = _k0 * _n0 * x, ny = _k0 * _n0 * y, y1 = _nrho0 - ny,
420 drho = den ? (_k0*x*nx - 2*_k0*y*_nrho0 + _k0*y*ny) / den : 0,
422 dsxia = - _scxi0 * (2 * _nrho0 + _n0 * drho) * drho /
424 txi = (_txi0 + dsxia) / sqrt(max(1 - dsxia * (2*_txi0 + dsxia), epsx2_)),
426 phi = _sign * atan(tphi),
427 theta = atan2(nx, y1),
428 lam = _n0 ? theta / (_k2 * _n0) : x / (y1 * _k0);
433 k = _k0 * (den ? (_nrho0 + _n0 * drho) * hyp(_fm * tphi) / _a : 1);
439 if (!(abs(lat) < 90))
440 throw GeographicErr(
"Latitude for SetScale not in (-90d, 90d)");
441 real x, y, gamma, kold;
442 Forward(0, lat, 0, x, y, gamma, kold);
static T AngNormalize(T x)
AlbersEqualArea(real a, real f, real stdlat, real k0)
static bool isfinite(T x)
static const AlbersEqualArea & CylindricalEqualArea()
Mathematical functions needed by GeographicLib.
Header for GeographicLib::AlbersEqualArea class.
Albers equal area conic projection.
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
static const AlbersEqualArea & AzimuthalEqualAreaNorth()
Namespace for GeographicLib.
static T AngDiff(T x, T y)
Exception handling for GeographicLib.
static const AlbersEqualArea & AzimuthalEqualAreaSouth()
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
#define GEOGRAPHICLIB_PANIC
void SetScale(real lat, real k=real(1))