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