GeographicLib  1.40
 All Classes Namespaces Files Functions Variables Typedefs Enumerations Enumerator Friends Macros Pages
TransverseMercator.cpp
Go to the documentation of this file.
1 /**
2  * \file TransverseMercator.cpp
3  * \brief Implementation for GeographicLib::TransverseMercator class
4  *
5  * Copyright (c) Charles Karney (2008-2014) <charles@karney.com> and licensed
6  * under the MIT/X11 License. For more information, see
7  * http://geographiclib.sourceforge.net/
8  *
9  * This implementation follows closely
10  * <a href="http://www.jhs-suositukset.fi/suomi/jhs154"> JHS 154, ETRS89 -
11  * j&auml;rjestelm&auml;&auml;n liittyv&auml;t karttaprojektiot,
12  * tasokoordinaatistot ja karttalehtijako</a> (Map projections, plane
13  * coordinates, and map sheet index for ETRS89), published by JUHTA, Finnish
14  * Geodetic Institute, and the National Land Survey of Finland (2006).
15  *
16  * The relevant section is available as the 2008 PDF file
17  * http://docs.jhs-suositukset.fi/jhs-suositukset/JHS154/JHS154_liite1.pdf
18  *
19  * This is a straight transcription of the formulas in this paper with the
20  * following exceptions:
21  * - use of 6th order series instead of 4th order series. This reduces the
22  * error to about 5nm for the UTM range of coordinates (instead of 200nm),
23  * with a speed penalty of only 1%;
24  * - use Newton's method instead of plain iteration to solve for latitude in
25  * terms of isometric latitude in the Reverse method;
26  * - use of Horner's representation for evaluating polynomials and Clenshaw's
27  * method for summing trigonometric series;
28  * - several modifications of the formulas to improve the numerical accuracy;
29  * - evaluating the convergence and scale using the expression for the
30  * projection or its inverse.
31  *
32  * If the preprocessor variable GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER is set
33  * to an integer between 4 and 8, then this specifies the order of the series
34  * used for the forward and reverse transformations. The default value is 6.
35  * (The series accurate to 12th order is given in \ref tmseries.)
36  *
37  * Other equivalent implementations are given in
38  * - http://www.ign.fr/DISPLAY/000/526/702/5267021/NTG_76.pdf
39  * - http://www.lantmateriet.se/upload/filer/kartor/geodesi_gps_och_detaljmatning/geodesi/Formelsamling/Gauss_Conformal_Projection.pdf
40  **********************************************************************/
41 
43 
44 #if defined(_MSC_VER)
45 // Squelch warnings about constant conditional expressions
46 # pragma warning (disable: 4127)
47 #endif
48 
49 namespace GeographicLib {
50 
51  using namespace std;
52 
53  TransverseMercator::TransverseMercator(real a, real f, real k0)
54  : tol_(real(0.1)*sqrt(numeric_limits<real>::epsilon()))
55  , _a(a)
56  , _f(f <= 1 ? f : 1/f)
57  , _k0(k0)
58  , _e2(_f * (2 - _f))
59  , _e(sqrt(abs(_e2)))
60  , _e2m(1 - _e2)
61  // _c = sqrt( pow(1 + _e, 1 + _e) * pow(1 - _e, 1 - _e) ) )
62  // See, for example, Lee (1976), p 100.
63  , _c( sqrt(_e2m) * exp(eatanhe(real(1))) )
64  , _n(_f / (2 - _f))
65  {
66  if (!(Math::isfinite(_a) && _a > 0))
67  throw GeographicErr("Major radius is not positive");
68  if (!(Math::isfinite(_f) && _f < 1))
69  throw GeographicErr("Minor radius is not positive");
70  if (!(Math::isfinite(_k0) && _k0 > 0))
71  throw GeographicErr("Scale is not positive");
72  // If coefficents might overflow_ an int, convert them to double (and they
73  // are all exactly representable as doubles).
74  real nx = Math::sq(_n);
75  switch (maxpow_) {
76  case 4:
77  _b1 = 1/(1+_n)*(nx*(nx+16)+64)/64;
78  _alp[1] = _n*(_n*(_n*(164*_n+225)-480)+360)/720;
79  _bet[1] = _n*(_n*((555-4*_n)*_n-960)+720)/1440;
80  _alp[2] = nx*(_n*(557*_n-864)+390)/1440;
81  _bet[2] = nx*((96-437*_n)*_n+30)/1440;
82  nx *= _n;
83  _alp[3] = (427-1236*_n)*nx/1680;
84  _bet[3] = (119-148*_n)*nx/3360;
85  nx *= _n;
86  _alp[4] = 49561*nx/161280;
87  _bet[4] = 4397*nx/161280;
88  break;
89  case 5:
90  _b1 = 1/(1+_n)*(nx*(nx+16)+64)/64;
91  _alp[1] = _n*(_n*(_n*((328-635*_n)*_n+450)-960)+720)/1440;
92  _bet[1] = _n*(_n*(_n*((-3645*_n-64)*_n+8880)-15360)+11520)/23040;
93  _alp[2] = nx*(_n*(_n*(4496*_n+3899)-6048)+2730)/10080;
94  _bet[2] = nx*(_n*(_n*(4416*_n-3059)+672)+210)/10080;
95  nx *= _n;
96  _alp[3] = nx*(_n*(15061*_n-19776)+6832)/26880;
97  _bet[3] = nx*((-627*_n-592)*_n+476)/13440;
98  nx *= _n;
99  _alp[4] = (49561-171840*_n)*nx/161280;
100  _bet[4] = (4397-3520*_n)*nx/161280;
101  nx *= _n;
102  _alp[5] = 34729*nx/80640;
103  _bet[5] = 4583*nx/161280;
104  break;
105  case 6:
106  _b1 = 1/(1+_n)*(nx*(nx*(nx+4)+64)+256)/256;
107  _alp[1] = _n*(_n*(_n*(_n*(_n*(31564*_n-66675)+34440)+47250)-100800)+
108  75600)/151200;
109  _bet[1] = _n*(_n*(_n*(_n*(_n*(384796*_n-382725)-6720)+932400)-1612800)+
110  1209600)/2419200;
111  _alp[2] = nx*(_n*(_n*((863232-1983433*_n)*_n+748608)-1161216)+524160)/
112  1935360;
113  _bet[2] = nx*(_n*(_n*((1695744-1118711*_n)*_n-1174656)+258048)+80640)/
114  3870720;
115  nx *= _n;
116  _alp[3] = nx*(_n*(_n*(670412*_n+406647)-533952)+184464)/725760;
117  _bet[3] = nx*(_n*(_n*(22276*_n-16929)-15984)+12852)/362880;
118  nx *= _n;
119  _alp[4] = nx*(_n*(6601661*_n-7732800)+2230245)/7257600;
120  _bet[4] = nx*((-830251*_n-158400)*_n+197865)/7257600;
121  nx *= _n;
122  _alp[5] = (3438171-13675556*_n)*nx/7983360;
123  _bet[5] = (453717-435388*_n)*nx/15966720;
124  nx *= _n;
125  _alp[6] = 212378941*nx/319334400;
126  _bet[6] = 20648693*nx/638668800;
127  break;
128  case 7:
129  _b1 = 1/(1+_n)*(nx*(nx*(nx+4)+64)+256)/256;
130  _alp[1] = _n*(_n*(_n*(_n*(_n*(_n*(1804025*_n+2020096)-4267200)+2204160)+
131  3024000)-6451200)+4838400)/9676800;
132  _bet[1] = _n*(_n*(_n*(_n*(_n*((6156736-5406467*_n)*_n-6123600)-107520)+
133  14918400)-25804800)+19353600)/38707200;
134  _alp[2] = nx*(_n*(_n*(_n*(_n*(4626384*_n-9917165)+4316160)+3743040)-
135  5806080)+2620800)/9676800;
136  _bet[2] = nx*(_n*(_n*(_n*(_n*(829456*_n-5593555)+8478720)-5873280)+
137  1290240)+403200)/19353600;
138  nx *= _n;
139  _alp[3] = nx*(_n*(_n*((26816480-67102379*_n)*_n+16265880)-21358080)+
140  7378560)/29030400;
141  _bet[3] = nx*(_n*(_n*(_n*(9261899*_n+3564160)-2708640)-2557440)+
142  2056320)/58060800;
143  nx *= _n;
144  _alp[4] = nx*(_n*(_n*(155912000*_n+72618271)-85060800)+24532695)/
145  79833600;
146  _bet[4] = nx*(_n*(_n*(14928352*_n-9132761)-1742400)+2176515)/79833600;
147  nx *= _n;
148  _alp[5] = nx*(_n*(102508609*_n-109404448)+27505368)/63866880;
149  _bet[5] = nx*((-8005831*_n-1741552)*_n+1814868)/63866880;
150  nx *= _n;
151  _alp[6] = (2760926233LL-12282192400LL*_n)*nx/4151347200LL;
152  _bet[6] = (268433009-261810608*_n)*nx/8302694400LL;
153  nx *= _n;
154  _alp[7] = 1522256789LL*nx/1383782400LL;
155  _bet[7] = 219941297*nx/5535129600LL;
156  break;
157  case 8:
158  _b1 = 1/(1+_n)*(nx*(nx*(nx*(25*nx+64)+256)+4096)+16384)/16384;
159  _alp[1] = _n*(_n*(_n*(_n*(_n*(_n*((37884525-75900428*_n)*_n+42422016)-
160  89611200)+46287360)+63504000)-135475200)+
161  101606400)/203212800;
162  _bet[1] = _n*(_n*(_n*(_n*(_n*(_n*(_n*(31777436*_n-37845269)+43097152)-
163  42865200)-752640)+104428800)-180633600)+
164  135475200)/270950400;
165  _alp[2] = nx*(_n*(_n*(_n*(_n*(_n*(148003883*_n+83274912)-178508970)+
166  77690880)+67374720)-104509440)+47174400)/
167  174182400;
168  _bet[2] = nx*(_n*(_n*(_n*(_n*(_n*(24749483*_n+14930208)-100683990)+
169  152616960)-105719040)+23224320)+7257600)/
170  348364800;
171  nx *= _n;
172  _alp[3] = nx*(_n*(_n*(_n*(_n*(318729724*_n-738126169)+294981280)+
173  178924680)-234938880)+81164160)/319334400;
174  _bet[3] = nx*(_n*(_n*(_n*((101880889-232468668*_n)*_n+39205760)-
175  29795040)-28131840)+22619520)/638668800;
176  nx *= _n;
177  _alp[4] = nx*(_n*(_n*((14967552000LL-40176129013LL*_n)*_n+6971354016LL)-
178  8165836800LL)+2355138720LL)/7664025600LL;
179  _bet[4] = nx*(_n*(_n*(_n*(324154477*_n+1433121792LL)-876745056)-
180  167270400)+208945440)/7664025600LL;
181  nx *= _n;
182  _alp[5] = nx*(_n*(_n*(10421654396LL*_n+3997835751LL)-4266773472LL)+
183  1072709352LL)/2490808320LL;
184  _bet[5] = nx*(_n*(_n*(457888660*_n-312227409)-67920528)+70779852)/
185  2490808320LL;
186  nx *= _n;
187  _alp[6] = nx*(_n*(175214326799LL*_n-171950693600LL)+38652967262LL)/
188  58118860800LL;
189  _bet[6] = nx*((-19841813847LL*_n-3665348512LL)*_n+3758062126LL)/
190  116237721600LL;
191  nx *= _n;
192  _alp[7] = (13700311101LL-67039739596LL*_n)*nx/12454041600LL;
193  _bet[7] = (1979471673LL-1989295244LL*_n)*nx/49816166400LL;
194  nx *= _n;
195  _alp[8] = 1424729850961LL*nx/743921418240LL;
196  _bet[8] = 191773887257LL*nx/3719607091200LL;
197  break;
198  default:
199  GEOGRAPHICLIB_STATIC_ASSERT(maxpow_ >= 4 && maxpow_ <= 8,
200  "Bad value of maxpow_");
201  }
202  // _a1 is the equivalent radius for computing the circumference of
203  // ellipse.
204  _a1 = _b1 * _a;
205  }
206 
208  static const TransverseMercator utm(Constants::WGS84_a(),
211  return utm;
212  }
213 
214  // Engsager and Poder (2007) use trigonometric series to convert between phi
215  // and phip. Here are the series...
216  //
217  // Conversion from phi to phip:
218  //
219  // phip = phi + sum(c[j] * sin(2*j*phi), j, 1, 6)
220  //
221  // c[1] = - 2 * n
222  // + 2/3 * n^2
223  // + 4/3 * n^3
224  // - 82/45 * n^4
225  // + 32/45 * n^5
226  // + 4642/4725 * n^6;
227  // c[2] = 5/3 * n^2
228  // - 16/15 * n^3
229  // - 13/9 * n^4
230  // + 904/315 * n^5
231  // - 1522/945 * n^6;
232  // c[3] = - 26/15 * n^3
233  // + 34/21 * n^4
234  // + 8/5 * n^5
235  // - 12686/2835 * n^6;
236  // c[4] = 1237/630 * n^4
237  // - 12/5 * n^5
238  // - 24832/14175 * n^6;
239  // c[5] = - 734/315 * n^5
240  // + 109598/31185 * n^6;
241  // c[6] = 444337/155925 * n^6;
242  //
243  // Conversion from phip to phi:
244  //
245  // phi = phip + sum(d[j] * sin(2*j*phip), j, 1, 6)
246  //
247  // d[1] = 2 * n
248  // - 2/3 * n^2
249  // - 2 * n^3
250  // + 116/45 * n^4
251  // + 26/45 * n^5
252  // - 2854/675 * n^6;
253  // d[2] = 7/3 * n^2
254  // - 8/5 * n^3
255  // - 227/45 * n^4
256  // + 2704/315 * n^5
257  // + 2323/945 * n^6;
258  // d[3] = 56/15 * n^3
259  // - 136/35 * n^4
260  // - 1262/105 * n^5
261  // + 73814/2835 * n^6;
262  // d[4] = 4279/630 * n^4
263  // - 332/35 * n^5
264  // - 399572/14175 * n^6;
265  // d[5] = 4174/315 * n^5
266  // - 144838/6237 * n^6;
267  // d[6] = 601676/22275 * n^6;
268  //
269  // In order to maintain sufficient relative accuracy close to the pole use
270  //
271  // S = sum(c[i]*sin(2*i*phi),i,1,6)
272  // taup = (tau + tan(S)) / (1 - tau * tan(S))
273 
274  // Here we evaluate the forward transform explicitly and solve the reverse
275  // one by Newton's method.
276  //
277  // taupf and tauf are adapted from TransverseMercatorExact (taup and
278  // taupinv). tau = tan(phi), taup = sinh(psi)
279  Math::real TransverseMercator::taupf(real tau) const {
280  if (!(abs(tau) < overflow()))
281  return tau;
282  real
283  tau1 = Math::hypot(real(1), tau),
284  sig = sinh( eatanhe(tau / tau1) );
285  return Math::hypot(real(1), sig) * tau - sig * tau1;
286  }
287 
288  Math::real TransverseMercator::tauf(real taup) const {
289  if (!(abs(taup) < overflow()))
290  return taup;
291  real
292  // To lowest order in e^2, taup = (1 - e^2) * tau = _e2m * tau; so use
293  // tau = taup/_e2m as a starting guess. (This starting guess is the
294  // geocentric latitude which, to first order in the flattening, is equal
295  // to the conformal latitude.) Only 1 iteration is needed for |lat| <
296  // 3.35 deg, otherwise 2 iterations are needed. If, instead, tau = taup
297  // is used the mean number of iterations increases to 1.99 (2 iterations
298  // are needed except near tau = 0).
299  tau = taup/_e2m,
300  stol = tol_ * max(real(1), abs(taup));
301  // min iterations = 1, max iterations = 2; mean = 1.94
302  for (int i = 0; i < numit_ || GEOGRAPHICLIB_PANIC; ++i) {
303  real
304  tau1 = Math::hypot(real(1), tau),
305  sig = sinh( eatanhe( tau / tau1 ) ),
306  taupa = Math::hypot(real(1), sig) * tau - sig * tau1,
307  dtau = (taup - taupa) * (1 + _e2m * Math::sq(tau)) /
308  ( _e2m * tau1 * Math::hypot(real(1), taupa) );
309  tau += dtau;
310  if (!(abs(dtau) >= stol))
311  break;
312  }
313  return tau;
314  }
315 
316  void TransverseMercator::Forward(real lon0, real lat, real lon,
317  real& x, real& y, real& gamma, real& k)
318  const {
320  // Explicitly enforce the parity
321  int
322  latsign = lat < 0 ? -1 : 1,
323  lonsign = lon < 0 ? -1 : 1;
324  lon *= lonsign;
325  lat *= latsign;
326  bool backside = lon > 90;
327  if (backside) {
328  if (lat == 0)
329  latsign = -1;
330  lon = 180 - lon;
331  }
332  real
333  phi = lat * Math::degree(),
334  lam = lon * Math::degree();
335  // phi = latitude
336  // phi' = conformal latitude
337  // psi = isometric latitude
338  // tau = tan(phi)
339  // tau' = tan(phi')
340  // [xi', eta'] = Gauss-Schreiber TM coordinates
341  // [xi, eta] = Gauss-Krueger TM coordinates
342  //
343  // We use
344  // tan(phi') = sinh(psi)
345  // sin(phi') = tanh(psi)
346  // cos(phi') = sech(psi)
347  // denom^2 = 1-cos(phi')^2*sin(lam)^2 = 1-sech(psi)^2*sin(lam)^2
348  // sin(xip) = sin(phi')/denom = tanh(psi)/denom
349  // cos(xip) = cos(phi')*cos(lam)/denom = sech(psi)*cos(lam)/denom
350  // cosh(etap) = 1/denom = 1/denom
351  // sinh(etap) = cos(phi')*sin(lam)/denom = sech(psi)*sin(lam)/denom
352  real etap, xip;
353  if (lat != 90) {
354  real
355  c = max(real(0), cos(lam)), // cos(pi/2) might be negative
356  tau = tan(phi),
357  taup = taupf(tau);
358  xip = atan2(taup, c);
359  // Used to be
360  // etap = Math::atanh(sin(lam) / cosh(psi));
361  etap = Math::asinh(sin(lam) / Math::hypot(taup, c));
362  // convergence and scale for Gauss-Schreiber TM (xip, etap) -- gamma0 =
363  // atan(tan(xip) * tanh(etap)) = atan(tan(lam) * sin(phi'));
364  // sin(phi') = tau'/sqrt(1 + tau'^2)
365  gamma = atan(tanx(lam) *
366  taup / Math::hypot(real(1), taup)); // Krueger p 22 (44)
367  // k0 = sqrt(1 - _e2 * sin(phi)^2) * (cos(phi') / cos(phi)) * cosh(etap)
368  // Note 1/cos(phi) = cosh(psip);
369  // and cos(phi') * cosh(etap) = 1/hypot(sinh(psi), cos(lam))
370  //
371  // This form has cancelling errors. This property is lost if cosh(psip)
372  // is replaced by 1/cos(phi), even though it's using "primary" data (phi
373  // instead of psip).
374  k = sqrt(_e2m + _e2 * Math::sq(cos(phi))) * Math::hypot(real(1), tau)
375  / Math::hypot(taup, c);
376  } else {
377  xip = Math::pi()/2;
378  etap = 0;
379  gamma = lam;
380  k = _c;
381  }
382  // {xi',eta'} is {northing,easting} for Gauss-Schreiber transverse Mercator
383  // (for eta' = 0, xi' = bet). {xi,eta} is {northing,easting} for transverse
384  // Mercator with constant scale on the central meridian (for eta = 0, xip =
385  // rectifying latitude). Define
386  //
387  // zeta = xi + i*eta
388  // zeta' = xi' + i*eta'
389  //
390  // The conversion from conformal to rectifying latitude can be expressed as
391  // a series in _n:
392  //
393  // zeta = zeta' + sum(h[j-1]' * sin(2 * j * zeta'), j = 1..maxpow_)
394  //
395  // where h[j]' = O(_n^j). The reversion of this series gives
396  //
397  // zeta' = zeta - sum(h[j-1] * sin(2 * j * zeta), j = 1..maxpow_)
398  //
399  // which is used in Reverse.
400  //
401  // Evaluate sums via Clenshaw method. See
402  // http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
403  //
404  // Let
405  //
406  // S = sum(c[k] * F[k](x), k = 0..N)
407  // F[n+1](x) = alpha(n,x) * F[n](x) + beta(n,x) * F[n-1](x)
408  //
409  // Evaluate S with
410  //
411  // y[N+2] = y[N+1] = 0
412  // y[k] = alpha(k,x) * y[k+1] + beta(k+1,x) * y[k+2] + c[k]
413  // S = c[0] * F[0](x) + y[1] * F[1](x) + beta(1,x) * F[0](x) * y[2]
414  //
415  // Here we have
416  //
417  // x = 2 * zeta'
418  // F[n](x) = sin(n * x)
419  // a(n, x) = 2 * cos(x)
420  // b(n, x) = -1
421  // [ sin(A+B) - 2*cos(B)*sin(A) + sin(A-B) = 0, A = n*x, B = x ]
422  // N = maxpow_
423  // c[k] = _alp[k]
424  // S = y[1] * sin(x)
425  //
426  // For the derivative we have
427  //
428  // x = 2 * zeta'
429  // F[n](x) = cos(n * x)
430  // a(n, x) = 2 * cos(x)
431  // b(n, x) = -1
432  // [ cos(A+B) - 2*cos(B)*cos(A) + cos(A-B) = 0, A = n*x, B = x ]
433  // c[0] = 1; c[k] = 2*k*_alp[k]
434  // S = (c[0] - y[2]) + y[1] * cos(x)
435  real
436  c0 = cos(2 * xip), ch0 = cosh(2 * etap),
437  s0 = sin(2 * xip), sh0 = sinh(2 * etap),
438  ar = 2 * c0 * ch0, ai = -2 * s0 * sh0; // 2 * cos(2*zeta')
439  int n = maxpow_;
440  real
441  xi0 = (n & 1 ? _alp[n] : 0), eta0 = 0,
442  xi1 = 0, eta1 = 0;
443  real // Accumulators for dzeta/dzeta'
444  yr0 = (n & 1 ? 2 * maxpow_ * _alp[n--] : 0), yi0 = 0,
445  yr1 = 0, yi1 = 0;
446  while (n) {
447  xi1 = ar * xi0 - ai * eta0 - xi1 + _alp[n];
448  eta1 = ai * xi0 + ar * eta0 - eta1;
449  yr1 = ar * yr0 - ai * yi0 - yr1 + 2 * n * _alp[n];
450  yi1 = ai * yr0 + ar * yi0 - yi1;
451  --n;
452  xi0 = ar * xi1 - ai * eta1 - xi0 + _alp[n];
453  eta0 = ai * xi1 + ar * eta1 - eta0;
454  yr0 = ar * yr1 - ai * yi1 - yr0 + 2 * n * _alp[n];
455  yi0 = ai * yr1 + ar * yi1 - yi0;
456  --n;
457  }
458  ar /= 2; ai /= 2; // cos(2*zeta')
459  yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
460  yi1 = - yi1 + ai * yr0 + ar * yi0;
461  ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta')
462  real
463  xi = xip + ar * xi0 - ai * eta0,
464  eta = etap + ai * xi0 + ar * eta0;
465  // Fold in change in convergence and scale for Gauss-Schreiber TM to
466  // Gauss-Krueger TM.
467  gamma -= atan2(yi1, yr1);
468  k *= _b1 * Math::hypot(yr1, yi1);
469  gamma /= Math::degree();
470  y = _a1 * _k0 * (backside ? Math::pi() - xi : xi) * latsign;
471  x = _a1 * _k0 * eta * lonsign;
472  if (backside)
473  gamma = 180 - gamma;
474  gamma *= latsign * lonsign;
475  k *= _k0;
476  }
477 
478  void TransverseMercator::Reverse(real lon0, real x, real y,
479  real& lat, real& lon, real& gamma, real& k)
480  const {
481  // This undoes the steps in Forward. The wrinkles are: (1) Use of the
482  // reverted series to express zeta' in terms of zeta. (2) Newton's method
483  // to solve for phi in terms of tan(phi).
484  real
485  xi = y / (_a1 * _k0),
486  eta = x / (_a1 * _k0);
487  // Explicitly enforce the parity
488  int
489  xisign = xi < 0 ? -1 : 1,
490  etasign = eta < 0 ? -1 : 1;
491  xi *= xisign;
492  eta *= etasign;
493  bool backside = xi > Math::pi()/2;
494  if (backside)
495  xi = Math::pi() - xi;
496  real
497  c0 = cos(2 * xi), ch0 = cosh(2 * eta),
498  s0 = sin(2 * xi), sh0 = sinh(2 * eta),
499  ar = 2 * c0 * ch0, ai = -2 * s0 * sh0; // 2 * cos(2*zeta)
500  int n = maxpow_;
501  real // Accumulators for zeta'
502  xip0 = (n & 1 ? -_bet[n] : 0), etap0 = 0,
503  xip1 = 0, etap1 = 0;
504  real // Accumulators for dzeta'/dzeta
505  yr0 = (n & 1 ? - 2 * maxpow_ * _bet[n--] : 0), yi0 = 0,
506  yr1 = 0, yi1 = 0;
507  while (n) {
508  xip1 = ar * xip0 - ai * etap0 - xip1 - _bet[n];
509  etap1 = ai * xip0 + ar * etap0 - etap1;
510  yr1 = ar * yr0 - ai * yi0 - yr1 - 2 * n * _bet[n];
511  yi1 = ai * yr0 + ar * yi0 - yi1;
512  --n;
513  xip0 = ar * xip1 - ai * etap1 - xip0 - _bet[n];
514  etap0 = ai * xip1 + ar * etap1 - etap0;
515  yr0 = ar * yr1 - ai * yi1 - yr0 - 2 * n * _bet[n];
516  yi0 = ai * yr1 + ar * yi1 - yi0;
517  --n;
518  }
519  ar /= 2; ai /= 2; // cos(2*zeta')
520  yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
521  yi1 = - yi1 + ai * yr0 + ar * yi0;
522  ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta)
523  real
524  xip = xi + ar * xip0 - ai * etap0,
525  etap = eta + ai * xip0 + ar * etap0;
526  // Convergence and scale for Gauss-Schreiber TM to Gauss-Krueger TM.
527  gamma = atan2(yi1, yr1);
528  k = _b1 / Math::hypot(yr1, yi1);
529  // JHS 154 has
530  //
531  // phi' = asin(sin(xi') / cosh(eta')) (Krueger p 17 (25))
532  // lam = asin(tanh(eta') / cos(phi')
533  // psi = asinh(tan(phi'))
534  real lam, phi;
535  real
536  s = sinh(etap),
537  c = max(real(0), cos(xip)), // cos(pi/2) might be negative
538  r = Math::hypot(s, c);
539  if (r != 0) {
540  lam = atan2(s, c); // Krueger p 17 (25)
541  // Use Newton's method to solve for tau
542  real
543  taup = sin(xip)/r,
544  tau = tauf(taup);
545  phi = atan(tau);
546  gamma += atan(tanx(xip) * tanh(etap)); // Krueger p 19 (31)
547  // Note cos(phi') * cosh(eta') = r
548  k *= sqrt(_e2m + _e2 * Math::sq(cos(phi))) *
549  Math::hypot(real(1), tau) * r;
550  } else {
551  phi = Math::pi()/2;
552  lam = 0;
553  k *= _c;
554  }
555  lat = phi / Math::degree() * xisign;
556  lon = lam / Math::degree();
557  if (backside)
558  lon = 180 - lon;
559  lon *= etasign;
560  lon = Math::AngNormalize(lon + Math::AngNormalize(lon0));
561  gamma /= Math::degree();
562  if (backside)
563  gamma = 180 - gamma;
564  gamma *= xisign * etasign;
565  k *= _k0;
566  }
567 
568 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:400
static T pi()
Definition: Math.hpp:214
GeographicLib::Math::real real
Definition: GeodSolve.cpp:32
static bool isfinite(T x)
Definition: Math.hpp:446
void Forward(real lon0, real lat, real lon, real &x, real &y, real &gamma, real &k) const
Transverse Mercator projection.
static const TransverseMercator & UTM()
Header for GeographicLib::TransverseMercator class.
static T asinh(T x)
Definition: Math.hpp:323
TransverseMercator(real a, real f, real k0)
static T hypot(T x, T y)
Definition: Math.hpp:255
static T sq(T x)
Definition: Math.hpp:244
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T degree()
Definition: Math.hpp:228
static T AngDiff(T x, T y)
Definition: Math.hpp:430
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
Exception handling for GeographicLib.
Definition: Constants.hpp:361
#define GEOGRAPHICLIB_PANIC
Definition: Math.hpp:87