GeographicLib  1.45
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-2015) <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 namespace GeographicLib {
45 
46  using namespace std;
47 
48  TransverseMercator::TransverseMercator(real a, real f, real k0)
49  : _a(a)
50  , _f(f <= 1 ? f : 1/f) // f > 1 behavior is DEPRECATED
51  , _k0(k0)
52  , _e2(_f * (2 - _f))
53  , _es((f < 0 ? -1 : 1) * sqrt(abs(_e2)))
54  , _e2m(1 - _e2)
55  // _c = sqrt( pow(1 + _e, 1 + _e) * pow(1 - _e, 1 - _e) ) )
56  // See, for example, Lee (1976), p 100.
57  , _c( sqrt(_e2m) * exp(Math::eatanhe(real(1), _es)) )
58  , _n(_f / (2 - _f))
59  {
60  if (!(Math::isfinite(_a) && _a > 0))
61  throw GeographicErr("Major radius is not positive");
62  if (!(Math::isfinite(_f) && _f < 1))
63  throw GeographicErr("Minor radius is not positive");
64  if (!(Math::isfinite(_k0) && _k0 > 0))
65  throw GeographicErr("Scale is not positive");
66 
67  // Generated by Maxima on 2015-05-14 22:55:13-04:00
68 #if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 2
69  static const real b1coeff[] = {
70  // b1*(n+1), polynomial in n2 of order 2
71  1, 16, 64, 64,
72  }; // count = 4
73 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 3
74  static const real b1coeff[] = {
75  // b1*(n+1), polynomial in n2 of order 3
76  1, 4, 64, 256, 256,
77  }; // count = 5
78 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER/2 == 4
79  static const real b1coeff[] = {
80  // b1*(n+1), polynomial in n2 of order 4
81  25, 64, 256, 4096, 16384, 16384,
82  }; // count = 6
83 #else
84 #error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER"
85 #endif
86 
87 #if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4
88  static const real alpcoeff[] = {
89  // alp[1]/n^1, polynomial in n of order 3
90  164, 225, -480, 360, 720,
91  // alp[2]/n^2, polynomial in n of order 2
92  557, -864, 390, 1440,
93  // alp[3]/n^3, polynomial in n of order 1
94  -1236, 427, 1680,
95  // alp[4]/n^4, polynomial in n of order 0
96  49561, 161280,
97  }; // count = 14
98 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5
99  static const real alpcoeff[] = {
100  // alp[1]/n^1, polynomial in n of order 4
101  -635, 328, 450, -960, 720, 1440,
102  // alp[2]/n^2, polynomial in n of order 3
103  4496, 3899, -6048, 2730, 10080,
104  // alp[3]/n^3, polynomial in n of order 2
105  15061, -19776, 6832, 26880,
106  // alp[4]/n^4, polynomial in n of order 1
107  -171840, 49561, 161280,
108  // alp[5]/n^5, polynomial in n of order 0
109  34729, 80640,
110  }; // count = 20
111 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6
112  static const real alpcoeff[] = {
113  // alp[1]/n^1, polynomial in n of order 5
114  31564, -66675, 34440, 47250, -100800, 75600, 151200,
115  // alp[2]/n^2, polynomial in n of order 4
116  -1983433, 863232, 748608, -1161216, 524160, 1935360,
117  // alp[3]/n^3, polynomial in n of order 3
118  670412, 406647, -533952, 184464, 725760,
119  // alp[4]/n^4, polynomial in n of order 2
120  6601661, -7732800, 2230245, 7257600,
121  // alp[5]/n^5, polynomial in n of order 1
122  -13675556, 3438171, 7983360,
123  // alp[6]/n^6, polynomial in n of order 0
124  212378941, 319334400,
125  }; // count = 27
126 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7
127  static const real alpcoeff[] = {
128  // alp[1]/n^1, polynomial in n of order 6
129  1804025, 2020096, -4267200, 2204160, 3024000, -6451200, 4838400, 9676800,
130  // alp[2]/n^2, polynomial in n of order 5
131  4626384, -9917165, 4316160, 3743040, -5806080, 2620800, 9676800,
132  // alp[3]/n^3, polynomial in n of order 4
133  -67102379, 26816480, 16265880, -21358080, 7378560, 29030400,
134  // alp[4]/n^4, polynomial in n of order 3
135  155912000, 72618271, -85060800, 24532695, 79833600,
136  // alp[5]/n^5, polynomial in n of order 2
137  102508609, -109404448, 27505368, 63866880,
138  // alp[6]/n^6, polynomial in n of order 1
139  -12282192400LL, 2760926233LL, 4151347200LL,
140  // alp[7]/n^7, polynomial in n of order 0
141  1522256789, 1383782400,
142  }; // count = 35
143 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8
144  static const real alpcoeff[] = {
145  // alp[1]/n^1, polynomial in n of order 7
146  -75900428, 37884525, 42422016, -89611200, 46287360, 63504000, -135475200,
147  101606400, 203212800,
148  // alp[2]/n^2, polynomial in n of order 6
149  148003883, 83274912, -178508970, 77690880, 67374720, -104509440,
150  47174400, 174182400,
151  // alp[3]/n^3, polynomial in n of order 5
152  318729724, -738126169, 294981280, 178924680, -234938880, 81164160,
153  319334400,
154  // alp[4]/n^4, polynomial in n of order 4
155  -40176129013LL, 14967552000LL, 6971354016LL, -8165836800LL, 2355138720LL,
156  7664025600LL,
157  // alp[5]/n^5, polynomial in n of order 3
158  10421654396LL, 3997835751LL, -4266773472LL, 1072709352, 2490808320LL,
159  // alp[6]/n^6, polynomial in n of order 2
160  175214326799LL, -171950693600LL, 38652967262LL, 58118860800LL,
161  // alp[7]/n^7, polynomial in n of order 1
162  -67039739596LL, 13700311101LL, 12454041600LL,
163  // alp[8]/n^8, polynomial in n of order 0
164  1424729850961LL, 743921418240LL,
165  }; // count = 44
166 #else
167 #error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER"
168 #endif
169 
170 #if GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 4
171  static const real betcoeff[] = {
172  // bet[1]/n^1, polynomial in n of order 3
173  -4, 555, -960, 720, 1440,
174  // bet[2]/n^2, polynomial in n of order 2
175  -437, 96, 30, 1440,
176  // bet[3]/n^3, polynomial in n of order 1
177  -148, 119, 3360,
178  // bet[4]/n^4, polynomial in n of order 0
179  4397, 161280,
180  }; // count = 14
181 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 5
182  static const real betcoeff[] = {
183  // bet[1]/n^1, polynomial in n of order 4
184  -3645, -64, 8880, -15360, 11520, 23040,
185  // bet[2]/n^2, polynomial in n of order 3
186  4416, -3059, 672, 210, 10080,
187  // bet[3]/n^3, polynomial in n of order 2
188  -627, -592, 476, 13440,
189  // bet[4]/n^4, polynomial in n of order 1
190  -3520, 4397, 161280,
191  // bet[5]/n^5, polynomial in n of order 0
192  4583, 161280,
193  }; // count = 20
194 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 6
195  static const real betcoeff[] = {
196  // bet[1]/n^1, polynomial in n of order 5
197  384796, -382725, -6720, 932400, -1612800, 1209600, 2419200,
198  // bet[2]/n^2, polynomial in n of order 4
199  -1118711, 1695744, -1174656, 258048, 80640, 3870720,
200  // bet[3]/n^3, polynomial in n of order 3
201  22276, -16929, -15984, 12852, 362880,
202  // bet[4]/n^4, polynomial in n of order 2
203  -830251, -158400, 197865, 7257600,
204  // bet[5]/n^5, polynomial in n of order 1
205  -435388, 453717, 15966720,
206  // bet[6]/n^6, polynomial in n of order 0
207  20648693, 638668800,
208  }; // count = 27
209 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 7
210  static const real betcoeff[] = {
211  // bet[1]/n^1, polynomial in n of order 6
212  -5406467, 6156736, -6123600, -107520, 14918400, -25804800, 19353600,
213  38707200,
214  // bet[2]/n^2, polynomial in n of order 5
215  829456, -5593555, 8478720, -5873280, 1290240, 403200, 19353600,
216  // bet[3]/n^3, polynomial in n of order 4
217  9261899, 3564160, -2708640, -2557440, 2056320, 58060800,
218  // bet[4]/n^4, polynomial in n of order 3
219  14928352, -9132761, -1742400, 2176515, 79833600,
220  // bet[5]/n^5, polynomial in n of order 2
221  -8005831, -1741552, 1814868, 63866880,
222  // bet[6]/n^6, polynomial in n of order 1
223  -261810608, 268433009, 8302694400LL,
224  // bet[7]/n^7, polynomial in n of order 0
225  219941297, 5535129600LL,
226  }; // count = 35
227 #elif GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER == 8
228  static const real betcoeff[] = {
229  // bet[1]/n^1, polynomial in n of order 7
230  31777436, -37845269, 43097152, -42865200, -752640, 104428800, -180633600,
231  135475200, 270950400,
232  // bet[2]/n^2, polynomial in n of order 6
233  24749483, 14930208, -100683990, 152616960, -105719040, 23224320, 7257600,
234  348364800,
235  // bet[3]/n^3, polynomial in n of order 5
236  -232468668, 101880889, 39205760, -29795040, -28131840, 22619520,
237  638668800,
238  // bet[4]/n^4, polynomial in n of order 4
239  324154477, 1433121792, -876745056, -167270400, 208945440, 7664025600LL,
240  // bet[5]/n^5, polynomial in n of order 3
241  457888660, -312227409, -67920528, 70779852, 2490808320LL,
242  // bet[6]/n^6, polynomial in n of order 2
243  -19841813847LL, -3665348512LL, 3758062126LL, 116237721600LL,
244  // bet[7]/n^7, polynomial in n of order 1
245  -1989295244, 1979471673, 49816166400LL,
246  // bet[8]/n^8, polynomial in n of order 0
247  191773887257LL, 3719607091200LL,
248  }; // count = 44
249 #else
250 #error "Bad value for GEOGRAPHICLIB_TRANSVERSEMERCATOR_ORDER"
251 #endif
252 
253  GEOGRAPHICLIB_STATIC_ASSERT(sizeof(b1coeff) / sizeof(real) ==
254  maxpow_/2 + 2,
255  "Coefficient array size mismatch for b1");
256  GEOGRAPHICLIB_STATIC_ASSERT(sizeof(alpcoeff) / sizeof(real) ==
257  (maxpow_ * (maxpow_ + 3))/2,
258  "Coefficient array size mismatch for alp");
259  GEOGRAPHICLIB_STATIC_ASSERT(sizeof(betcoeff) / sizeof(real) ==
260  (maxpow_ * (maxpow_ + 3))/2,
261  "Coefficient array size mismatch for bet");
262  int m = maxpow_/2;
263  _b1 = Math::polyval(m, b1coeff, Math::sq(_n)) / (b1coeff[m + 1] * (1+_n));
264  // _a1 is the equivalent radius for computing the circumference of
265  // ellipse.
266  _a1 = _b1 * _a;
267  int o = 0;
268  real d = _n;
269  for (int l = 1; l <= maxpow_; ++l) {
270  m = maxpow_ - l;
271  _alp[l] = d * Math::polyval(m, alpcoeff + o, _n) / alpcoeff[o + m + 1];
272  _bet[l] = d * Math::polyval(m, betcoeff + o, _n) / betcoeff[o + m + 1];
273  o += m + 2;
274  d *= _n;
275  }
276  // Post condition: o == sizeof(alpcoeff) / sizeof(real) &&
277  // o == sizeof(betcoeff) / sizeof(real)
278  }
279 
281  static const TransverseMercator utm(Constants::WGS84_a(),
284  return utm;
285  }
286 
287  // Engsager and Poder (2007) use trigonometric series to convert between phi
288  // and phip. Here are the series...
289  //
290  // Conversion from phi to phip:
291  //
292  // phip = phi + sum(c[j] * sin(2*j*phi), j, 1, 6)
293  //
294  // c[1] = - 2 * n
295  // + 2/3 * n^2
296  // + 4/3 * n^3
297  // - 82/45 * n^4
298  // + 32/45 * n^5
299  // + 4642/4725 * n^6;
300  // c[2] = 5/3 * n^2
301  // - 16/15 * n^3
302  // - 13/9 * n^4
303  // + 904/315 * n^5
304  // - 1522/945 * n^6;
305  // c[3] = - 26/15 * n^3
306  // + 34/21 * n^4
307  // + 8/5 * n^5
308  // - 12686/2835 * n^6;
309  // c[4] = 1237/630 * n^4
310  // - 12/5 * n^5
311  // - 24832/14175 * n^6;
312  // c[5] = - 734/315 * n^5
313  // + 109598/31185 * n^6;
314  // c[6] = 444337/155925 * n^6;
315  //
316  // Conversion from phip to phi:
317  //
318  // phi = phip + sum(d[j] * sin(2*j*phip), j, 1, 6)
319  //
320  // d[1] = 2 * n
321  // - 2/3 * n^2
322  // - 2 * n^3
323  // + 116/45 * n^4
324  // + 26/45 * n^5
325  // - 2854/675 * n^6;
326  // d[2] = 7/3 * n^2
327  // - 8/5 * n^3
328  // - 227/45 * n^4
329  // + 2704/315 * n^5
330  // + 2323/945 * n^6;
331  // d[3] = 56/15 * n^3
332  // - 136/35 * n^4
333  // - 1262/105 * n^5
334  // + 73814/2835 * n^6;
335  // d[4] = 4279/630 * n^4
336  // - 332/35 * n^5
337  // - 399572/14175 * n^6;
338  // d[5] = 4174/315 * n^5
339  // - 144838/6237 * n^6;
340  // d[6] = 601676/22275 * n^6;
341  //
342  // In order to maintain sufficient relative accuracy close to the pole use
343  //
344  // S = sum(c[i]*sin(2*i*phi),i,1,6)
345  // taup = (tau + tan(S)) / (1 - tau * tan(S))
346 
347  // In Math::taupf and Math::tauf we evaluate the forward transform explicitly
348  // and solve the reverse one by Newton's method.
349  //
350  // There are adapted from TransverseMercatorExact (taup and taupinv). tau =
351  // tan(phi), taup = sinh(psi)
352 
353  void TransverseMercator::Forward(real lon0, real lat, real lon,
354  real& x, real& y, real& gamma, real& k)
355  const {
356  lat = Math::LatFix(lat);
357  lon = Math::AngDiff(lon0, lon);
358  // Explicitly enforce the parity
359  int
360  latsign = lat < 0 ? -1 : 1,
361  lonsign = lon < 0 ? -1 : 1;
362  lon *= lonsign;
363  lat *= latsign;
364  bool backside = lon > 90;
365  if (backside) {
366  if (lat == 0)
367  latsign = -1;
368  lon = 180 - lon;
369  }
370  real sphi, cphi, slam, clam;
371  Math::sincosd(lat, sphi, cphi);
372  Math::sincosd(lon, slam, clam);
373  // phi = latitude
374  // phi' = conformal latitude
375  // psi = isometric latitude
376  // tau = tan(phi)
377  // tau' = tan(phi')
378  // [xi', eta'] = Gauss-Schreiber TM coordinates
379  // [xi, eta] = Gauss-Krueger TM coordinates
380  //
381  // We use
382  // tan(phi') = sinh(psi)
383  // sin(phi') = tanh(psi)
384  // cos(phi') = sech(psi)
385  // denom^2 = 1-cos(phi')^2*sin(lam)^2 = 1-sech(psi)^2*sin(lam)^2
386  // sin(xip) = sin(phi')/denom = tanh(psi)/denom
387  // cos(xip) = cos(phi')*cos(lam)/denom = sech(psi)*cos(lam)/denom
388  // cosh(etap) = 1/denom = 1/denom
389  // sinh(etap) = cos(phi')*sin(lam)/denom = sech(psi)*sin(lam)/denom
390  real etap, xip;
391  if (lat != 90) {
392  real
393  tau = sphi / cphi,
394  taup = Math::taupf(tau, _es);
395  xip = atan2(taup, clam);
396  // Used to be
397  // etap = Math::atanh(sin(lam) / cosh(psi));
398  etap = Math::asinh(slam / Math::hypot(taup, clam));
399  // convergence and scale for Gauss-Schreiber TM (xip, etap) -- gamma0 =
400  // atan(tan(xip) * tanh(etap)) = atan(tan(lam) * sin(phi'));
401  // sin(phi') = tau'/sqrt(1 + tau'^2)
402  // Krueger p 22 (44)
403  gamma = Math::atan2d(slam * taup, clam * Math::hypot(real(1), taup));
404  // k0 = sqrt(1 - _e2 * sin(phi)^2) * (cos(phi') / cos(phi)) * cosh(etap)
405  // Note 1/cos(phi) = cosh(psip);
406  // and cos(phi') * cosh(etap) = 1/hypot(sinh(psi), cos(lam))
407  //
408  // This form has cancelling errors. This property is lost if cosh(psip)
409  // is replaced by 1/cos(phi), even though it's using "primary" data (phi
410  // instead of psip).
411  k = sqrt(_e2m + _e2 * Math::sq(cphi)) * Math::hypot(real(1), tau)
412  / Math::hypot(taup, clam);
413  } else {
414  xip = Math::pi()/2;
415  etap = 0;
416  gamma = lon;
417  k = _c;
418  }
419  // {xi',eta'} is {northing,easting} for Gauss-Schreiber transverse Mercator
420  // (for eta' = 0, xi' = bet). {xi,eta} is {northing,easting} for transverse
421  // Mercator with constant scale on the central meridian (for eta = 0, xip =
422  // rectifying latitude). Define
423  //
424  // zeta = xi + i*eta
425  // zeta' = xi' + i*eta'
426  //
427  // The conversion from conformal to rectifying latitude can be expressed as
428  // a series in _n:
429  //
430  // zeta = zeta' + sum(h[j-1]' * sin(2 * j * zeta'), j = 1..maxpow_)
431  //
432  // where h[j]' = O(_n^j). The reversion of this series gives
433  //
434  // zeta' = zeta - sum(h[j-1] * sin(2 * j * zeta), j = 1..maxpow_)
435  //
436  // which is used in Reverse.
437  //
438  // Evaluate sums via Clenshaw method. See
439  // http://mathworld.wolfram.com/ClenshawRecurrenceFormula.html
440  //
441  // Let
442  //
443  // S = sum(c[k] * F[k](x), k = 0..N)
444  // F[n+1](x) = alpha(n,x) * F[n](x) + beta(n,x) * F[n-1](x)
445  //
446  // Evaluate S with
447  //
448  // y[N+2] = y[N+1] = 0
449  // y[k] = alpha(k,x) * y[k+1] + beta(k+1,x) * y[k+2] + c[k]
450  // S = c[0] * F[0](x) + y[1] * F[1](x) + beta(1,x) * F[0](x) * y[2]
451  //
452  // Here we have
453  //
454  // x = 2 * zeta'
455  // F[n](x) = sin(n * x)
456  // a(n, x) = 2 * cos(x)
457  // b(n, x) = -1
458  // [ sin(A+B) - 2*cos(B)*sin(A) + sin(A-B) = 0, A = n*x, B = x ]
459  // N = maxpow_
460  // c[k] = _alp[k]
461  // S = y[1] * sin(x)
462  //
463  // For the derivative we have
464  //
465  // x = 2 * zeta'
466  // F[n](x) = cos(n * x)
467  // a(n, x) = 2 * cos(x)
468  // b(n, x) = -1
469  // [ cos(A+B) - 2*cos(B)*cos(A) + cos(A-B) = 0, A = n*x, B = x ]
470  // c[0] = 1; c[k] = 2*k*_alp[k]
471  // S = (c[0] - y[2]) + y[1] * cos(x)
472  real
473  c0 = cos(2 * xip), ch0 = cosh(2 * etap),
474  s0 = sin(2 * xip), sh0 = sinh(2 * etap),
475  ar = 2 * c0 * ch0, ai = -2 * s0 * sh0; // 2 * cos(2*zeta')
476  int n = maxpow_;
477  real
478  xi0 = (n & 1 ? _alp[n] : 0), eta0 = 0,
479  xi1 = 0, eta1 = 0;
480  real // Accumulators for dzeta/dzeta'
481  yr0 = (n & 1 ? 2 * maxpow_ * _alp[n--] : 0), yi0 = 0,
482  yr1 = 0, yi1 = 0;
483  while (n) {
484  xi1 = ar * xi0 - ai * eta0 - xi1 + _alp[n];
485  eta1 = ai * xi0 + ar * eta0 - eta1;
486  yr1 = ar * yr0 - ai * yi0 - yr1 + 2 * n * _alp[n];
487  yi1 = ai * yr0 + ar * yi0 - yi1;
488  --n;
489  xi0 = ar * xi1 - ai * eta1 - xi0 + _alp[n];
490  eta0 = ai * xi1 + ar * eta1 - eta0;
491  yr0 = ar * yr1 - ai * yi1 - yr0 + 2 * n * _alp[n];
492  yi0 = ai * yr1 + ar * yi1 - yi0;
493  --n;
494  }
495  ar /= 2; ai /= 2; // cos(2*zeta')
496  yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
497  yi1 = - yi1 + ai * yr0 + ar * yi0;
498  ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta')
499  real
500  xi = xip + ar * xi0 - ai * eta0,
501  eta = etap + ai * xi0 + ar * eta0;
502  // Fold in change in convergence and scale for Gauss-Schreiber TM to
503  // Gauss-Krueger TM.
504  gamma -= Math::atan2d(yi1, yr1);
505  k *= _b1 * Math::hypot(yr1, yi1);
506  y = _a1 * _k0 * (backside ? Math::pi() - xi : xi) * latsign;
507  x = _a1 * _k0 * eta * lonsign;
508  if (backside)
509  gamma = 180 - gamma;
510  gamma *= latsign * lonsign;
511  gamma = Math::AngNormalize(gamma);
512  k *= _k0;
513  }
514 
515  void TransverseMercator::Reverse(real lon0, real x, real y,
516  real& lat, real& lon, real& gamma, real& k)
517  const {
518  // This undoes the steps in Forward. The wrinkles are: (1) Use of the
519  // reverted series to express zeta' in terms of zeta. (2) Newton's method
520  // to solve for phi in terms of tan(phi).
521  real
522  xi = y / (_a1 * _k0),
523  eta = x / (_a1 * _k0);
524  // Explicitly enforce the parity
525  int
526  xisign = xi < 0 ? -1 : 1,
527  etasign = eta < 0 ? -1 : 1;
528  xi *= xisign;
529  eta *= etasign;
530  bool backside = xi > Math::pi()/2;
531  if (backside)
532  xi = Math::pi() - xi;
533  real
534  c0 = cos(2 * xi), ch0 = cosh(2 * eta),
535  s0 = sin(2 * xi), sh0 = sinh(2 * eta),
536  ar = 2 * c0 * ch0, ai = -2 * s0 * sh0; // 2 * cos(2*zeta)
537  int n = maxpow_;
538  real // Accumulators for zeta'
539  xip0 = (n & 1 ? -_bet[n] : 0), etap0 = 0,
540  xip1 = 0, etap1 = 0;
541  real // Accumulators for dzeta'/dzeta
542  yr0 = (n & 1 ? - 2 * maxpow_ * _bet[n--] : 0), yi0 = 0,
543  yr1 = 0, yi1 = 0;
544  while (n) {
545  xip1 = ar * xip0 - ai * etap0 - xip1 - _bet[n];
546  etap1 = ai * xip0 + ar * etap0 - etap1;
547  yr1 = ar * yr0 - ai * yi0 - yr1 - 2 * n * _bet[n];
548  yi1 = ai * yr0 + ar * yi0 - yi1;
549  --n;
550  xip0 = ar * xip1 - ai * etap1 - xip0 - _bet[n];
551  etap0 = ai * xip1 + ar * etap1 - etap0;
552  yr0 = ar * yr1 - ai * yi1 - yr0 - 2 * n * _bet[n];
553  yi0 = ai * yr1 + ar * yi1 - yi0;
554  --n;
555  }
556  ar /= 2; ai /= 2; // cos(2*zeta')
557  yr1 = 1 - yr1 + ar * yr0 - ai * yi0;
558  yi1 = - yi1 + ai * yr0 + ar * yi0;
559  ar = s0 * ch0; ai = c0 * sh0; // sin(2*zeta)
560  real
561  xip = xi + ar * xip0 - ai * etap0,
562  etap = eta + ai * xip0 + ar * etap0;
563  // Convergence and scale for Gauss-Schreiber TM to Gauss-Krueger TM.
564  gamma = Math::atan2d(yi1, yr1);
565  k = _b1 / Math::hypot(yr1, yi1);
566  // JHS 154 has
567  //
568  // phi' = asin(sin(xi') / cosh(eta')) (Krueger p 17 (25))
569  // lam = asin(tanh(eta') / cos(phi')
570  // psi = asinh(tan(phi'))
571  real
572  s = sinh(etap),
573  c = max(real(0), cos(xip)), // cos(pi/2) might be negative
574  r = Math::hypot(s, c);
575  if (r != 0) {
576  lon = Math::atan2d(s, c); // Krueger p 17 (25)
577  // Use Newton's method to solve for tau
578  real
579  sxip = sin(xip),
580  tau = Math::tauf(sxip/r, _es);
581  gamma += Math::atan2d(sxip * tanh(etap), c); // Krueger p 19 (31)
582  lat = Math::atand(tau);
583  // Note cos(phi') * cosh(eta') = r
584  k *= sqrt(_e2m + _e2 / (1 + Math::sq(tau))) *
585  Math::hypot(real(1), tau) * r;
586  } else {
587  lat = 90;
588  lon = 0;
589  k *= _c;
590  }
591  lat *= xisign;
592  if (backside)
593  lon = 180 - lon;
594  lon *= etasign;
595  lon = Math::AngNormalize(lon + lon0);
596  if (backside)
597  gamma = 180 - gamma;
598  gamma *= xisign * etasign;
599  gamma = Math::AngNormalize(gamma);
600  k *= _k0;
601  }
602 
603 } // namespace GeographicLib
static T AngNormalize(T x)
Definition: Math.hpp:451
static T pi()
Definition: Math.hpp:216
static T atand(T x)
Definition: Math.hpp:708
static bool isfinite(T x)
Definition: Math.hpp:768
static T LatFix(T x)
Definition: Math.hpp:482
Mathematical functions needed by GeographicLib.
Definition: Math.hpp:102
static void sincosd(T x, T &sinx, T &cosx)
Definition: Math.hpp:559
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:325
TransverseMercator(real a, real f, real k0)
static T hypot(T x, T y)
Definition: Math.hpp:257
static T sq(T x)
Definition: Math.hpp:246
static T atan2d(T y, T x)
Definition: Math.hpp:676
static T polyval(int N, const T p[], T x)
Definition: Math.hpp:439
Namespace for GeographicLib.
Definition: Accumulator.cpp:12
static T AngDiff(T x, T y)
Definition: Math.hpp:499
void Reverse(real lon0, real x, real y, real &lat, real &lon, real &gamma, real &k) const
static T tauf(T taup, T es)
Exception handling for GeographicLib.
Definition: Constants.hpp:386
static T taupf(T tau, T es)