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