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