00001
00002
00003
00004
00005
00006
00007
00008
00009
00010 #include <GeographicLib/LambertConformalConic.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 LambertConformalConic::LambertConformalConic(real a, real f,
00022 real stdlat, real k0)
00023 : eps_(numeric_limits<real>::epsilon())
00024 , epsx_(Math::sq(eps_))
00025 , tol_(real(0.1) * sqrt(eps_))
00026 , ahypover_(Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
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 {
00034 if (!(Math::isfinite(_a) && _a > 0))
00035 throw GeographicErr("Major radius is not positive");
00036 if (!(Math::isfinite(_f) && _f < 1))
00037 throw GeographicErr("Minor radius is not positive");
00038 if (!(Math::isfinite(k0) && k0 > 0))
00039 throw GeographicErr("Scale is not positive");
00040 if (!(abs(stdlat) <= 90))
00041 throw GeographicErr("Standard latitude not in [-90d, 90d]");
00042 real
00043 phi = stdlat * Math::degree(),
00044 sphi = sin(phi),
00045 cphi = abs(stdlat) != 90 ? cos(phi) : 0;
00046 Init(sphi, cphi, sphi, cphi, k0);
00047 }
00048
00049 LambertConformalConic::LambertConformalConic(real a, real f,
00050 real stdlat1, real stdlat2,
00051 real k1)
00052 : eps_(numeric_limits<real>::epsilon())
00053 , epsx_(Math::sq(eps_))
00054 , tol_(real(0.1) * sqrt(eps_))
00055 , ahypover_(Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
00056 , _a(a)
00057 , _f(f <= 1 ? f : 1/f)
00058 , _fm(1 - _f)
00059 , _e2(_f * (2 - _f))
00060 , _e(sqrt(abs(_e2)))
00061 , _e2m(1 - _e2)
00062 {
00063 if (!(Math::isfinite(_a) && _a > 0))
00064 throw GeographicErr("Major radius is not positive");
00065 if (!(Math::isfinite(_f) && _f < 1))
00066 throw GeographicErr("Minor radius is not positive");
00067 if (!(Math::isfinite(k1) && k1 > 0))
00068 throw GeographicErr("Scale is not positive");
00069 if (!(abs(stdlat1) <= 90))
00070 throw GeographicErr("Standard latitude 1 not in [-90d, 90d]");
00071 if (!(abs(stdlat2) <= 90))
00072 throw GeographicErr("Standard latitude 2 not in [-90d, 90d]");
00073 real
00074 phi1 = stdlat1 * Math::degree(),
00075 phi2 = stdlat2 * Math::degree();
00076 Init(sin(phi1), abs(stdlat1) != 90 ? cos(phi1) : 0,
00077 sin(phi2), abs(stdlat2) != 90 ? cos(phi2) : 0, k1);
00078 }
00079
00080 LambertConformalConic::LambertConformalConic(real a, real f,
00081 real sinlat1, real coslat1,
00082 real sinlat2, real coslat2,
00083 real k1)
00084 : eps_(numeric_limits<real>::epsilon())
00085 , epsx_(Math::sq(eps_))
00086 , tol_(real(0.1) * sqrt(eps_))
00087 , ahypover_(Math::digits() * log(real(numeric_limits<real>::radix)) + 2)
00088 , _a(a)
00089 , _f(f <= 1 ? f : 1/f)
00090 , _fm(1 - _f)
00091 , _e2(_f * (2 - _f))
00092 , _e(sqrt(abs(_e2)))
00093 , _e2m(1 - _e2)
00094 {
00095 if (!(Math::isfinite(_a) && _a > 0))
00096 throw GeographicErr("Major radius is not positive");
00097 if (!(Math::isfinite(_f) && _f < 1))
00098 throw GeographicErr("Minor radius is not positive");
00099 if (!(Math::isfinite(k1) && k1 > 0))
00100 throw GeographicErr("Scale is not positive");
00101 if (!(coslat1 >= 0))
00102 throw GeographicErr("Standard latitude 1 not in [-90d, 90d]");
00103 if (!(coslat2 >= 0))
00104 throw GeographicErr("Standard latitude 2 not in [-90d, 90d]");
00105 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0))
00106 throw GeographicErr("Bad sine/cosine of standard latitude 1");
00107 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0))
00108 throw GeographicErr("Bad sine/cosine of standard latitude 2");
00109 if (coslat1 == 0 || coslat2 == 0)
00110 if (!(coslat1 == coslat2 && sinlat1 == sinlat2))
00111 throw GeographicErr
00112 ("Standard latitudes must be equal is either is a pole");
00113 Init(sinlat1, coslat1, sinlat2, coslat2, k1);
00114 }
00115
00116 void LambertConformalConic::Init(real sphi1, real cphi1,
00117 real sphi2, real cphi2, real k1) {
00118 {
00119 real r;
00120 r = Math::hypot(sphi1, cphi1);
00121 sphi1 /= r; cphi1 /= r;
00122 r = Math::hypot(sphi2, cphi2);
00123 sphi2 /= r; cphi2 /= r;
00124 }
00125 bool polar = (cphi1 == 0);
00126 cphi1 = max(epsx_, cphi1);
00127 cphi2 = max(epsx_, cphi2);
00128
00129 _sign = sphi1 + sphi2 >= 0 ? 1 : -1;
00130
00131 sphi1 *= _sign; sphi2 *= _sign;
00132 if (sphi1 > sphi2) {
00133 swap(sphi1, sphi2); swap(cphi1, cphi2);
00134 }
00135 real
00136 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0;
00137
00138
00139
00140
00141
00142
00143
00144
00145
00146
00147
00148
00149
00150
00151
00152
00153
00154
00155 real
00156 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1),
00157 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2);
00158 real
00159 scphi1 = 1/cphi1,
00160 xi1 = eatanhe(sphi1), shxi1 = sinh(xi1), chxi1 = hyp(shxi1),
00161 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1),
00162 scphi2 = 1/cphi2,
00163 xi2 = eatanhe(sphi2), shxi2 = sinh(xi2), chxi2 = hyp(shxi2),
00164 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2),
00165 psi1 = Math::asinh(tchi1);
00166 if (tphi2 - tphi1 != 0) {
00167
00168 real num = Dlog1p(Math::sq(tbet2)/(1 + scbet2),
00169 Math::sq(tbet1)/(1 + scbet1))
00170 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm;
00171
00172 real den = Dasinh(tphi2, tphi1, scphi2, scphi1)
00173 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1);
00174 _n = num/den;
00175
00176 if (_n < 0.25)
00177 _nc = sqrt((1 - _n) * (1 + _n));
00178 else {
00179
00180
00181
00182
00183
00184
00185
00186
00187
00188
00189
00190
00191
00192
00193
00194
00195 real t;
00196 {
00197 real
00198
00199 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) -
00200 Math::sq(shxi1) * (1 + 2 * Math::sq(tphi1))),
00201 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) -
00202 Math::sq(shxi2) * (1 + 2 * Math::sq(tphi2))),
00203
00204 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1),
00205 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2),
00206 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2),
00207 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1);
00208 t = Dlog1p(a2, a1) / den;
00209 }
00210
00211 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) +
00212 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) /
00213 (4 * scbet1 * scbet2) ) * _fm;
00214
00215
00216
00217
00218
00219
00220 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) +
00221 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) /
00222 (scbet1+scbet2) );
00223
00224
00225
00226
00227
00228
00229
00230
00231 real
00232
00233 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1),
00234
00235 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) +
00236 1 / (scbet1 + _fm * scphi1) );
00237
00238
00239
00240
00241
00242
00243
00244
00245
00246
00247
00248
00249
00250 real
00251 xiZ = eatanhe(real(1)), shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ),
00252
00253
00254 dxiZ1 = Deatanhe(real(1), sphi1)/(scphi1*(tphi1+scphi1)),
00255 dxiZ2 = Deatanhe(real(1), sphi2)/(scphi2*(tphi2+scphi2)),
00256 dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1,
00257 dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2,
00258 dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1,
00259 dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2,
00260
00261 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1
00262 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2),
00263
00264 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1),
00265
00266 dnu12 =
00267 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ?
00268
00269 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2)
00270 - ( (scphi1 + scphi2)/2
00271 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) :
00272
00273 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1))
00274 + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2)
00275 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi )
00276 - (dchxiZ1 + dchxiZ2)/2 ),
00277
00278 dchia = (amu12 - dnu12 * (scphi2 + scphi1)),
00279 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2);
00280 t *= tbm - tam;
00281 _nc = sqrt(max(real(0), t) * (1 + _n));
00282 }
00283 {
00284 real r = Math::hypot(_n, _nc);
00285 _n /= r;
00286 _nc /= r;
00287 }
00288 tphi0 = _n / _nc;
00289 } else {
00290 tphi0 = tphi1;
00291 _nc = 1/hyp(tphi0);
00292 _n = tphi0 * _nc;
00293 if (polar)
00294 _nc = 0;
00295 }
00296
00297 _scbet0 = hyp(_fm * tphi0);
00298 real shxi0 = sinh(eatanhe(_n));
00299 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0);
00300 _psi0 = Math::asinh(_tchi0);
00301
00302 _lat0 = atan(_sign * tphi0) / Math::degree();
00303 _t0nm1 = Math::expm1(- _n * _psi0);
00304
00305
00306 _scale = _a * k1 / scbet1 *
00307
00308
00309 exp( - (Math::sq(_nc)/(1 + _n)) * psi1 )
00310 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1));
00311
00312
00313
00314 _k0 = k1 * (_scbet0/scbet1) *
00315 exp( - (Math::sq(_nc)/(1 + _n)) *
00316 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0))
00317 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) /
00318 (_scchi0 + _tchi0);
00319 _nrho0 = polar ? 0 : _a * _k0 / _scbet0;
00320 {
00321
00322 real
00323 sphi = -1, cphi = epsx_,
00324 tphi = sphi/cphi,
00325 scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
00326 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
00327 psi = Math::asinh(tchi),
00328 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0);
00329 _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ?
00330 (exp(Math::sq(_nc)/(1 + _n) * psi ) *
00331 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
00332 - (_t0nm1 + 1))/(-_n) :
00333 Dexp(-_n * psi, -_n * _psi0) * dpsi);
00334 }
00335 }
00336
00337 const LambertConformalConic& LambertConformalConic::Mercator() {
00338 static const LambertConformalConic mercator(Constants::WGS84_a(),
00339 Constants::WGS84_f(),
00340 real(0), real(1));
00341 return mercator;
00342 }
00343
00344 void LambertConformalConic::Forward(real lon0, real lat, real lon,
00345 real& x, real& y, real& gamma, real& k)
00346 const {
00347 lon = Math::AngDiff(Math::AngNormalize(lon0), Math::AngNormalize(lon));
00348
00349
00350
00351
00352
00353
00354
00355
00356
00357
00358 real
00359 lam = lon * Math::degree(),
00360 phi = _sign * lat * Math::degree(),
00361 sphi = sin(phi), cphi = abs(lat) != 90 ? cos(phi) : epsx_,
00362 tphi = sphi/cphi, scbet = hyp(_fm * tphi),
00363 scphi = 1/cphi, shxi = sinh(eatanhe(sphi)),
00364 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi),
00365 psi = Math::asinh(tchi),
00366 theta = _n * lam, stheta = sin(theta), ctheta = cos(theta),
00367 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0),
00368 drho = - _scale * (2 * _nc < 1 && dpsi != 0 ?
00369 (exp(Math::sq(_nc)/(1 + _n) * psi ) *
00370 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi))
00371 - (_t0nm1 + 1))/(-_n) :
00372 Dexp(-_n * psi, -_n * _psi0) * dpsi);
00373 x = (_nrho0 + _n * drho) * (_n ? stheta / _n : lam);
00374 y = _nrho0 *
00375 (_n ?
00376 (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n : 0)
00377 - drho * ctheta;
00378 k = _k0 * (scbet/_scbet0) /
00379 (exp( - (Math::sq(_nc)/(1 + _n)) * dpsi )
00380 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
00381 y *= _sign;
00382 gamma = _sign * theta / Math::degree();
00383 }
00384
00385 void LambertConformalConic::Reverse(real lon0, real x, real y,
00386 real& lat, real& lon,
00387 real& gamma, real& k)
00388 const {
00389
00390
00391
00392
00393
00394
00395
00396
00397
00398
00399
00400
00401 y *= _sign;
00402 real
00403
00404 nx = _n * x, ny = _n ? _n * y : 0, y1 = _nrho0 - ny,
00405 den = Math::hypot(nx, y1) + _nrho0,
00406
00407 drho = ((den != 0 && Math::isfinite(den))
00408 ? (x*nx + y * (ny - 2*_nrho0)) / den
00409 : den);
00410 drho = min(drho, _drhomax);
00411 if (_n == 0)
00412 drho = max(drho, -_drhomax);
00413 real
00414 tnm1 = _t0nm1 + _n * drho/_scale,
00415 dpsi = (den == 0 ? 0 :
00416 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale :
00417 ahypover_));
00418 real tchi;
00419 if (2 * _n <= 1) {
00420
00421 real
00422 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia),
00423 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi;
00424 tchi = _tchi0 + dtchi;
00425 } else {
00426
00427
00428
00429
00430
00431
00432 real
00433 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1,
00434 sh = sinh( -Math::sq(_nc)/(_n * (1 + _n)) *
00435 (2 * tn > 1 ? Math::log1p(tnm1) : log(tn)) );
00436 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2;
00437 }
00438
00439
00440 real
00441
00442 tphi = tchi/_e2m,
00443 stol = tol_ * max(real(1), abs(tchi));
00444
00445 for (int i = 0; i < numit_ || GEOGRAPHICLIB_PANIC; ++i) {
00446 real
00447 scphi = hyp(tphi),
00448 shxi = sinh( eatanhe( tphi / scphi ) ),
00449 tchia = hyp(shxi) * tphi - shxi * scphi,
00450 dtphi = (tchi - tchia) * (1 + _e2m * Math::sq(tphi)) /
00451 ( _e2m * scphi * hyp(tchia) );
00452 tphi += dtphi;
00453 if (!(abs(dtphi) >= stol))
00454 break;
00455 }
00456
00457 gamma = atan2(nx, y1);
00458 real
00459 phi = _sign * atan(tphi),
00460 scbet = hyp(_fm * tphi), scchi = hyp(tchi),
00461 lam = _n ? gamma / _n : x / y1;
00462 lat = phi / Math::degree();
00463 lon = lam / Math::degree();
00464 lon = Math::AngNormalize(lon + Math::AngNormalize(lon0));
00465 k = _k0 * (scbet/_scbet0) /
00466 (exp(_nc ? - (Math::sq(_nc)/(1 + _n)) * dpsi : 0)
00467 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0));
00468 gamma /= _sign * Math::degree();
00469 }
00470
00471 void LambertConformalConic::SetScale(real lat, real k) {
00472 if (!(Math::isfinite(k) && k > 0))
00473 throw GeographicErr("Scale is not positive");
00474 if (!(abs(lat) <= 90))
00475 throw GeographicErr("Latitude for SetScale not in [-90d, 90d]");
00476 if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0))
00477 throw GeographicErr("Incompatible polar latitude in SetScale");
00478 real x, y, gamma, kold;
00479 Forward(0, lat, 0, x, y, gamma, kold);
00480 k /= kold;
00481 _scale *= k;
00482 _k0 *= k;
00483 }
00484
00485 }