1 /** 2 * \file LambertConformalConic.cpp 3 * \brief Implementation for GeographicLib::LambertConformalConic class 4 * 5 * Copyright (c) Charles Karney (2010-2020) <charles@karney.com> and licensed 6 * under the MIT/X11 License. For more information, see 7 * https://geographiclib.sourceforge.io/ 8 **********************************************************************/ 9 10 #include <GeographicLib/LambertConformalConic.hpp> 11 12 namespace GeographicLib { 13 14 using namespace std; 15 LambertConformalConic(real a,real f,real stdlat,real k0)16 LambertConformalConic::LambertConformalConic(real a, real f, 17 real stdlat, real k0) 18 : eps_(numeric_limits<real>::epsilon()) 19 , epsx_(Math::sq(eps_)) 20 , ahypover_(Math::digits() * log(real(numeric_limits<real>::radix)) + 2) 21 , _a(a) 22 , _f(f) 23 , _fm(1 - _f) 24 , _e2(_f * (2 - _f)) 25 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2))) 26 { 27 if (!(isfinite(_a) && _a > 0)) 28 throw GeographicErr("Equatorial radius is not positive"); 29 if (!(isfinite(_f) && _f < 1)) 30 throw GeographicErr("Polar semi-axis is not positive"); 31 if (!(isfinite(k0) && k0 > 0)) 32 throw GeographicErr("Scale is not positive"); 33 if (!(abs(stdlat) <= 90)) 34 throw GeographicErr("Standard latitude not in [-90d, 90d]"); 35 real sphi, cphi; 36 Math::sincosd(stdlat, sphi, cphi); 37 Init(sphi, cphi, sphi, cphi, k0); 38 } 39 LambertConformalConic(real a,real f,real stdlat1,real stdlat2,real k1)40 LambertConformalConic::LambertConformalConic(real a, real f, 41 real stdlat1, real stdlat2, 42 real k1) 43 : eps_(numeric_limits<real>::epsilon()) 44 , epsx_(Math::sq(eps_)) 45 , ahypover_(Math::digits() * log(real(numeric_limits<real>::radix)) + 2) 46 , _a(a) 47 , _f(f) 48 , _fm(1 - _f) 49 , _e2(_f * (2 - _f)) 50 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2))) 51 { 52 if (!(isfinite(_a) && _a > 0)) 53 throw GeographicErr("Equatorial radius is not positive"); 54 if (!(isfinite(_f) && _f < 1)) 55 throw GeographicErr("Polar semi-axis is not positive"); 56 if (!(isfinite(k1) && k1 > 0)) 57 throw GeographicErr("Scale is not positive"); 58 if (!(abs(stdlat1) <= 90)) 59 throw GeographicErr("Standard latitude 1 not in [-90d, 90d]"); 60 if (!(abs(stdlat2) <= 90)) 61 throw GeographicErr("Standard latitude 2 not in [-90d, 90d]"); 62 real sphi1, cphi1, sphi2, cphi2; 63 Math::sincosd(stdlat1, sphi1, cphi1); 64 Math::sincosd(stdlat2, sphi2, cphi2); 65 Init(sphi1, cphi1, sphi2, cphi2, k1); 66 } 67 LambertConformalConic(real a,real f,real sinlat1,real coslat1,real sinlat2,real coslat2,real k1)68 LambertConformalConic::LambertConformalConic(real a, real f, 69 real sinlat1, real coslat1, 70 real sinlat2, real coslat2, 71 real k1) 72 : eps_(numeric_limits<real>::epsilon()) 73 , epsx_(Math::sq(eps_)) 74 , ahypover_(Math::digits() * log(real(numeric_limits<real>::radix)) + 2) 75 , _a(a) 76 , _f(f) 77 , _fm(1 - _f) 78 , _e2(_f * (2 - _f)) 79 , _es((_f < 0 ? -1 : 1) * sqrt(abs(_e2))) 80 { 81 if (!(isfinite(_a) && _a > 0)) 82 throw GeographicErr("Equatorial radius is not positive"); 83 if (!(isfinite(_f) && _f < 1)) 84 throw GeographicErr("Polar semi-axis is not positive"); 85 if (!(isfinite(k1) && k1 > 0)) 86 throw GeographicErr("Scale is not positive"); 87 if (!(coslat1 >= 0)) 88 throw GeographicErr("Standard latitude 1 not in [-90d, 90d]"); 89 if (!(coslat2 >= 0)) 90 throw GeographicErr("Standard latitude 2 not in [-90d, 90d]"); 91 if (!(abs(sinlat1) <= 1 && coslat1 <= 1) || (coslat1 == 0 && sinlat1 == 0)) 92 throw GeographicErr("Bad sine/cosine of standard latitude 1"); 93 if (!(abs(sinlat2) <= 1 && coslat2 <= 1) || (coslat2 == 0 && sinlat2 == 0)) 94 throw GeographicErr("Bad sine/cosine of standard latitude 2"); 95 if (coslat1 == 0 || coslat2 == 0) 96 if (!(coslat1 == coslat2 && sinlat1 == sinlat2)) 97 throw GeographicErr 98 ("Standard latitudes must be equal is either is a pole"); 99 Init(sinlat1, coslat1, sinlat2, coslat2, k1); 100 } 101 Init(real sphi1,real cphi1,real sphi2,real cphi2,real k1)102 void LambertConformalConic::Init(real sphi1, real cphi1, 103 real sphi2, real cphi2, real k1) { 104 { 105 real r; 106 r = hypot(sphi1, cphi1); 107 sphi1 /= r; cphi1 /= r; 108 r = hypot(sphi2, cphi2); 109 sphi2 /= r; cphi2 /= r; 110 } 111 bool polar = (cphi1 == 0); 112 cphi1 = max(epsx_, cphi1); // Avoid singularities at poles 113 cphi2 = max(epsx_, cphi2); 114 // Determine hemisphere of tangent latitude 115 _sign = sphi1 + sphi2 >= 0 ? 1 : -1; 116 // Internally work with tangent latitude positive 117 sphi1 *= _sign; sphi2 *= _sign; 118 if (sphi1 > sphi2) { 119 swap(sphi1, sphi2); swap(cphi1, cphi2); // Make phi1 < phi2 120 } 121 real 122 tphi1 = sphi1/cphi1, tphi2 = sphi2/cphi2, tphi0; 123 // 124 // Snyder: 15-8: n = (log(m1) - log(m2))/(log(t1)-log(t2)) 125 // 126 // m = cos(bet) = 1/sec(bet) = 1/sqrt(1+tan(bet)^2) 127 // bet = parametric lat, tan(bet) = (1-f)*tan(phi) 128 // 129 // t = tan(pi/4-chi/2) = 1/(sec(chi) + tan(chi)) = sec(chi) - tan(chi) 130 // log(t) = -asinh(tan(chi)) = -psi 131 // chi = conformal lat 132 // tan(chi) = tan(phi)*cosh(xi) - sinh(xi)*sec(phi) 133 // xi = eatanhe(sin(phi)), eatanhe(x) = e * atanh(e*x) 134 // 135 // n = (log(sec(bet2))-log(sec(bet1)))/(asinh(tan(chi2))-asinh(tan(chi1))) 136 // 137 // Let log(sec(bet)) = b(tphi), asinh(tan(chi)) = c(tphi) 138 // Then n = Db(tphi2, tphi1)/Dc(tphi2, tphi1) 139 // In limit tphi2 -> tphi1, n -> sphi1 140 // 141 real 142 tbet1 = _fm * tphi1, scbet1 = hyp(tbet1), 143 tbet2 = _fm * tphi2, scbet2 = hyp(tbet2); 144 real 145 scphi1 = 1/cphi1, 146 xi1 = Math::eatanhe(sphi1, _es), shxi1 = sinh(xi1), chxi1 = hyp(shxi1), 147 tchi1 = chxi1 * tphi1 - shxi1 * scphi1, scchi1 = hyp(tchi1), 148 scphi2 = 1/cphi2, 149 xi2 = Math::eatanhe(sphi2, _es), shxi2 = sinh(xi2), chxi2 = hyp(shxi2), 150 tchi2 = chxi2 * tphi2 - shxi2 * scphi2, scchi2 = hyp(tchi2), 151 psi1 = asinh(tchi1); 152 if (tphi2 - tphi1 != 0) { 153 // Db(tphi2, tphi1) 154 real num = Dlog1p(Math::sq(tbet2)/(1 + scbet2), 155 Math::sq(tbet1)/(1 + scbet1)) 156 * Dhyp(tbet2, tbet1, scbet2, scbet1) * _fm; 157 // Dc(tphi2, tphi1) 158 real den = Dasinh(tphi2, tphi1, scphi2, scphi1) 159 - Deatanhe(sphi2, sphi1) * Dsn(tphi2, tphi1, sphi2, sphi1); 160 _n = num/den; 161 162 if (_n < 0.25) 163 _nc = sqrt((1 - _n) * (1 + _n)); 164 else { 165 // Compute nc = cos(phi0) = sqrt((1 - n) * (1 + n)), evaluating 1 - n 166 // carefully. First write 167 // 168 // Dc(tphi2, tphi1) * (tphi2 - tphi1) 169 // = log(tchi2 + scchi2) - log(tchi1 + scchi1) 170 // 171 // then den * (1 - n) = 172 // (log((tchi2 + scchi2)/(2*scbet2)) - 173 // log((tchi1 + scchi1)/(2*scbet1))) / (tphi2 - tphi1) 174 // = Dlog1p(a2, a1) * (tchi2+scchi2 + tchi1+scchi1)/(4*scbet1*scbet2) 175 // * fm * Q 176 // 177 // where 178 // a1 = ( (tchi1 - scbet1) + (scchi1 - scbet1) ) / (2 * scbet1) 179 // Q = ((scbet2 + scbet1)/fm)/((scchi2 + scchi1)/D(tchi2, tchi1)) 180 // - (tbet2 + tbet1)/(scbet2 + scbet1) 181 real t; 182 { 183 real 184 // s1 = (scbet1 - scchi1) * (scbet1 + scchi1) 185 s1 = (tphi1 * (2 * shxi1 * chxi1 * scphi1 - _e2 * tphi1) - 186 Math::sq(shxi1) * (1 + 2 * Math::sq(tphi1))), 187 s2 = (tphi2 * (2 * shxi2 * chxi2 * scphi2 - _e2 * tphi2) - 188 Math::sq(shxi2) * (1 + 2 * Math::sq(tphi2))), 189 // t1 = scbet1 - tchi1 190 t1 = tchi1 < 0 ? scbet1 - tchi1 : (s1 + 1)/(scbet1 + tchi1), 191 t2 = tchi2 < 0 ? scbet2 - tchi2 : (s2 + 1)/(scbet2 + tchi2), 192 a2 = -(s2 / (scbet2 + scchi2) + t2) / (2 * scbet2), 193 a1 = -(s1 / (scbet1 + scchi1) + t1) / (2 * scbet1); 194 t = Dlog1p(a2, a1) / den; 195 } 196 // multiply by (tchi2 + scchi2 + tchi1 + scchi1)/(4*scbet1*scbet2) * fm 197 t *= ( ( (tchi2 >= 0 ? scchi2 + tchi2 : 1/(scchi2 - tchi2)) + 198 (tchi1 >= 0 ? scchi1 + tchi1 : 1/(scchi1 - tchi1)) ) / 199 (4 * scbet1 * scbet2) ) * _fm; 200 201 // Rewrite 202 // Q = (1 - (tbet2 + tbet1)/(scbet2 + scbet1)) - 203 // (1 - ((scbet2 + scbet1)/fm)/((scchi2 + scchi1)/D(tchi2, tchi1))) 204 // = tbm - tam 205 // where 206 real tbm = ( ((tbet1 > 0 ? 1/(scbet1+tbet1) : scbet1 - tbet1) + 207 (tbet2 > 0 ? 1/(scbet2+tbet2) : scbet2 - tbet2)) / 208 (scbet1+scbet2) ); 209 210 // tam = (1 - ((scbet2+scbet1)/fm)/((scchi2+scchi1)/D(tchi2, tchi1))) 211 // 212 // Let 213 // (scbet2 + scbet1)/fm = scphi2 + scphi1 + dbet 214 // (scchi2 + scchi1)/D(tchi2, tchi1) = scphi2 + scphi1 + dchi 215 // then 216 // tam = D(tchi2, tchi1) * (dchi - dbet) / (scchi1 + scchi2) 217 real 218 // D(tchi2, tchi1) 219 dtchi = den / Dasinh(tchi2, tchi1, scchi2, scchi1), 220 // (scbet2 + scbet1)/fm - (scphi2 + scphi1) 221 dbet = (_e2/_fm) * ( 1 / (scbet2 + _fm * scphi2) + 222 1 / (scbet1 + _fm * scphi1) ); 223 224 // dchi = (scchi2 + scchi1)/D(tchi2, tchi1) - (scphi2 + scphi1) 225 // Let 226 // tzet = chxiZ * tphi - shxiZ * scphi 227 // tchi = tzet + nu 228 // scchi = sczet + mu 229 // where 230 // xiZ = eatanhe(1), shxiZ = sinh(xiZ), chxiZ = cosh(xiZ) 231 // nu = scphi * (shxiZ - shxi) - tphi * (chxiZ - chxi) 232 // mu = - scphi * (chxiZ - chxi) + tphi * (shxiZ - shxi) 233 // then 234 // dchi = ((mu2 + mu1) - D(nu2, nu1) * (scphi2 + scphi1)) / 235 // D(tchi2, tchi1) 236 real 237 xiZ = Math::eatanhe(real(1), _es), 238 shxiZ = sinh(xiZ), chxiZ = hyp(shxiZ), 239 // These are differences not divided differences 240 // dxiZ1 = xiZ - xi1; dshxiZ1 = shxiZ - shxi; dchxiZ1 = chxiZ - chxi 241 dxiZ1 = Deatanhe(real(1), sphi1)/(scphi1*(tphi1+scphi1)), 242 dxiZ2 = Deatanhe(real(1), sphi2)/(scphi2*(tphi2+scphi2)), 243 dshxiZ1 = Dsinh(xiZ, xi1, shxiZ, shxi1, chxiZ, chxi1) * dxiZ1, 244 dshxiZ2 = Dsinh(xiZ, xi2, shxiZ, shxi2, chxiZ, chxi2) * dxiZ2, 245 dchxiZ1 = Dhyp(shxiZ, shxi1, chxiZ, chxi1) * dshxiZ1, 246 dchxiZ2 = Dhyp(shxiZ, shxi2, chxiZ, chxi2) * dshxiZ2, 247 // mu1 + mu2 248 amu12 = (- scphi1 * dchxiZ1 + tphi1 * dshxiZ1 249 - scphi2 * dchxiZ2 + tphi2 * dshxiZ2), 250 // D(xi2, xi1) 251 dxi = Deatanhe(sphi1, sphi2) * Dsn(tphi2, tphi1, sphi2, sphi1), 252 // D(nu2, nu1) 253 dnu12 = 254 ( (_f * 4 * scphi2 * dshxiZ2 > _f * scphi1 * dshxiZ1 ? 255 // Use divided differences 256 (dshxiZ1 + dshxiZ2)/2 * Dhyp(tphi1, tphi2, scphi1, scphi2) 257 - ( (scphi1 + scphi2)/2 258 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) : 259 // Use ratio of differences 260 (scphi2 * dshxiZ2 - scphi1 * dshxiZ1)/(tphi2 - tphi1)) 261 + ( (tphi1 + tphi2)/2 * Dhyp(shxi1, shxi2, chxi1, chxi2) 262 * Dsinh(xi1, xi2, shxi1, shxi2, chxi1, chxi2) * dxi ) 263 - (dchxiZ1 + dchxiZ2)/2 ), 264 // dtchi * dchi 265 dchia = (amu12 - dnu12 * (scphi2 + scphi1)), 266 tam = (dchia - dtchi * dbet) / (scchi1 + scchi2); 267 t *= tbm - tam; 268 _nc = sqrt(max(real(0), t) * (1 + _n)); 269 } 270 { 271 real r = hypot(_n, _nc); 272 _n /= r; 273 _nc /= r; 274 } 275 tphi0 = _n / _nc; 276 } else { 277 tphi0 = tphi1; 278 _nc = 1/hyp(tphi0); 279 _n = tphi0 * _nc; 280 if (polar) 281 _nc = 0; 282 } 283 284 _scbet0 = hyp(_fm * tphi0); 285 real shxi0 = sinh(Math::eatanhe(_n, _es)); 286 _tchi0 = tphi0 * hyp(shxi0) - shxi0 * hyp(tphi0); _scchi0 = hyp(_tchi0); 287 _psi0 = asinh(_tchi0); 288 289 _lat0 = atan(_sign * tphi0) / Math::degree(); 290 _t0nm1 = expm1(- _n * _psi0); // Snyder's t0^n - 1 291 // a * k1 * m1/t1^n = a * k1 * m2/t2^n = a * k1 * n * (Snyder's F) 292 // = a * k1 / (scbet1 * exp(-n * psi1)) 293 _scale = _a * k1 / scbet1 * 294 // exp(n * psi1) = exp(- (1 - n) * psi1) * exp(psi1) 295 // with (1-n) = nc^2/(1+n) and exp(-psi1) = scchi1 + tchi1 296 exp( - (Math::sq(_nc)/(1 + _n)) * psi1 ) 297 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)); 298 // Scale at phi0 = k0 = k1 * (scbet0*exp(-n*psi0))/(scbet1*exp(-n*psi1)) 299 // = k1 * scbet0/scbet1 * exp(n * (psi1 - psi0)) 300 // psi1 - psi0 = Dasinh(tchi1, tchi0) * (tchi1 - tchi0) 301 _k0 = k1 * (_scbet0/scbet1) * 302 exp( - (Math::sq(_nc)/(1 + _n)) * 303 Dasinh(tchi1, _tchi0, scchi1, _scchi0) * (tchi1 - _tchi0)) 304 * (tchi1 >= 0 ? scchi1 + tchi1 : 1 / (scchi1 - tchi1)) / 305 (_scchi0 + _tchi0); 306 _nrho0 = polar ? 0 : _a * _k0 / _scbet0; 307 { 308 // Figure _drhomax using code at beginning of Forward with lat = -90 309 real 310 sphi = -1, cphi = epsx_, 311 tphi = sphi/cphi, 312 scphi = 1/cphi, shxi = sinh(Math::eatanhe(sphi, _es)), 313 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi), 314 psi = asinh(tchi), 315 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0); 316 _drhomax = - _scale * (2 * _nc < 1 && dpsi != 0 ? 317 (exp(Math::sq(_nc)/(1 + _n) * psi ) * 318 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi)) 319 - (_t0nm1 + 1))/(-_n) : 320 Dexp(-_n * psi, -_n * _psi0) * dpsi); 321 } 322 } 323 Mercator()324 const LambertConformalConic& LambertConformalConic::Mercator() { 325 static const LambertConformalConic mercator(Constants::WGS84_a(), 326 Constants::WGS84_f(), 327 real(0), real(1)); 328 return mercator; 329 } 330 Forward(real lon0,real lat,real lon,real & x,real & y,real & gamma,real & k) const331 void LambertConformalConic::Forward(real lon0, real lat, real lon, 332 real& x, real& y, 333 real& gamma, real& k) const { 334 lon = Math::AngDiff(lon0, lon); 335 // From Snyder, we have 336 // 337 // theta = n * lambda 338 // x = rho * sin(theta) 339 // = (nrho0 + n * drho) * sin(theta)/n 340 // y = rho0 - rho * cos(theta) 341 // = nrho0 * (1-cos(theta))/n - drho * cos(theta) 342 // 343 // where nrho0 = n * rho0, drho = rho - rho0 344 // and drho is evaluated with divided differences 345 real sphi, cphi; 346 Math::sincosd(Math::LatFix(lat) * _sign, sphi, cphi); 347 cphi = max(epsx_, cphi); 348 real 349 lam = lon * Math::degree(), 350 tphi = sphi/cphi, scbet = hyp(_fm * tphi), 351 scphi = 1/cphi, shxi = sinh(Math::eatanhe(sphi, _es)), 352 tchi = hyp(shxi) * tphi - shxi * scphi, scchi = hyp(tchi), 353 psi = asinh(tchi), 354 theta = _n * lam, stheta = sin(theta), ctheta = cos(theta), 355 dpsi = Dasinh(tchi, _tchi0, scchi, _scchi0) * (tchi - _tchi0), 356 drho = - _scale * (2 * _nc < 1 && dpsi != 0 ? 357 (exp(Math::sq(_nc)/(1 + _n) * psi ) * 358 (tchi > 0 ? 1/(scchi + tchi) : (scchi - tchi)) 359 - (_t0nm1 + 1))/(-_n) : 360 Dexp(-_n * psi, -_n * _psi0) * dpsi); 361 x = (_nrho0 + _n * drho) * (_n != 0 ? stheta / _n : lam); 362 y = _nrho0 * 363 (_n != 0 ? 364 (ctheta < 0 ? 1 - ctheta : Math::sq(stheta)/(1 + ctheta)) / _n : 0) 365 - drho * ctheta; 366 k = _k0 * (scbet/_scbet0) / 367 (exp( - (Math::sq(_nc)/(1 + _n)) * dpsi ) 368 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0)); 369 y *= _sign; 370 gamma = _sign * theta / Math::degree(); 371 } 372 Reverse(real lon0,real x,real y,real & lat,real & lon,real & gamma,real & k) const373 void LambertConformalConic::Reverse(real lon0, real x, real y, 374 real& lat, real& lon, 375 real& gamma, real& k) const { 376 // From Snyder, we have 377 // 378 // x = rho * sin(theta) 379 // rho0 - y = rho * cos(theta) 380 // 381 // rho = hypot(x, rho0 - y) 382 // drho = (n*x^2 - 2*y*nrho0 + n*y^2)/(hypot(n*x, nrho0-n*y) + nrho0) 383 // theta = atan2(n*x, nrho0-n*y) 384 // 385 // From drho, obtain t^n-1 386 // psi = -log(t), so 387 // dpsi = - Dlog1p(t^n-1, t0^n-1) * drho / scale 388 y *= _sign; 389 real 390 // Guard against 0 * inf in computation of ny 391 nx = _n * x, ny = _n != 0 ? _n * y : 0, y1 = _nrho0 - ny, 392 den = hypot(nx, y1) + _nrho0, // 0 implies origin with polar aspect 393 // isfinite test is to avoid inf/inf 394 drho = ((den != 0 && isfinite(den)) 395 ? (x*nx + y * (ny - 2*_nrho0)) / den 396 : den); 397 drho = min(drho, _drhomax); 398 if (_n == 0) 399 drho = max(drho, -_drhomax); 400 real 401 tnm1 = _t0nm1 + _n * drho/_scale, 402 dpsi = (den == 0 ? 0 : 403 (tnm1 + 1 != 0 ? - Dlog1p(tnm1, _t0nm1) * drho / _scale : 404 ahypover_)); 405 real tchi; 406 if (2 * _n <= 1) { 407 // tchi = sinh(psi) 408 real 409 psi = _psi0 + dpsi, tchia = sinh(psi), scchi = hyp(tchia), 410 dtchi = Dsinh(psi, _psi0, tchia, _tchi0, scchi, _scchi0) * dpsi; 411 tchi = _tchi0 + dtchi; // Update tchi using divided difference 412 } else { 413 // tchi = sinh(-1/n * log(tn)) 414 // = sinh((1-1/n) * log(tn) - log(tn)) 415 // = + sinh((1-1/n) * log(tn)) * cosh(log(tn)) 416 // - cosh((1-1/n) * log(tn)) * sinh(log(tn)) 417 // (1-1/n) = - nc^2/(n*(1+n)) 418 // cosh(log(tn)) = (tn + 1/tn)/2; sinh(log(tn)) = (tn - 1/tn)/2 419 real 420 tn = tnm1 + 1 == 0 ? epsx_ : tnm1 + 1, 421 sh = sinh( -Math::sq(_nc)/(_n * (1 + _n)) * 422 (2 * tn > 1 ? log1p(tnm1) : log(tn)) ); 423 tchi = sh * (tn + 1/tn)/2 - hyp(sh) * (tnm1 * (tn + 1)/tn)/2; 424 } 425 426 // log(t) = -asinh(tan(chi)) = -psi 427 gamma = atan2(nx, y1); 428 real 429 tphi = Math::tauf(tchi, _es), 430 scbet = hyp(_fm * tphi), scchi = hyp(tchi), 431 lam = _n != 0 ? gamma / _n : x / y1; 432 lat = Math::atand(_sign * tphi); 433 lon = lam / Math::degree(); 434 lon = Math::AngNormalize(lon + Math::AngNormalize(lon0)); 435 k = _k0 * (scbet/_scbet0) / 436 (exp(_nc != 0 ? - (Math::sq(_nc)/(1 + _n)) * dpsi : 0) 437 * (tchi >= 0 ? scchi + tchi : 1 / (scchi - tchi)) / (_scchi0 + _tchi0)); 438 gamma /= _sign * Math::degree(); 439 } 440 SetScale(real lat,real k)441 void LambertConformalConic::SetScale(real lat, real k) { 442 if (!(isfinite(k) && k > 0)) 443 throw GeographicErr("Scale is not positive"); 444 if (!(abs(lat) <= 90)) 445 throw GeographicErr("Latitude for SetScale not in [-90d, 90d]"); 446 if (abs(lat) == 90 && !(_nc == 0 && lat * _n > 0)) 447 throw GeographicErr("Incompatible polar latitude in SetScale"); 448 real x, y, gamma, kold; 449 Forward(0, lat, 0, x, y, gamma, kold); 450 k /= kold; 451 _scale *= k; 452 _k0 *= k; 453 } 454 455 } // namespace GeographicLib 456