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