GeographicLib
1.21
|
00001 /** 00002 * \file LambertConformalConic.cpp 00003 * \brief Implementation for GeographicLib::LambertConformalConic class 00004 * 00005 * Copyright (c) Charles Karney (2010, 2011) <charles@karney.com> and licensed 00006 * under the MIT/X11 License. For more information, see 00007 * http://geographiclib.sourceforge.net/ 00008 **********************************************************************/ 00009 00010 #include <GeographicLib/LambertConformalConic.hpp> 00011 00012 #define GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_CPP \ 00013 "$Id: da8f6ce89092006a26946d671edca1a7836e7ce6 $" 00014 00015 RCSID_DECL(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_CPP) 00016 RCSID_DECL(GEOGRAPHICLIB_LAMBERTCONFORMALCONIC_HPP) 00017 00018 namespace GeographicLib { 00019 00020 using namespace std; 00021 00022 const Math::real LambertConformalConic::eps_ = 00023 numeric_limits<real>::epsilon(); 00024 const Math::real LambertConformalConic::epsx_ = Math::sq(eps_); 00025 const Math::real LambertConformalConic::tol_ = real(0.1) * sqrt(eps_); 00026 const Math::real LambertConformalConic::ahypover_ = 00027 real(numeric_limits<real>::digits) * log(real(numeric_limits<real>::radix)) 00028 + 2; 00029 00030 LambertConformalConic::LambertConformalConic(real a, real f, 00031 real stdlat, real k0) 00032 : _a(a) 00033 , _f(f <= 1 ? f : 1/f) 00034 , _fm(1 - _f) 00035 , _e2(_f * (2 - _f)) 00036 , _e(sqrt(abs(_e2))) 00037 , _e2m(1 - _e2) 00038 { 00039 if (!(Math::isfinite(_a) && _a > 0)) 00040 throw GeographicErr("Major radius is not positive"); 00041 if (!(Math::isfinite(_f) && _f < 1)) 00042 throw GeographicErr("Minor radius is not positive"); 00043 if (!(Math::isfinite(k0) && k0 > 0)) 00044 throw GeographicErr("Scale is not positive"); 00045 if (!(abs(stdlat) <= 90)) 00046 throw GeographicErr("Standard latitude not in [-90, 90]"); 00047 real 00048 phi = stdlat * Math::degree<real>(), 00049 sphi = sin(phi), 00050 cphi = abs(stdlat) != 90 ? cos(phi) : 0; 00051 Init(sphi, cphi, sphi, cphi, k0); 00052 } 00053 00054 LambertConformalConic::LambertConformalConic(real a, real f, 00055 real stdlat1, real stdlat2, 00056 real k1) 00057 : _a(a) 00058 , _f(f <= 1 ? f : 1/f) 00059 , _fm(1 - _f) 00060 , _e2(_f * (2 - _f)) 00061 , _e(sqrt(abs(_e2))) 00062 , _e2m(1 - _e2) 00063 { 00064 if (!(Math::isfinite(_a) && _a > 0)) 00065 throw GeographicErr("Major radius is not positive"); 00066 if (!(Math::isfinite(_f) && _f < 1)) 00067 throw GeographicErr("Minor radius is not positive"); 00068 if (!(Math::isfinite(k1) && 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 f, 00082 real sinlat1, real coslat1, 00083 real sinlat2, real coslat2, 00084 real k1) 00085 : _a(a) 00086 , _f(f <= 1 ? f : 1/f) 00087 , _fm(1 - _f) 00088 , _e2(_f * (2 - _f)) 00089 , _e(sqrt(abs(_e2))) 00090 , _e2m(1 - _e2) 00091 { 00092 if (!(Math::isfinite(_a) && _a > 0)) 00093 throw GeographicErr("Major radius is not positive"); 00094 if (!(Math::isfinite(_f) && _f < 1)) 00095 throw GeographicErr("Minor radius is not positive"); 00096 if (!(Math::isfinite(k1) && k1 > 0)) 00097 throw GeographicErr("Scale is not positive"); 00098 if (!(coslat1 >= 0)) 00099 throw GeographicErr("Standard latitude 1 not in [-90, 90]"); 00100 if (!(coslat2 >= 0)) 00101 throw GeographicErr("Standard latitude 2 not in [-90, 90]"); 00102 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0)) 00103 throw GeographicErr("Bad sine/cosine of standard latitude 1"); 00104 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0)) 00105 throw GeographicErr("Bad sine/cosine of standard latitude 2"); 00106 if (coslat1 == 0 || coslat2 == 0) 00107 if (!(coslat1 == coslat2 && sinlat1 == sinlat2)) 00108 throw GeographicErr 00109 ("Standard latitudes must be equal is either is a pole"); 00110 Init(sinlat1, coslat1, sinlat2, coslat2, k1); 00111 } 00112 00113 void LambertConformalConic::Init(real sphi1, real cphi1, 00114 real sphi2, real cphi2, real k1) throw() { 00115 { 00116 real r; 00117 r = Math::hypot(sphi1, cphi1); 00118 sphi1 /= r; cphi1 /= r; 00119 r = Math::hypot(sphi2, cphi2); 00120 sphi2 /= r; cphi2 /= r; 00121 } 00122 bool polar = (cphi1 == 0); 00123 cphi1 = max(epsx_, cphi1); // Avoid singularities at poles 00124 cphi2 = max(epsx_, cphi2); 00125 // Determine hemisphere of tangent latitude 00126 _sign = sphi1 + sphi2 >= 0 ? 1 : -1; 00127 // Internally work with tangent latitude positive 00128 sphi1 *= _sign; sphi2 *= _sign; 00129 if (sphi1 > sphi2) { 00130 swap(sphi1, sphi2); swap(cphi1, cphi2); // Make phi1 < phi2 00131 } 00132 real 00133 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0; 00134 // 00135 // Snyder: 15-8: n = (log(m1) - log(m2))/(log(t1)-log(t2)) 00136 // 00137 // m = cos(bet) = 1/sec(bet) = 1/sqrt(1+tan(bet)^2) 00138 // bet = parametric lat, tan(bet) = (1-f)*tan(phi) 00139 // 00140 // t = tan(pi/4-chi/2) = 1/(sec(chi) + tan(chi)) = sec(chi) - tan(chi) 00141 // log(t) = -asinh(tan(chi)) = -psi 00142 // chi = conformal lat 00143 // tan(chi) = tan(phi)*cosh(xi) - sinh(xi)*sec(phi) 00144 // xi = eatanhe(sin(phi)), eatanhe(x) = e * atanh(e*x) 00145 // 00146 // n = (log(sec(bet2))-log(sec(bet1)))/(asinh(tan(chi2))-asinh(tan(chi1))) 00147 // 00148 // Let log(sec(bet)) = b(tphi), asinh(tan(chi)) = c(tphi) 00149 // Then n = Db(tphi2, tphi1)/Dc(tphi2, tphi1) 00150 // In limit tphi2 -> tphi1, n -> sphi1 00151 // 00152 real 00153 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1), 00154 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2); 00155 real 00156 scphi1 = 1/cphi1, 00157 xi1 = eatanhe(sphi1), shxi1 = sinh(xi1), chxi1 = hyp(shxi1), 00158 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1), 00159 scphi2 = 1/cphi2, 00160 xi2 = eatanhe(sphi2), shxi2 = sinh(xi2), chxi2 = hyp(shxi2), 00161 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2), 00162 psi1 = Math::asinh(tchi1); 00163 if (tphi2 - tphi1 != 0) { 00164 // Db(tphi2, tphi1) 00165 real num = Dlog1p(Math::sq(tbet2)/(1 + scbet2), 00166 Math::sq(tbet1)/(1 + scbet1)) 00167 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm; 00168 // Dc(tphi2, tphi1) 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 // Compute nc = cos(phi0) = sqrt((1 - n) * (1 + n)), evaluating 1 - n 00177 // carefully. First write 00178 // 00179 // Dc(tphi2, tphi1) * (tphi2 - tphi1) 00180 // = log(tchi2 + scchi2) - log(tchi1 + scchi1) 00181 // 00182 // then den * (1 - n) = 00183 // (log((tchi2 + scchi2)/(2*scbet2)) - log((tchi1 + scchi1)/(2*scbet1))) 00184 // / (tphi2 - tphi1) 00185 // = Dlog1p(a2, a1) * (tchi2+scchi2 + tchi1+scchi1)/(4*scbet1*scbet2) 00186 // * fm * Q 00187 // 00188 // where 00189 // a1 = ( (tchi1 - scbet1) + (scchi1 - scbet1) ) / (2 * scbet1) 00190 // Q = ((scbet2 + scbet1)/fm)/((scchi2 + scchi1)/D(tchi2, tchi1)) 00191 // - (tbet2 + tbet1)/(scbet2 + scbet1) 00192 real t; 00193 { 00194 real 00195 // s1 = (scbet1 - scchi1) * (scbet1 + scchi1) 00196 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) - 00197 Math::sq(shxi1) * (1 + 2 * Math::sq(tphi1))), 00198 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) - 00199 Math::sq(shxi2) * (1 + 2 * Math::sq(tphi2))), 00200 // t1 = scbet1 - tchi1 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 // multiply by (tchi2 + scchi2 + tchi1 + scchi1)/(4*scbet1*scbet2) * fm 00208 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) + 00209 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) / 00210 (4 * scbet1 * scbet2) ) * _fm; 00211 00212 // Rewrite 00213 // Q = (1 - (tbet2 + tbet1)/(scbet2 + scbet1)) - 00214 // (1 - ((scbet2 + scbet1)/fm)/((scchi2 + scchi1)/D(tchi2, tchi1))) 00215 // = tbm - tam 00216 // where 00217 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) + 00218 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) / 00219 (scbet1+scbet2) ); 00220 00221 // tam = (1 - ((scbet2+scbet1)/fm)/((scchi2+scchi1)/D(tchi2, tchi1))) 00222 // 00223 // Let 00224 // (scbet2 + scbet1)/fm = scphi2 + scphi1 + dbet 00225 // (scchi2 + scchi1)/D(tchi2, tchi1) = scphi2 + scphi1 + dchi 00226 // then 00227 // tam = D(tchi2, tchi1) * (dchi - dbet) / (scchi1 + scchi2) 00228 real 00229 // D(tchi2, tchi1) 00230 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1), 00231 // (scbet2 + scbet1)/fm - (scphi2 + scphi1) 00232 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) + 00233 1 / (scbet1 + _fm * scphi1) ); 00234 00235 // dchi = (scchi2 + scchi1)/D(tchi2, tchi1) - (scphi2 + scphi1) 00236 // Let 00237 // tzet = chxiZ * tphi - shxiZ * scphi 00238 // tchi = tzet + nu 00239 // scchi = sczet + mu 00240 // where 00241 // xiZ = eatanhe(1), shxiZ = sinh(xiZ), chxiZ = cosh(xiZ) 00242 // nu = scphi * (shxiZ - shxi) - tphi * (chxiZ - chxi) 00243 // mu = - scphi * (chxiZ - chxi) + tphi * (shxiZ - shxi) 00244 // then 00245 // dchi = ((mu2 + mu1) - D(nu2, nu1) * (scphi2 + scphi1)) / 00246 // D(tchi2, tchi1) 00247 real 00248 xiZ = eatanhe(real(1)), shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ), 00249 // These are differences not divided differences 00250 // dxiZ1 = xiZ - xi1; dshxiZ1 = shxiZ - shxi; dchxiZ1 = chxiZ - chxi 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 // mu1 + mu2 00258 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1 00259 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2), 00260 // D(xi2, xi1) 00261 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1), 00262 // D(nu2, nu1) 00263 dnu12 = 00264 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ? 00265 // Use divided differences 00266 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2) 00267 - ( (scphi1 + scphi2)/2 00268 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) : 00269 // Use ratio of differences 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 // dtchi * dchi 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); // Snyder's t0^n - 1 00301 // a * k1 * m1/t1^n = a * k1 * m2/t2^n = a * k1 * n * (Snyder's F) 00302 // = a * k1 / (scbet1 * exp(-n * psi1)) 00303 _scale = _a * k1 / scbet1 * 00304 // exp(n * psi1) = exp(- (1 - n) * psi1) * exp(psi1) 00305 // with (1-n) = nc^2/(1+n) and exp(-psi1) = scchi1 + tchi1 00306 exp( - (Math::sq(_nc)/(1 + _n)) * psi1 ) 00307 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)); 00308 // Scale at phi0 = k0 = k1 * (scbet0*exp(-n*psi0))/(scbet1*exp(-n*psi1)) 00309 // = k1 * scbet0/scbet1 * exp(n * (psi1 - psi0)) 00310 // psi1 - psi0 = Dasinh(tchi1, tchi0) * (tchi1 - tchi0) 00311 _k0 = k1 * (_scbet0/scbet1) * 00312 exp( - (Math::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_f<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 // From Snyder, we have 00335 // 00336 // theta = n * lambda 00337 // x = rho * sin(theta) 00338 // = (nrho0 + n * drho) * sin(theta)/n 00339 // y = rho0 - rho * cos(theta) 00340 // = nrho0 * (1-cos(theta))/n - drho * cos(theta) 00341 // 00342 // where nrho0 = n * rho0, drho = rho - rho0 00343 // and drho is evaluated with divided differences 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(Math::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 ? 00362 (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n : 0) 00363 - drho * ctheta; 00364 k = _k0 * (scbet/_scbet0) / 00365 (exp( - (Math::sq(_nc)/(1 + _n)) * dpsi ) 00366 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0)); 00367 y *= _sign; 00368 gamma = _sign * theta / Math::degree<real>(); 00369 } 00370 00371 void LambertConformalConic::Reverse(real lon0, real x, real y, 00372 real& lat, real& lon, 00373 real& gamma, real& k) 00374 const throw() { 00375 // From Snyder, we have 00376 // 00377 // x = rho * sin(theta) 00378 // rho0 - y = rho * cos(theta) 00379 // 00380 // rho = hypot(x, rho0 - y) 00381 // drho = (n*x^2 - 2*y*nrho0 + n*y^2)/(hypot(n*x, nrho0-n*y) + nrho0) 00382 // theta = atan2(n*x, nrho0-n*y) 00383 // 00384 // From drho, obtain t^n-1 00385 // psi = -log(t), so 00386 // dpsi = - Dlog1p(t^n-1, t0^n-1) * drho / scale 00387 y *= _sign; 00388 real 00389 nx = _n * x, ny = _n * y, y1 = _nrho0 - ny, 00390 den = Math::hypot(nx, y1) + _nrho0, // 0 implies origin with polar aspect 00391 drho = den != 0 ? (x*nx - 2*y*_nrho0 + y*ny) / den : 0, 00392 tnm1 = _t0nm1 + _n * drho/_scale, 00393 dpsi = (den == 0 ? 0 : 00394 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale : 00395 ahypover_)); 00396 real tchi; 00397 if (2 * _n <= 1) { 00398 // tchi = sinh(psi) 00399 real 00400 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia), 00401 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi; 00402 tchi = _tchi0 + dtchi; // Update tchi using divided difference 00403 } else { 00404 // tchi = sinh(-1/n * log(tn)) 00405 // = sinh((1-1/n) * log(tn) - log(tn)) 00406 // = + sinh((1-1/n) * log(tn)) * cosh(log(tn)) 00407 // - cosh((1-1/n) * log(tn)) * sinh(log(tn)) 00408 // (1-1/n) = - nc^2/(n*(1+n)) 00409 // cosh(log(tn)) = (tn + 1/tn)/2; sinh(log(tn)) = (tn - 1/tn)/2 00410 real 00411 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1, 00412 sh = sinh( -Math::sq(_nc)/(_n * (1 + _n)) * 00413 (2 * tn > 1 ? Math::log1p(tnm1) : log(tn)) ); 00414 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2; 00415 } 00416 00417 // Use Newton's method to solve for tphi 00418 real 00419 // See comment in TransverseMercator.cpp about the initial guess 00420 tphi = tchi/_e2m, 00421 stol = tol_ * max(real(1), abs(tchi)); 00422 // min iterations = 1, max iterations = 2; mean = 1.94 00423 for (int i = 0; i < numit_; ++i) { 00424 real 00425 scphi = hyp(tphi), 00426 shxi = sinh( eatanhe( tphi / scphi ) ), 00427 tchia = hyp(shxi) * tphi - shxi * scphi, 00428 dtphi = (tchi - tchia) * (1 + _e2m * Math::sq(tphi)) / 00429 ( _e2m * scphi * hyp(tchia) ); 00430 tphi += dtphi; 00431 if (!(abs(dtphi) >= stol)) 00432 break; 00433 } 00434 // log(t) = -asinh(tan(chi)) = -psi 00435 gamma = atan2(nx, y1); 00436 real 00437 phi = _sign * atan(tphi), 00438 scbet = hyp(_fm * tphi), scchi = hyp(tchi), 00439 lam = _n != 0 ? gamma / _n : x / y1; 00440 lat = phi / Math::degree<real>(); 00441 lon = lam / Math::degree<real>(); 00442 // Avoid losing a bit of accuracy in lon (assuming lon0 is an integer) 00443 if (lon + lon0 >= 180) 00444 lon += lon0 - 360; 00445 else if (lon + lon0 < -180) 00446 lon += lon0 + 360; 00447 else 00448 lon += lon0; 00449 k = _k0 * (scbet/_scbet0) / 00450 (exp(_nc != 0 ? - (Math::sq(_nc)/(1 + _n)) * dpsi : 0) 00451 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0)); 00452 gamma /= _sign * Math::degree<real>(); 00453 } 00454 00455 void LambertConformalConic::SetScale(real lat, real k) { 00456 if (!(Math::isfinite(k) && k > 0)) 00457 throw GeographicErr("Scale is not positive"); 00458 if (!(abs(lat) <= 90)) 00459 throw GeographicErr("Latitude for SetScale not in [-90, 90]"); 00460 if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0)) 00461 throw GeographicErr("Incompatible polar latitude in SetScale"); 00462 real x, y, gamma, kold; 00463 Forward(0, lat, 0, x, y, gamma, kold); 00464 k /= kold; 00465 _scale *= k; 00466 _k0 *= k; 00467 } 00468 00469 } // namespace GeographicLib