00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include "GeographicLib/AlbersEqualArea.hpp"
00011
00012 #define GEOGRAPHICLIB_ALBERSEQUALAREA_CPP "$Id: AlbersEqualArea.cpp 6937 2011-02-01 20:17:13Z karney $"
00013
00014 RCSID_DECL(GEOGRAPHICLIB_ALBERSEQUALAREA_CPP)
00015 RCSID_DECL(GEOGRAPHICLIB_ALBERSEQUALAREA_HPP)
00016
00017 namespace GeographicLib {
00018
00019 using namespace std;
00020
00021 const Math::real AlbersEqualArea::eps = numeric_limits<real>::epsilon();
00022 const Math::real AlbersEqualArea::epsx = sq(eps);
00023 const Math::real AlbersEqualArea::epsx2 = sq(epsx);
00024 const Math::real AlbersEqualArea::tol = sqrt(eps);
00025 const Math::real AlbersEqualArea::tol0 = tol * sqrt(sqrt(eps));
00026 const Math::real AlbersEqualArea::ahypover =
00027 real(numeric_limits<real>::digits) * log(real(numeric_limits<real>::radix))
00028 + 2;
00029
00030 AlbersEqualArea::AlbersEqualArea(real a, real r,
00031 real stdlat, real k0)
00032 : _a(a)
00033 , _r(r)
00034 , _f(_r != 0 ? 1 / _r : 0)
00035 , _fm(1 - _f)
00036 , _e2(_f * (2 - _f))
00037 , _e(sqrt(abs(_e2)))
00038 , _e2m(1 - _e2)
00039 , _qZ(1 + _e2m * atanhee(real(1)))
00040 , _qx(_qZ / ( 2 * _e2m ))
00041 {
00042 if (!(_a > 0))
00043 throw GeographicErr("Major radius is not positive");
00044 if (!(_f < 1))
00045 throw GeographicErr("Minor radius is not positive");
00046 if (!(k0 > 0))
00047 throw GeographicErr("Scale is not positive");
00048 if (!(abs(stdlat) <= 90))
00049 throw GeographicErr("Standard latitude not in [-90, 90]");
00050 real
00051 phi = stdlat * Math::degree<real>(),
00052 sphi = sin(phi),
00053 cphi = abs(stdlat) != 90 ? cos(phi) : 0;
00054 Init(sphi, cphi, sphi, cphi, k0);
00055 }
00056
00057 AlbersEqualArea::AlbersEqualArea(real a, real r,
00058 real stdlat1, real stdlat2,
00059 real k1)
00060 : _a(a)
00061 , _r(r)
00062 , _f(_r != 0 ? 1 / _r : 0)
00063 , _fm(1 - _f)
00064 , _e2(_f * (2 - _f))
00065 , _e(sqrt(abs(_e2)))
00066 , _e2m(1 - _e2)
00067 , _qZ(1 + _e2m * atanhee(real(1)))
00068 , _qx(_qZ / ( 2 * _e2m ))
00069 {
00070 if (!(_a > 0))
00071 throw GeographicErr("Major radius is not positive");
00072 if (!(_f < 1))
00073 throw GeographicErr("Minor radius is not positive");
00074 if (!(k1 > 0))
00075 throw GeographicErr("Scale is not positive");
00076 if (!(abs(stdlat1) <= 90))
00077 throw GeographicErr("Standard latitude 1 not in [-90, 90]");
00078 if (!(abs(stdlat2) <= 90))
00079 throw GeographicErr("Standard latitude 2 not in [-90, 90]");
00080 real
00081 phi1 = stdlat1 * Math::degree<real>(),
00082 phi2 = stdlat2 * Math::degree<real>();
00083 Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
00084 sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
00085 }
00086
00087 AlbersEqualArea::AlbersEqualArea(real a, real r,
00088 real sinlat1, real coslat1,
00089 real sinlat2, real coslat2,
00090 real k1)
00091 : _a(a)
00092 , _r(r)
00093 , _f(_r != 0 ? 1 / _r : 0)
00094 , _fm(1 - _f)
00095 , _e2(_f * (2 - _f))
00096 , _e(sqrt(abs(_e2)))
00097 , _e2m(1 - _e2)
00098 , _qZ(1 + _e2m * atanhee(real(1)))
00099 , _qx(_qZ / ( 2 * _e2m ))
00100 {
00101 if (!(_a > 0))
00102 throw GeographicErr("Major radius is not positive");
00103 if (!(_f < 1))
00104 throw GeographicErr("Minor radius is not positive");
00105 if (!(k1 > 0))
00106 throw GeographicErr("Scale is not positive");
00107 if (!(coslat1 >= 0))
00108 throw GeographicErr("Standard latitude 1 not in [-90, 90]");
00109 if (!(coslat2 >= 0))
00110 throw GeographicErr("Standard latitude 2 not in [-90, 90]");
00111 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
00112 throw GeographicErr("Bad sine/cosine of standard latitude 1");
00113 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
00114 throw GeographicErr("Bad sine/cosine of standard latitude 2");
00115 if (coslat1 == 0 && coslat2 == 0 && sinlat1 * sinlat2 <= 0)
00116 throw GeographicErr
00117 ("Standard latitudes cannot be opposite poles");
00118 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
00119 }
00120
00121 void AlbersEqualArea::Init(real sphi1, real cphi1,
00122 real sphi2, real cphi2, real k1) throw() {
00123 {
00124 real r;
00125 r = Math::hypot(sphi1, cphi1);
00126 sphi1 /= r; cphi1 /= r;
00127 r = Math::hypot(sphi2, cphi2);
00128 sphi2 /= r; cphi2 /= r;
00129 }
00130 bool polar = (cphi1 == 0);
00131 cphi1 = max(epsx, cphi1);
00132 cphi2 = max(epsx, cphi2);
00133
00134 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
00135
00136 sphi1 *= _sign; sphi2 *= _sign;
00137 if (sphi1 > sphi2) {
00138 swap(sphi1, sphi2); swap(cphi1, cphi2);
00139 }
00140 real
00141 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2;
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155
00156
00157
00158
00159
00160
00161
00162
00163
00164
00165 real tphi0, C;
00166 if (polar || tphi1 == tphi2) {
00167 tphi0 = tphi2;
00168 C = 1;
00169 } else {
00170 real
00171 tbet1 = _fm * tphi1, scbet12 = 1 + sq(tbet1),
00172 tbet2 = _fm * tphi2, scbet22 = 1 + sq(tbet2),
00173 txi1 = txif(tphi1), cxi1 = 1/hyp(txi1), sxi1 = txi1 * cxi1,
00174 txi2 = txif(tphi2), cxi2 = 1/hyp(txi2), sxi2 = txi2 * cxi2,
00175 dtbet2 = _fm * (tbet1 + tbet2),
00176 es1 = 1 - _e2 * sq(sphi1), es2 = 1 - _e2 * sq(sphi2),
00177
00178
00179
00180
00181
00182 dsxi = ( (1 + _e2 * sphi1 * sphi2) / (es2 * es1) +
00183 Datanhee(sphi2, sphi1) ) * Dsn(tphi2, tphi1, sphi2, sphi1) /
00184 ( 2 * _qx ),
00185 den = (sxi2 + sxi1) * dtbet2 + (scbet22 + scbet12) * dsxi,
00186
00187 s = 2 * dtbet2 / den,
00188
00189
00190
00191
00192 sm1 = -Dsn(tphi2, tphi1, sphi2, sphi1) *
00193 ( -( ((sphi2 <= 0 ? (1 - sxi2) / (1 - sphi2) :
00194 sq(cxi2/cphi2) * (1 + sphi2) / (1 + sxi2)) +
00195 (sphi1 <= 0 ? (1 - sxi1) / (1 - sphi1) :
00196 sq(cxi1/cphi1) * (1 + sphi1) / (1 + sxi1))) ) *
00197 (1 + _e2 * (sphi1 + sphi2 + sphi1 * sphi2)) /
00198 (1 + (sphi1 + sphi2 + sphi1 * sphi2)) +
00199 (scbet22 * (sphi2 <= 0 ? 1 - sphi2 : sq(cphi2) / ( 1 + sphi2)) +
00200 scbet12 * (sphi1 <= 0 ? 1 - sphi1 : sq(cphi1) / ( 1 + sphi1))) *
00201 (_e2 * (1 + sphi1 + sphi2 + _e2 * sphi1 * sphi2)/(es1 * es2)
00202 +_e2m * DDatanhee(sphi1, sphi2) ) / _qZ ) / den;
00203
00204 C = den / (2 * scbet12 * scbet22 * dsxi);
00205 tphi0 = (tphi2 + tphi1)/2;
00206 real stol = tol0 * max(real(1), abs(tphi0));
00207 for (int i = 0; i < 2*numit0; ++i) {
00208
00209
00210
00211
00212
00213
00214
00215
00216
00217
00218
00219
00220
00221
00222
00223
00224
00225
00226
00227
00228
00229
00230
00231
00232
00233
00234
00235
00236
00237
00238
00239 real
00240 scphi02 = 1 + sq(tphi0), scphi0 = sqrt(scphi02),
00241
00242 sphi0 = tphi0 / scphi0, sphi0m = 1/(scphi0 * (tphi0 + scphi0)),
00243
00244 g = (1 + sq( _fm * tphi0 )) * sphi0,
00245
00246 dg = _e2m * scphi02 * (1 + 2 * sq(tphi0)) + _e2,
00247 D = sphi0m * (1 - _e2*(1 + 2*sphi0*(1+sphi0))) / (_e2m * (1+sphi0)),
00248
00249 dD = -2 * (1 - _e2*sq(sphi0) * (2*sphi0+3)) / (_e2m * sq(1+sphi0)),
00250 A = -_e2 * sq(sphi0m) * (2+(1+_e2)*sphi0) / (_e2m*(1-_e2*sq(sphi0))),
00251 B = (sphi0m * _e2m / (1 - _e2*sphi0) *
00252 (atanhxm1(_e2 * sq(sphi0m / (1-_e2*sphi0))) - _e2*sphi0m/_e2m)),
00253
00254 dAB = (2 * _e2 * (2 - _e2 * (1 + sq(sphi0))) /
00255 (_e2m * sq(1 - _e2*sq(sphi0)) * scphi02)),
00256 u = sm1 * g - s/_qZ * ( D - g * (A + B) ),
00257
00258 du = sm1 * dg - s/_qZ * (dD - dg * (A + B) - g * dAB),
00259 dtu = -u/du * (scphi0 * scphi02);
00260 tphi0 += dtu;
00261 if (!(abs(dtu) >= stol))
00262 break;
00263 }
00264 }
00265 _txi0 = txif(tphi0); _scxi0 = hyp(_txi0); _sxi0 = _txi0 / _scxi0;
00266 _n0 = tphi0/hyp(tphi0);
00267 _m02 = 1 / (1 + sq(_fm * tphi0));
00268 _nrho0 = polar ? 0 : _a * sqrt(_m02);
00269 _k0 = sqrt(tphi1 == tphi2 ? 1 : C / (_m02 + _n0 * _qZ * _sxi0)) * k1;
00270 _k2 = sq(_k0);
00271 _lat0 = _sign * atan(tphi0)/Math::degree<real>();
00272 }
00273
00274 const AlbersEqualArea
00275 AlbersEqualArea::CylindricalEqualArea(Constants::WGS84_a<real>(),
00276 Constants::WGS84_r<real>(),
00277 real(0), real(1), real(0), real(1),
00278 real(1));
00279
00280 const AlbersEqualArea
00281 AlbersEqualArea::AzimuthalEqualAreaNorth(Constants::WGS84_a<real>(),
00282 Constants::WGS84_r<real>(),
00283 real(1), real(0), real(1), real(0),
00284 real(1));
00285
00286 const AlbersEqualArea
00287 AlbersEqualArea::AzimuthalEqualAreaSouth(Constants::WGS84_a<real>(),
00288 Constants::WGS84_r<real>(),
00289 real(-1), real(0), real(-1), real(0),
00290 real(1));
00291
00292 Math::real AlbersEqualArea::txif(real tphi) const throw() {
00293
00294
00295
00296
00297
00298
00299
00300
00301
00302
00303 int s = tphi < 0 ? -1 : 1;
00304 tphi *= s;
00305 real
00306 cphi2 = 1 / (1 + sq(tphi)),
00307 sphi = tphi * sqrt(cphi2),
00308 es1 = _e2 * sphi,
00309 es2m1 = 1 - es1 * sphi,
00310 sp1 = 1 + sphi,
00311 es1m1 = (1 - es1) * sp1,
00312 es2m1a = _e2m * es2m1,
00313 es1p1 = sp1 / (1 + es1);
00314 return s * ( sphi / es2m1 + atanhee(sphi) ) /
00315 sqrt( ( cphi2 / (es1p1 * es2m1a) + atanhee(cphi2 / es1m1) ) *
00316 ( es1m1 / es2m1a + atanhee(es1p1) ) );
00317 }
00318
00319 Math::real AlbersEqualArea::tphif(real txi) const throw() {
00320 real
00321 tphi = txi,
00322 stol = tol * max(real(1), abs(txi));
00323
00324 for (int i = 0; i < numit; ++i) {
00325
00326 real
00327 txia = txif(tphi),
00328 tphi2 = sq(tphi),
00329 scphi2 = 1 + tphi2,
00330 scterm = scphi2/(1 + sq(txia)),
00331 dtphi = (txi - txia) * scterm * sqrt(scterm) *
00332 _qx * sq(1 - _e2 * tphi2 / scphi2);
00333 tphi += dtphi;
00334 if (!(abs(dtphi) >= stol))
00335 break;
00336 }
00337 return tphi;
00338 }
00339
00340
00341
00342 Math::real AlbersEqualArea::atanhxm1(real x) throw() {
00343 real s = 0;
00344 if (abs(x) < real(0.5)) {
00345 real os = -1, y = 1, k = 1;
00346 while (os != s) {
00347 os = s;
00348 y *= x;
00349 k += 2;
00350 s += y/k;
00351 }
00352 } else {
00353 real xs = sqrt(abs(x));
00354 s = (x > 0 ? Math::atanh(xs) : atan(xs)) / xs - 1;
00355 }
00356 return s;
00357 }
00358
00359
00360 Math::real AlbersEqualArea::DDatanhee(real x, real y) const throw() {
00361 real s = 0;
00362 if (_e2 * (abs(x) + abs(y)) < real(0.5)) {
00363 real os = -1, z = 1, k = 1, t = 0, c = 0, en = 1;
00364 while (os != s) {
00365 os = s;
00366 t = y * t + z; c += t; z *= x;
00367 t = y * t + z; c += t; z *= x;
00368 k += 2; en *= _e2;
00369
00370
00371 s += en * c / k;
00372 }
00373
00374
00375 } else
00376 s = (Datanhee(1, y) - Datanhee(x, y))/(1 - x);
00377 return s;
00378 }
00379
00380 void AlbersEqualArea::Forward(real lon0, real lat, real lon,
00381 real& x, real& y, real& gamma, real& k)
00382 const throw() {
00383 if (lon - lon0 >= 180)
00384 lon -= lon0 + 360;
00385 else if (lon - lon0 < -180)
00386 lon -= lon0 - 360;
00387 else
00388 lon -= lon0;
00389 lat *= _sign;
00390 real
00391 lam = lon * Math::degree<real>(),
00392 phi = lat * Math::degree<real>(),
00393 sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx,
00394 tphi = sphi/cphi, txi = txif(tphi), sxi = txi/hyp(txi),
00395 dq = _qZ * Dsn(txi, _txi0, sxi, _sxi0) * (txi - _txi0),
00396 drho = - _a * dq / (sqrt(_m02 - _n0 * dq) + _nrho0 / _a),
00397 theta = _k2 * _n0 * lam, stheta = sin(theta), ctheta = cos(theta),
00398 t = _nrho0 + _n0 * drho;
00399 x = t * (_n0 != 0 ? stheta / _n0 : _k2 * lam) / _k0;
00400 y = (_nrho0 *
00401 (_n0 != 0 ?
00402 (ctheta < 0 ? 1 - ctheta : sq(stheta)/(1 + ctheta)) / _n0 :
00403 0)
00404 - drho * ctheta) / _k0;
00405 k = _k0 * (t != 0 ? t * hyp(_fm * tphi) / _a : 1);
00406 y *= _sign;
00407 gamma = _sign * theta / Math::degree<real>();
00408 }
00409
00410 void AlbersEqualArea::Reverse(real lon0, real x, real y,
00411 real& lat, real& lon,
00412 real& gamma, real& k)
00413 const throw() {
00414 y *= _sign;
00415 real
00416 nx = _k0 * _n0 * x, ny = _k0 * _n0 * y, y1 = _nrho0 - ny,
00417 den = Math::hypot(nx, y1) + _nrho0,
00418 drho = den != 0 ? (_k0*x*nx - 2*_k0*y*_nrho0 + _k0*y*ny) / den : 0,
00419
00420 dsxia = - _scxi0 * (2 * _nrho0 + _n0 * drho) * drho / (sq(_a) * _qZ),
00421 txi = (_txi0 + dsxia) / sqrt(max(1 - dsxia * (2*_txi0 + dsxia), epsx2)),
00422 tphi = tphif(txi),
00423 phi = _sign * atan(tphi),
00424 theta = atan2(nx, y1),
00425 lam = _n0 != 0 ? theta / (_k2 * _n0) : x / (y1 * _k0);
00426 gamma = _sign * theta / Math::degree<real>();
00427 lat = phi / Math::degree<real>();
00428 lon = lam / Math::degree<real>();
00429 k = _k0 * (den != 0 ? (_nrho0 + _n0 * drho) * hyp(_fm * tphi) / _a : 1);
00430 }
00431
00432 void AlbersEqualArea::SetScale(real lat, real k) {
00433 if (!(k > 0))
00434 throw GeographicErr("Scale is not positive");
00435 if (!(abs(lat) < 90))
00436 throw GeographicErr("Latitude for SetScale not in (-90, 90)");
00437 real x, y, gamma, kold;
00438 Forward(0, lat, 0, x, y, gamma, kold);
00439 k /= kold;
00440 _k0 *= k;
00441 _k2 = sq(_k0);
00442 }
00443
00444 }