UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
TransverseMercator.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* FILE: TransverseMercator.cpp
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude) and Transverse Mercator projection coordinates
10  * (easting and northing).
11  *
12  * MODIFICATIONS
13  *
14  * Date Description
15  * ---- -----------
16  * 2-26-07 Original C++ Code
17  * 7/19/10 N. Lundgren BAEts27271 Correct the test for TRANMERC_LON_WARNING
18  * in convertToGeodetic, by removing the multiply by cos(latitude)
19  * 3/23/11 N. Lundgren BAEts28583 Updated memory leak checks for
20  * code consistency.
21  * 7/02/14 Updated to algorithm in NGA Transverse Mercator document.
22  *
23  * 1/19/16 A. Layne Updated algorithm for R4 MSP_DR30344
24  *
25  * 1/19/16 A. Layne MSP_DR30125 Updated generateCoefficients to use book values for
26  * coefficients based on supplied ellipsoid code. If user defined, use old default computation.
27  *
28  * 10/12/17 S. Holt LSC-16618 Initialized last two elements of aCoeff[] and bCoeff[] to 0.0
29  * 10/12/17 S. Holt LSC-9451 Swapped values for TranMerc_Delta_Northing and TranMerc_Delta_Easting
30  */
31 #include <iostream>
32 
33 #include <cmath>
34 #include <math.h>
35 #include "TransverseMercator.h"
36 #include "MapProjection5Parameters.h"
37 #include "MapProjectionCoordinates.h"
38 #include "GeodeticCoordinates.h"
39 #include "CoordinateConversionException.h"
40 #include "ErrorMessages.h"
41 
42 using MSP::CCS::TransverseMercator;
43 // Terms in series A and B coefficients
44 #define N_TERMS 6
45 #define MAX_TERMS 8
46 
47 // DEFINES
48 #define PI 3.14159265358979323e0
49 #define PI_OVER_2 (PI/2.0e0)
50 #define MAX_DELTA_LONG ((PI * 70)/180.0)
51 #define MIN_SCALE_FACTOR 0.1
52 #define MAX_SCALE_FACTOR 10.0
53 
54 
55 TransverseMercator::TransverseMercator(
56  double ellipsoidSemiMajorAxis,
57  double ellipsoidFlattening,
58  double centralMeridian,
59  double latitudeOfTrueScale,
60  double falseEasting,
61  double falseNorthing,
62  double scaleFactor,
63  char *ellipsoidCode) :
64  CoordinateSystem( ellipsoidSemiMajorAxis, ellipsoidFlattening ),
65  TranMerc_Origin_Long( centralMeridian ),
66  TranMerc_Origin_Lat( latitudeOfTrueScale ),
67  TranMerc_False_Easting( falseEasting ),
68  TranMerc_False_Northing( falseNorthing ),
69  TranMerc_Scale_Factor( scaleFactor ),
70  TranMerc_Delta_Northing( 20000000.0 ), // LSC-9451: Swapped values for Northing / Easting
71  TranMerc_Delta_Easting( 10000000.0 ) // LSC-9451: Swapped values for Northing / Easting
72 {
73  double TranMerc_b; // Semi-minor axis of ellipsoid, in meters
74  double invFlattening = 1.0 / ellipsoidFlattening;
75 
76  strcpy( ellipsCode, ellipsoidCode );
77 
78  if (ellipsCode[0] == '\0')
79  { // no ellipsoid code
80  throw CoordinateConversionException ( ErrorMessages::invalidEllipsoidCode );
81  }
82  if (ellipsoidSemiMajorAxis <= 0.0)
83  { // Semi-major axis must be greater than zero
84  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
85  }
86  if ( invFlattening < 150 )
87  {
88  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
89  }
90  if ((latitudeOfTrueScale < -PI_OVER_2) || (latitudeOfTrueScale > PI_OVER_2))
91  { // latitudeOfTrueScale out of range
92  throw CoordinateConversionException( ErrorMessages::originLatitude );
93  }
94  if ((centralMeridian < -PI) || (centralMeridian > (2*PI)))
95  { // centralMeridian out of range
96  throw CoordinateConversionException( ErrorMessages::centralMeridian );
97  }
98  if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
99  {
100  throw CoordinateConversionException( ErrorMessages::scaleFactor );
101  }
102 
103  if (TranMerc_Origin_Long > PI)
104  TranMerc_Origin_Long -= (2*PI);
105 
106  // Eccentricity
107  TranMerc_eps = sqrt( 2 * flattening - flattening * flattening );
108 
109  double n1, R4oa;
110  //added ellipsoid code as part of DR30125
111  generateCoefficients(
112  invFlattening, n1, TranMerc_aCoeff, TranMerc_bCoeff, R4oa, ellipsCode );
113 
114  TranMerc_K0R4 = R4oa * TranMerc_Scale_Factor * ellipsoidSemiMajorAxis;
115  TranMerc_K0R4inv = 1.0 / TranMerc_K0R4;
116 }
117 
118 
119 TransverseMercator::TransverseMercator( const TransverseMercator &tm )
120 {
121  *this = tm;
122 }
123 
124 
125 TransverseMercator::~TransverseMercator()
126 {
127 }
128 
129 TransverseMercator& TransverseMercator::operator=(
130  const TransverseMercator &tm )
131 {
132  if( this != &tm )
133  {
134  semiMajorAxis = tm.semiMajorAxis;
135  flattening = tm.flattening;
136  TranMerc_eps = tm.TranMerc_eps;
137 
138  TranMerc_K0R4 = tm.TranMerc_K0R4;
139  TranMerc_K0R4inv = tm.TranMerc_K0R4inv;
140 
141  for( int i = 0; i < MAX_TERMS; i++ )
142  {
143  TranMerc_aCoeff[i] = tm.TranMerc_aCoeff[i];
144  TranMerc_bCoeff[i] = tm.TranMerc_bCoeff[i];
145  }
146 
147  TranMerc_Origin_Long = tm.TranMerc_Origin_Long;
148  TranMerc_Origin_Lat = tm.TranMerc_Origin_Lat;
149  TranMerc_False_Easting = tm.TranMerc_False_Easting;
150  TranMerc_False_Northing = tm.TranMerc_False_Northing;
151  TranMerc_Scale_Factor = tm.TranMerc_Scale_Factor;
152 
153  TranMerc_Delta_Easting = tm.TranMerc_Delta_Easting;
154  TranMerc_Delta_Northing = tm.TranMerc_Delta_Northing;
155  }
156 
157  return *this;
158 }
159 
160 
161 MSP::CCS::MapProjection5Parameters* TransverseMercator::getParameters() const
162 {
163  return new MapProjection5Parameters(
164  CoordinateType::transverseMercator,
165  TranMerc_Origin_Long, TranMerc_Origin_Lat, TranMerc_Scale_Factor,
166  TranMerc_False_Easting, TranMerc_False_Northing );
167 }
168 
169 
170 MSP::CCS::MapProjectionCoordinates* TransverseMercator::convertFromGeodetic(
171  MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
172 {
173  double longitude = geodeticCoordinates->longitude();
174  double latitude = geodeticCoordinates->latitude();
175 
176  if (longitude > PI)
177  longitude -= (2 * PI);
178  if (longitude < -PI)
179  longitude += (2 * PI);
180 
181  // Convert longitude (Greenwhich) to longitude from the central meridian
182  // (-Pi, Pi] equivalent needed for checkLatLon.
183  // Compute its cosine and sine.
184  double lambda = longitude - TranMerc_Origin_Long;
185  if (lambda > PI)
186  lambda -= (2 * PI);
187  if (lambda < -PI)
188  lambda += (2 * PI);
189  checkLatLon( latitude, lambda );
190 
191  double easting, northing;
192  latLonToNorthingEasting( latitude, longitude, northing, easting );
193 
194  // The origin may move form (0,0) and this is represented by
195  // a change in the false Northing/Easting values.
196  double falseEasting, falseNorthing;
197  latLonToNorthingEasting(
198  TranMerc_Origin_Lat, TranMerc_Origin_Long, falseNorthing, falseEasting );
199 
200  easting += TranMerc_False_Easting - falseEasting;
201  northing += TranMerc_False_Northing - falseNorthing;
202 
203  char warning[256] = "";
204  warning[0] = '\0';
205  double invFlattening = 1.0 / flattening;
206  if( invFlattening < 290.0 || invFlattening > 301.0 )
207  strcat( warning,
208  "Eccentricity is outside range that algorithm accuracy has been tested." );
209 
210  return new MapProjectionCoordinates(
211  CoordinateType::transverseMercator, warning, easting, northing );
212 }
213 
214 
215 void TransverseMercator::latLonToNorthingEasting(
216  const double &latitude,
217  const double &longitude,
218  double &northing,
219  double &easting )
220 {
221  // Convert longitude (Greenwhich) to longitude from the central meridian
222  // (-Pi, Pi] equivalent needed for checkLatLon.
223  // Compute its cosine and sine.
224  double lambda = longitude - TranMerc_Origin_Long;
225  if (lambda > PI)
226  lambda -= (2 * PI);
227  if (lambda < -PI)
228  lambda += (2 * PI);
229  checkLatLon( latitude, lambda );
230 
231  double cosLam = cos(lambda);
232  double sinLam = sin(lambda);
233  double cosPhi = cos(latitude);
234  double sinPhi = sin(latitude);
235 
236  double P, part1, part2, denom, cosChi, sinChi;
237  double U, V;
238  double c2ku[MAX_TERMS], s2ku[MAX_TERMS];
239  double c2kv[MAX_TERMS], s2kv[MAX_TERMS];
240 
241  // Ellipsoid to sphere
242  // --------- -- ------
243 
244  // Convert geodetic latitude, Phi, to conformal latitude, Chi
245  // Only the cosine and sine of Chi are actually needed.
246  P = exp(TranMerc_eps * aTanH(TranMerc_eps * sinPhi));
247  part1 = (1 + sinPhi) / P;
248  part2 = (1 - sinPhi) * P;
249  denom = part1 + part2;
250  cosChi = 2 * cosPhi / denom;
251  sinChi = (part1 - part2) / denom;
252 
253  // Sphere to first plane
254  // ------ -- ----- -----
255 
256  // Apply spherical theory of transverse Mercator to get (u,v) coord.s
257  U = aTanH(cosChi * sinLam);
258  V = atan2(sinChi, cosChi * cosLam);
259 
260  // Use trig identities to compute cosh(2kU), sinh(2kU), cos(2kV), sin(2kV)
261  computeHyperbolicSeries( 2.0 * U, c2ku, s2ku );
262  computeTrigSeries( 2.0 * V, c2kv, s2kv );
263 
264  // First plane to second plane
265  // Accumulate terms for X and Y
266  double xStar = 0;
267  double yStar = 0;
268 
269  for (int k = N_TERMS - 1; k >= 0; k--)
270  {
271  xStar += TranMerc_aCoeff[k] * s2ku[k] * c2kv[k];
272  yStar += TranMerc_aCoeff[k] * c2ku[k] * s2kv[k];
273  }
274 
275  xStar += U;
276  yStar += V;
277 
278  // Apply isoperimetric radius, scale adjustment, and offsets
279  easting = (TranMerc_K0R4 * xStar);
280  northing = (TranMerc_K0R4 * yStar);
281 }
282 
283 
284 MSP::CCS::GeodeticCoordinates* TransverseMercator::convertToGeodetic(
285  MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
286 {
287  double easting = mapProjectionCoordinates->easting();
288  double northing = mapProjectionCoordinates->northing();
289 
290  if ( (easting < (TranMerc_False_Easting - TranMerc_Delta_Easting))
291  ||(easting > (TranMerc_False_Easting + TranMerc_Delta_Easting)))
292  { // easting out of range
293  throw CoordinateConversionException( ErrorMessages::easting );
294  }
295 
296  if ( (northing < (TranMerc_False_Northing - TranMerc_Delta_Northing))
297  || (northing > (TranMerc_False_Northing + TranMerc_Delta_Northing)))
298  { // northing out of range
299  throw CoordinateConversionException( ErrorMessages::northing );
300  }
301 
302  double longitude, latitude;
303  // The origin may move form (0,0) and this is represented by
304  // a change in the false Northing/Easting values.
305  double falseEasting, falseNorthing;
306  latLonToNorthingEasting(
307  TranMerc_Origin_Lat, TranMerc_Origin_Long, falseNorthing, falseEasting );
308 
309  easting -= (TranMerc_False_Easting - falseEasting);
310  northing -= (TranMerc_False_Northing - falseNorthing);
311 
312  northingEastingToLatLon( northing, easting, latitude, longitude );
313 
314  longitude = (longitude > PI) ? longitude - (2 * PI): longitude;
315  longitude = (longitude <= -PI) ? longitude + (2 * PI): longitude;
316 
317  if(fabs(latitude) > (90.0 * PI / 180.0))
318  {
319  throw CoordinateConversionException( ErrorMessages::northing );
320  }
321  if((longitude) > (PI))
322  {
323  longitude -= (2 * PI);
324  if(fabs(longitude) > PI)
325  throw CoordinateConversionException( ErrorMessages::easting );
326  }
327  else if((longitude) < (-PI))
328  {
329  longitude += (2 * PI);
330  if(fabs(longitude) > PI)
331  throw CoordinateConversionException( ErrorMessages::easting );
332  }
333 
334  char warning[256];
335  warning[0] = '\0';
336  double invFlattening = 1.0 / flattening;
337  if( invFlattening < 290.0 || invFlattening > 301.0 )
338  strcat( warning,
339  "Eccentricity is outside range that algorithm accuracy has been tested." );
340 
341  return new GeodeticCoordinates(
342  CoordinateType::geodetic, warning, longitude, latitude );
343 }
344 
345 void TransverseMercator::northingEastingToLatLon(
346  const double &northing,
347  const double &easting,
348  double &latitude,
349  double &longitude )
350 {
351  double c2kx[MAX_TERMS], s2kx[MAX_TERMS], c2ky[MAX_TERMS], s2ky[MAX_TERMS];
352  double U, V;
353  double lambda;
354  double sinChi;
355 
356  // Undo offsets, scale change, and factor R4
357  // ---- ------- ----- ------ --- ------ --
358  double xStar = TranMerc_K0R4inv * (easting);
359  double yStar = TranMerc_K0R4inv * (northing);
360 
361  // Use trig identities to compute cosh(2kU), sinh(2kU), cos(2kV), sin(2kV)
362  computeHyperbolicSeries( 2.0 * xStar, c2kx, s2kx );
363  computeTrigSeries( 2.0 * yStar, c2ky, s2ky );
364 
365  // Second plane (x*, y*) to first plane (u, v)
366  // ------ ----- -------- -- ----- ----- ------
367  U = 0;
368  V = 0;
369 
370  for (int k = N_TERMS - 1; k >= 0; k--)
371  {
372  U += TranMerc_bCoeff[k] * s2kx[k] * c2ky[k];
373  V += TranMerc_bCoeff[k] * c2kx[k] * s2ky[k];
374  }
375 
376  U += xStar;
377  V += yStar;
378 
379  // First plane to sphere
380  // ----- ----- -- ------
381  double coshU = cosh(U);
382  double sinhU = sinh(U);
383  double cosV = cos(V);
384  double sinV = sin(V);
385 
386  // Longitude from central meridian
387  if ((fabs(cosV) < 10E-12) && (fabs(coshU) < 10E-12))
388  lambda = 0;
389  else
390  lambda = atan2(sinhU, cosV);
391 
392  // Conformal latitude
393  sinChi = sinV / coshU;
394  latitude = geodeticLat( sinChi, TranMerc_eps );
395 
396  // Longitude from Greenwich
397  // -------- ---- ---------
398  longitude = TranMerc_Origin_Long + lambda;
399 }
400 
401 // PRIVATE FUNCTIONS
402 // added ellipsoidCode to function arguments as part of DR30125
403 void TransverseMercator::generateCoefficients(
404  double invfla,
405  double &n1,
406  double aCoeff[MAX_TERMS],
407  double bCoeff[MAX_TERMS],
408  double &R4oa,
409  char *ellipsoidCode)
410 {
411  /* Generate Coefficients for Transverse Mercator algorithms
412  ===----- ===---------
413  Algorithm developed by: C. Rollins April 18, 2006
414 
415  INPUT
416  -----
417  invfla Inverse flattening (reciprocal flattening)
418 
419  OUTPUT
420  ------
421  n1 Helmert's "n"
422  aCoeff Coefficients for omega as a trig series in chi
423  bBoeff Coefficients for chi as a trig series in omega
424  R4oa Ratio "R4 over a", i.e. R4/a
425 
426  EXPLANATIONS
427  ------------
428  omega is rectifying latitude
429  chi is conformal latitude
430  psi is geocentric latitude
431  phi is geodetic latitude, commonly, "the latitude"
432  R4 is the meridional isoperimetric radius
433  "a" is the semi-major axis of the ellipsoid
434  "b" is the semi-minor axis of the ellipsoid
435  Helmert's n = (a - b)/(a + b)
436 
437  This calculation depends only on the shape of the ellipsoid and are
438  independent of the ellipsoid size.
439 
440  The array Acoeff(8) stores eight coefficients corresponding
441  to k = 2, 4, 6, 8, 10, 12, 14, 16 in the notation "a sub k".
442  Likewise Bcoeff(8) etc.
443 */
444 
445  double n2, n3, n4, n5, n6, n7, n8, n9, n10, coeff;
446 
447  n1 = 1.0 / (2*invfla - 1.0);
448 
449  n2 = n1 * n1;
450  n3 = n2 * n1;
451  n4 = n3 * n1;
452  n5 = n4 * n1;
453  n6 = n5 * n1;
454  n7 = n6 * n1;
455  n8 = n7 * n1;
456  n9 = n8 * n1;
457  n10 = n9 * n1;
458 
459  // checks ellipsoid code and assigns values for corresponding coefficients.
460  // Uses default computation if ellipsoid code isn't found. This will be
461  // for user defined ellipsoids.
462 
463  if (( strcmp( ellipsoidCode, "AA") == 0) || (strcmp(ellipsoidCode, "AM") == 0))
464  {
465  aCoeff[0] = 8.3474517669594013740e-04;
466  aCoeff[1] = 7.554352936725572895e-07;
467  aCoeff[2] = 1.18487391005135489e-09;
468  aCoeff[3] = 2.3946872955703565e-12;
469  aCoeff[4] = 5.610633978440270e-15;
470  aCoeff[5] = 1.44858956458553e-17;
471  aCoeff[6] = 0.0; // LSC-16618
472  aCoeff[7] = 0.0; // LSC-16618
473 
474  bCoeff[0] = -8.3474551646761162264e-04;
475  bCoeff[1] = -5.863630361809676570e-08;
476  bCoeff[2] = -1.65562038746920803e-10;
477  bCoeff[3] = -2.1340335537652749e-13;
478  bCoeff[4] = -3.720760760132477e-16;
479  bCoeff[5] = -7.08304328877781e-19;
480  bCoeff[6] = 0.0; // LSC-16618
481  bCoeff[7] = 0.0; // LSC-16618
482  }
483  else if (( strcmp( ellipsoidCode, "EA") == 0) || ( strcmp( ellipsoidCode, "EB") == 0) ||
484  ( strcmp( ellipsoidCode, "EC") == 0) || ( strcmp( ellipsoidCode, "ED") == 0) ||
485  ( strcmp( ellipsoidCode, "EE") == 0))
486  {
487  aCoeff[0] = 8.3064943111192510534e-04;
488  aCoeff[1] = 7.480375027595025021e-07;
489  aCoeff[2] = 1.16750772278215999e-09;
490  aCoeff[3] = 2.3479972304395461e-12;
491  aCoeff[4] = 5.474212231879573e-15;
492  aCoeff[5] = 1.40642257446745e-17;
493  aCoeff[6] = 0.0; // LSC-16618
494  aCoeff[7] = 0.0; // LSC-16618
495 
496  bCoeff[0] = -8.3064976590443772201e-04;
497  bCoeff[1] = -5.805953517555717859e-08;
498  bCoeff[2] = -1.63133251663416522e-10;
499  bCoeff[3] = -2.0923797199593389e-13;
500  bCoeff[4] = -3.630200927775259e-16;
501  bCoeff[5] = -6.87666654919219e-19;
502  bCoeff[6] = 0.0; // LSC-16618
503  bCoeff[7] = 0.0; // LSC-16618
504  }
505  else if (( strcmp( ellipsoidCode, "BN") == 0) || ( strcmp( ellipsoidCode, "BR") == 0))
506  {
507  aCoeff[0] = 8.3522527226849818552e-04;
508  aCoeff[1] = 7.563048340614894422e-07;
509  aCoeff[2] = 1.18692075307408346e-09;
510  aCoeff[3] = 2.4002054791393298e-12;
511  aCoeff[4] = 5.626801597980756e-15;
512  aCoeff[5] = 1.45360057224474e-17;
513  aCoeff[6] = 0.0; // LSC-16618
514  aCoeff[7] = 0.0; // LSC-16618
515 
516  bCoeff[0] = -8.3522561262703079182e-04;
517  bCoeff[1] = -5.870409978661008580e-08;
518  bCoeff[2] = -1.65848307463131468e-10;
519  bCoeff[3] = -2.1389565927064571e-13;
520  bCoeff[4] = -3.731493368666479e-16;
521  bCoeff[5] = -7.10756898071999e-19;
522  bCoeff[6] = 0.0; // LSC-16618
523  bCoeff[7] = 0.0; // LSC-16618
524  }
525  else if (( strcmp( ellipsoidCode, "KA") == 0) || ( strcmp( ellipsoidCode, "HE") == 0) ||
526  ( strcmp( ellipsoidCode, "FA") == 0))
527  {
528  aCoeff[0] = 8.3761175713442343106e-04;
529  aCoeff[1] = 7.606346200814720197e-07;
530  aCoeff[2] = 1.19713032035541037e-09;
531  aCoeff[3] = 2.4277772986483520e-12;
532  aCoeff[4] = 5.707722772225013e-15;
533  aCoeff[5] = 1.47872454335773e-17;
534  aCoeff[6] = 0.0; // LSC-16618
535  aCoeff[7] = 0.0; // LSC-16618
536 
537  bCoeff[0] = -8.3761210042019176501e-04;
538  bCoeff[1] = -5.904169154078546237e-08;
539  bCoeff[2] = -1.67276212891429215e-10;
540  bCoeff[3] = -2.1635549847939549e-13;
541  bCoeff[4] = -3.785212121016612e-16;
542  bCoeff[5] = -7.23053625983667e-19;
543  bCoeff[6] = 0.0; // LSC-16618
544  bCoeff[7] = 0.0; // LSC-16618
545  }
546  else if ( strcmp( ellipsoidCode, "WD") == 0)
547  {
548  aCoeff[0] = 8.3772481044362217923e-04;
549  aCoeff[1] = 7.608400388863560936e-07;
550  aCoeff[2] = 1.19761541904924067e-09;
551  aCoeff[3] = 2.4290893081322466e-12;
552  aCoeff[4] = 5.711579173743133e-15;
553  aCoeff[5] = 1.47992364667635e-17;
554  aCoeff[6] = 0.0; // LSC-16618
555  aCoeff[7] = 0.0; // LSC-16618
556 
557  bCoeff[0] = -8.3772515386847544554e-04;
558  bCoeff[1] = -5.905770828762463028e-08;
559  bCoeff[2] = -1.67344058948464124e-10;
560  bCoeff[3] = -2.1647255130188214e-13;
561  bCoeff[4] = -3.787772179729998e-16;
562  bCoeff[5] = -7.23640523525528e-19;
563  bCoeff[6] = 0.0; // LSC-16618
564  bCoeff[7] = 0.0; // LSC-16618
565  }
566  else if ( strcmp( ellipsoidCode, "WE") == 0)
567  {
568  aCoeff[0] = 8.3773182062446983032e-04;
569  aCoeff[1] = 7.608527773572489156e-07;
570  aCoeff[2] = 1.19764550324249210e-09;
571  aCoeff[3] = 2.4291706803973131e-12;
572  aCoeff[4] = 5.711818369154105e-15;
573  aCoeff[5] = 1.47999802705262e-17;
574  aCoeff[6] = 0.0; // LSC-16618
575  aCoeff[7] = 0.0; // LSC-16618
576 
577  bCoeff[0] = -8.3773216405794867707e-04;
578  bCoeff[1] = -5.905870152220365181e-08;
579  bCoeff[2] = -1.67348266534382493e-10;
580  bCoeff[3] = -2.1647981104903862e-13;
581  bCoeff[4] = -3.787930968839601e-16;
582  bCoeff[5] = -7.23676928796690e-19;
583  bCoeff[6] = 0.0; // LSC-16618
584  bCoeff[7] = 0.0; // LSC-16618
585  }
586  else if ( strcmp( ellipsoidCode, "RF") == 0)
587  {
588  aCoeff[0] = 8.3773182472855134012e-04;
589  aCoeff[1] = 7.608527848149655006e-07;
590  aCoeff[2] = 1.19764552085530681e-09;
591  aCoeff[3] = 2.4291707280369697e-12;
592  aCoeff[4] = 5.711818509192422e-15;
593  aCoeff[5] = 1.47999807059922e-17;
594  aCoeff[6] = 0.0; // LSC-16618
595  aCoeff[7] = 0.0; // LSC-16618
596 
597  bCoeff[0] = -8.3773216816203523672e-04;
598  bCoeff[1] = -5.905870210369121594e-08;
599  bCoeff[2] = -1.67348268997717031e-10;
600  bCoeff[3] = -2.1647981529928124e-13;
601  bCoeff[4] = -3.787931061803592e-16;
602  bCoeff[5] = -7.23676950110361e-19;
603  bCoeff[6] = 0.0; // LSC-16618
604  bCoeff[7] = 0.0; // LSC-16618
605  }
606  else if (( strcmp( ellipsoidCode, "SA") == 0) || ( strcmp( ellipsoidCode, "AN") == 0))
607  {
608  aCoeff[0] = 8.3775209887947194075e-04;
609  aCoeff[1] = 7.608896263599627157e-07;
610  aCoeff[2] = 1.19773253021831769e-09;
611  aCoeff[3] = 2.4294060763606098e-12;
612  aCoeff[4] = 5.712510331613028e-15;
613  aCoeff[5] = 1.48021320370432e-17;
614  aCoeff[6] = 0.0; // LSC-16618
615  aCoeff[7] = 0.0; // LSC-16618
616 
617  bCoeff[0] = -8.3775244233790270051e-04;
618  bCoeff[1] = -5.906157468586898015e-08;
619  bCoeff[2] = -1.67360438158764851e-10;
620  bCoeff[3] = -2.1650081225048788e-13;
621  bCoeff[4] = -3.788390325953455e-16;
622  bCoeff[5] = -7.23782246429908e-19;
623  bCoeff[6] = 0.0; // LSC-16618
624  bCoeff[7] = 0.0; // LSC-16618
625  }
626  else if ( strcmp( ellipsoidCode, "ID") == 0)
627  {
628  aCoeff[0] = 8.3776052087969078729e-04;
629  aCoeff[1] = 7.609049308144604484e-07;
630  aCoeff[2] = 1.19776867565343872e-09;
631  aCoeff[3] = 2.4295038464530901e-12;
632  aCoeff[4] = 5.712797738386076e-15;
633  aCoeff[5] = 1.48030257891140e-17;
634  aCoeff[6] = 0.0; // LSC-16618
635  aCoeff[7] = 0.0; // LSC-16618
636 
637  bCoeff[0] = -8.3776086434848497443e-04;
638  bCoeff[1] = -5.906276799395007586e-08;
639  bCoeff[2] = -1.67365493472742884e-10;
640  bCoeff[3] = -2.1650953495573773e-13;
641  bCoeff[4] = -3.788581120060625e-16;
642  bCoeff[5] = -7.23825990889693e-19;
643  bCoeff[6] = 0.0; // LSC-16618
644  bCoeff[7] = 0.0; // LSC-16618
645  }
646  else if (( strcmp( ellipsoidCode, "IN") == 0) || ( strcmp( ellipsoidCode, "HO") == 0))
647  {
648  aCoeff[0] = 8.4127599100356448089e-04;
649  aCoeff[1] = 7.673066923431950296e-07;
650  aCoeff[2] = 1.21291995794281190e-09;
651  aCoeff[3] = 2.4705731165688123e-12;
652  aCoeff[4] = 5.833780550286833e-15;
653  aCoeff[5] = 1.51800420867708e-17;
654  aCoeff[6] = 0.0; // LSC-16618
655  aCoeff[7] = 0.0; // LSC-16618
656 
657  bCoeff[0] = -8.4127633881644851945e-04;
658  bCoeff[1] = -5.956193574768780571e-08;
659  bCoeff[2] = -1.69484573979154433e-10;
660  bCoeff[3] = -2.2017363465021880e-13;
661  bCoeff[4] = -3.868896221495780e-16;
662  bCoeff[5] = -7.42279219864412e-19;
663  bCoeff[6] = 0.0; // LSC-16618
664  bCoeff[7] = 0.0; // LSC-16618
665  }
666  else if ( strcmp( ellipsoidCode, "WO") == 0)
667  {
668  aCoeff[0] = 8.4411652150600103279e-04;
669  aCoeff[1] = 7.724989750172583427e-07;
670  aCoeff[2] = 1.22525529789972041e-09;
671  aCoeff[3] = 2.5041361775549209e-12;
672  aCoeff[4] = 5.933026083631383e-15;
673  aCoeff[5] = 1.54904908794521e-17;
674  aCoeff[6] = 0.0; // LSC-16618
675  aCoeff[7] = 0.0; // LSC-16618
676 
677  bCoeff[0] = -8.4411687285559594196e-04;
678  bCoeff[1] = -5.996681687064322548e-08;
679  bCoeff[2] = -1.71209836918814857e-10;
680  bCoeff[3] = -2.2316811233502163e-13;
681  bCoeff[4] = -3.934782433323038e-16;
682  bCoeff[5] = -7.57474665717687e-19;
683  bCoeff[6] = 0.0; // LSC-16618
684  bCoeff[7] = 0.0; // LSC-16618
685  }
686  else if ( strcmp( ellipsoidCode, "CC") == 0)
687  {
688  aCoeff[0] = 8.4703742793654652315e-04;
689  aCoeff[1] = 7.778564517658115212e-07;
690  aCoeff[2] = 1.23802665917879731e-09;
691  aCoeff[3] = 2.5390045684252928e-12;
692  aCoeff[4] = 6.036484469753319e-15;
693  aCoeff[5] = 1.58152259295850e-17;
694  aCoeff[6] = 0.0; // LSC-16618
695  aCoeff[7] = 0.0; // LSC-16618
696 
697  bCoeff[0] = -8.4703778294785813001e-04;
698  bCoeff[1] = -6.038459874600183555e-08;
699  bCoeff[2] = -1.72996106059227725e-10;
700  bCoeff[3] = -2.2627911073545072e-13;
701  bCoeff[4] = -4.003466873888566e-16;
702  bCoeff[5] = -7.73369749524777e-19;
703  bCoeff[6] = 0.0; // LSC-16618
704  bCoeff[7] = 0.0; // LSC-16618
705  }
706  else if ( strcmp( ellipsoidCode, "CG") == 0)
707  {
708  aCoeff[0] = 8.5140099460764136776e-04;
709  aCoeff[1] = 7.858945456038187774e-07;
710  aCoeff[2] = 1.25727085106103462e-09;
711  aCoeff[3] = 2.5917718627340128e-12;
712  aCoeff[4] = 6.193726879043722e-15;
713  aCoeff[5] = 1.63109098395549e-17;
714  aCoeff[6] = 0.0; // LSC-16618
715  aCoeff[7] = 0.0; // LSC-16618
716 
717  bCoeff[0] = -8.5140135513650084564e-04;
718  bCoeff[1] = -6.101145475063033499e-08;
719  bCoeff[2] = -1.75687742410879760e-10;
720  bCoeff[3] = -2.3098718484594067e-13;
721  bCoeff[4] = -4.107860472919190e-16;
722  bCoeff[5] = -7.97633133452512e-19;
723  bCoeff[6] = 0.0; // LSC-16618
724  bCoeff[7] = 0.0; // LSC-16618
725  }
726  else if ( strcmp( ellipsoidCode, "CD") == 0)
727  {
728  aCoeff[0] = 8.5140395445291970541e-04;
729  aCoeff[1] = 7.859000119464140978e-07;
730  aCoeff[2] = 1.25728397182445579e-09;
731  aCoeff[3] = 2.5918079321459932e-12;
732  aCoeff[4] = 6.193834639108787e-15;
733  aCoeff[5] = 1.63112504092335e-17;
734  aCoeff[6] = 0.0; // LSC-16618
735  aCoeff[7] = 0.0; // LSC-16618
736 
737  bCoeff[0] = -8.5140431498554106268e-04;
738  bCoeff[1] = -6.101188106187092184e-08;
739  bCoeff[2] = -1.75689577596504470e-10;
740  bCoeff[3] = -2.3099040312610703e-13;
741  bCoeff[4] = -4.107932016207395e-16;
742  bCoeff[5] = -7.97649804397335e-19;
743  bCoeff[6] = 0.0; // LSC-16618
744  bCoeff[7] = 0.0; // LSC-16618
745  }
746  else
747  {
748 
749  // computation below is for user defined ellipsoids
750  // Computation of coefficient a2
751  coeff = 0.0;
752  coeff += (-18975107.0) * n8 / 50803200.0;
753  coeff += (72161.0) * n7 / 387072.0;
754  coeff += (7891.0) * n6 / 37800.0;
755  coeff += (-127.0) * n5 / 288.0;
756  coeff += (41.0) * n4 / 180.0;
757  coeff += (5.0) * n3 / 16.0;
758  coeff += (-2.0) * n2 / 3.0;
759  coeff += (1.0) * n1 / 2.0;
760 
761  aCoeff[0] = coeff;
762 
763  // Computation of coefficient a4
764  coeff = 0.0;
765  coeff += (148003883.0) * n8 / 174182400.0;
766  coeff += (13769.0) * n7 / 28800.0;
767  coeff += (-1983433.0) * n6 / 1935360.0;
768  coeff += (281.0) * n5 / 630.0;
769  coeff += (557.0) * n4 / 1440.0;
770  coeff += (-3.0) * n3 / 5.0;
771  coeff += (13.0) * n2 / 48.0;
772 
773  aCoeff[1] = coeff;
774 
775  // Computation of coefficient a6
776  coeff = 0.0;
777  coeff += (79682431.0) * n8 / 79833600.0;
778  coeff += (-67102379.0) * n7 / 29030400.0;
779  coeff += (167603.0) * n6 / 181440.0;
780  coeff += (15061.0) * n5 / 26880.0;
781  coeff += (-103.0) * n4 / 140.0;
782  coeff += (61.0) * n3 / 240.0;
783 
784  aCoeff[2] = coeff;
785 
786  // Computation of coefficient a8
787  coeff = 0.0;
788  coeff += (-40176129013.0) * n8 / 7664025600.0;
789  coeff += (97445.0) * n7 / 49896.0;
790  coeff += (6601661.0) * n6 / 7257600.0;
791  coeff += (-179.0) * n5 / 168.0;
792  coeff += (49561.0) * n4 / 161280.0;
793 
794  aCoeff[3] = coeff;
795 
796  // Computation of coefficient a10
797  coeff = 0.0;
798  coeff += (2605413599.0) * n8 / 622702080.0;
799  coeff += (14644087.0) * n7 / 9123840.0;
800  coeff += (-3418889.0) * n6 / 1995840.0;
801  coeff += (34729.0) * n5 / 80640.0;
802 
803  aCoeff[4] = coeff;
804 
805  // Computation of coefficient a12
806  coeff = 0.0;
807  coeff += (175214326799.0) * n8 / 58118860800.0;
808  coeff += (-30705481.0) * n7 / 10378368.0;
809  coeff += (212378941.0) * n6 / 319334400.0;
810 
811  aCoeff[5] = coeff;
812 
813  // Computation of coefficient a14
814  coeff = 0.0;
815  coeff += (-16759934899.0) * n8 / 3113510400.0;
816  coeff += (1522256789.0) * n7 / 1383782400.0;
817 
818  aCoeff[6] = coeff;
819 
820  // Computation of coefficient a16
821  coeff = 0.0;
822  coeff += (1424729850961.0) * n8 / 743921418240.0;
823 
824  aCoeff[7] = coeff;
825 
826  // Computation of coefficient b2
827  coeff = 0.0;
828  coeff += (-7944359.0) * n8 / 67737600.0;
829  coeff += (5406467.0) * n7 / 38707200.0;
830  coeff += (-96199.0) * n6 / 604800.0;
831  coeff += (81.0) * n5 / 512.0;
832  coeff += (1.0) * n4 / 360.0;
833  coeff += (-37.0) * n3 / 96.0;
834  coeff += (2.0) * n2 / 3.0;
835  coeff += (-1.0) * n1 / 2.0;
836 
837  bCoeff[0] = coeff;
838 
839  // Computation of coefficient b4
840  coeff = 0.0;
841  coeff += (-24749483.0) * n8 / 348364800.0;
842  coeff += (-51841.0) * n7 / 1209600.0;
843  coeff += (1118711.0) * n6 / 3870720.0;
844  coeff += (-46.0) * n5 / 105.0;
845  coeff += (437.0) * n4 / 1440.0;
846  coeff += (-1.0) * n3 / 15.0;
847  coeff += (-1.0) * n2 / 48.0;
848 
849  bCoeff[1] = coeff;
850 
851  // Computation of coefficient b6
852  coeff = 0.0;
853  coeff += (6457463.0) * n8 / 17740800.0;
854  coeff += (-9261899.0) * n7 / 58060800.0;
855  coeff += (-5569.0) * n6 / 90720.0;
856  coeff += (209.0) * n5 / 4480.0;
857  coeff += (37.0) * n4 / 840.0;
858  coeff += (-17.0) * n3 / 480.0;
859 
860  bCoeff[2] = coeff;
861 
862  // Computation of coefficient b8
863  coeff = 0.0;
864  coeff += (-324154477.0) * n8 / 7664025600.0;
865  coeff += (-466511.0) * n7 / 2494800.0;
866  coeff += (830251.0) * n6 / 7257600.0;
867  coeff += (11.0) * n5 / 504.0;
868  coeff += (-4397.0) * n4 / 161280.0;
869 
870  bCoeff[3] = coeff;
871 
872  // Computation of coefficient b10
873  coeff = 0.0;
874  coeff += (-22894433.0) * n8 / 124540416.0;
875  coeff += (8005831.0) * n7 / 63866880.0;
876  coeff += (108847.0) * n6 / 3991680.0;
877  coeff += (-4583.0) * n5 / 161280.0;
878 
879  bCoeff[4] = coeff;
880 
881  // Computation of coefficient b12
882  coeff = 0.0;
883  coeff += (2204645983.0) * n8 / 12915302400.0;
884  coeff += (16363163.0) * n7 / 518918400.0;
885  coeff += (-20648693.0) * n6 / 638668800.0;
886 
887  bCoeff[5] = coeff;
888 
889  // Computation of coefficient b14
890  coeff = 0.0;
891  coeff += (497323811.0) * n8 / 12454041600.0;
892  coeff += (-219941297.0) * n7 / 5535129600.0;
893 
894  bCoeff[6] = coeff;
895 
896  // Computation of coefficient b16
897  coeff = 0.0;
898  coeff += (-191773887257.0) * n8 / 3719607091200.0;
899 
900  bCoeff[7] = coeff;
901  }
902 
903  coeff = 0.0;
904  coeff += 49 * n10 / 65536.0;
905  coeff += 25 * n8 / 16384.0;
906  coeff += n6 / 256.0;
907  coeff += n4 / 64.0;
908  coeff += n2 / 4;
909  coeff += 1;
910  R4oa = coeff / (1 + n1);
911 }
912 
913 
914 void TransverseMercator::checkLatLon( double latitude, double deltaLon )
915 {
916  // test is based on distance from central meridian = deltaLon
917  if (deltaLon > PI)
918  deltaLon -= (2 * PI);
919  if (deltaLon < -PI)
920  deltaLon += (2 * PI);
921 
922  double testAngle = fabs( deltaLon );
923 
924  double delta = fabs( deltaLon - PI );
925  if( delta < testAngle )
926  testAngle = delta;
927 
928  delta = fabs( deltaLon + PI );
929  if( delta < testAngle )
930  testAngle = delta;
931 
932  // Away from the equator, is also valid
933  delta = PI_OVER_2 - latitude;
934  if( delta < testAngle )
935  testAngle = delta;
936 
937  delta = PI_OVER_2 + latitude;
938  if( delta < testAngle )
939  testAngle = delta;
940 
941  if( testAngle > MAX_DELTA_LONG )
942  {
943  throw CoordinateConversionException( ErrorMessages::longitude );
944  }
945 }
946 
947 
948 double TransverseMercator::aTanH(double x)
949 {
950  return(0.5 * log((1 + x) / (1 - x)));
951 }
952 
953 
954 double TransverseMercator::geodeticLat(
955  double sinChi,
956  double e )
957 {
958  double p;
959  double pSq;
960  double s_old = 1.0e99;
961  double s = sinChi;
962  double onePlusSinChi = 1.0+sinChi;
963  double oneMinusSinChi = 1.0-sinChi;
964 
965  for( int n = 0; n < 30; n++ )
966  {
967  p = exp( e * aTanH( e * s ) );
968  pSq = p * p;
969  s = ( onePlusSinChi * pSq - oneMinusSinChi )
970  /( onePlusSinChi * pSq + oneMinusSinChi );
971 
972  if( fabs( s - s_old ) < 1.0e-12 )
973  {
974  break;
975  }
976  s_old = s;
977  }
978  return asin(s);
979 }
980 
981 void TransverseMercator::computeHyperbolicSeries(
982  double twoX,
983  double c2kx[],
984  double s2kx[])
985 {
986  // Use trig identities to compute
987  // c2kx[k] = cosh(2kX), s2kx[k] = sinh(2kX) for k = 0 .. 8
988  c2kx[0] = cosh(twoX);
989  s2kx[0] = sinh(twoX);
990  c2kx[1] = 2.0 * c2kx[0] * c2kx[0] - 1.0;
991  s2kx[1] = 2.0 * c2kx[0] * s2kx[0];
992  c2kx[2] = c2kx[0] * c2kx[1] + s2kx[0] * s2kx[1];
993  s2kx[2] = c2kx[1] * s2kx[0] + c2kx[0] * s2kx[1];
994  c2kx[3] = 2.0 * c2kx[1] * c2kx[1] - 1.0;
995  s2kx[3] = 2.0 * c2kx[1] * s2kx[1];
996  c2kx[4] = c2kx[0] * c2kx[3] + s2kx[0] * s2kx[3];
997  s2kx[4] = c2kx[3] * s2kx[0] + c2kx[0] * s2kx[3];
998  c2kx[5] = 2.0 * c2kx[2] * c2kx[2] - 1.0;
999  s2kx[5] = 2.0 * c2kx[2] * s2kx[2];
1000  c2kx[6] = c2kx[0] * c2kx[5] + s2kx[0] * s2kx[5];
1001  s2kx[6] = c2kx[5] * s2kx[0] + c2kx[0] * s2kx[5];
1002  c2kx[7] = 2.0 * c2kx[3] * c2kx[3] - 1.0;
1003  s2kx[7] = 2.0 * c2kx[3] * s2kx[3];
1004 }
1005 
1006 void TransverseMercator::computeTrigSeries(
1007  double twoY,
1008  double c2ky[],
1009  double s2ky[])
1010 {
1011  // Use trig identities to compute
1012  // c2ky[k] = cos(2kY), s2ky[k] = sin(2kY) for k = 0 .. 8
1013  c2ky[0] = cos(twoY);
1014  s2ky[0] = sin(twoY);
1015  c2ky[1] = 2.0 * c2ky[0] * c2ky[0] - 1.0;
1016  s2ky[1] = 2.0 * c2ky[0] * s2ky[0];
1017  c2ky[2] = c2ky[1] * c2ky[0] - s2ky[1] * s2ky[0];
1018  s2ky[2] = c2ky[1] * s2ky[0] + c2ky[0] * s2ky[1];
1019  c2ky[3] = 2.0 * c2ky[1] * c2ky[1] - 1.0;
1020  s2ky[3] = 2.0 * c2ky[1] * s2ky[1];
1021  c2ky[4] = c2ky[3] * c2ky[0] - s2ky[3] * s2ky[0];
1022  s2ky[4] = c2ky[3] * s2ky[0] + c2ky[0] * s2ky[3];
1023  c2ky[5] = 2.0 * c2ky[2] * c2ky[2] - 1.0;
1024  s2ky[5] = 2.0 * c2ky[2] * s2ky[2];
1025  c2ky[6] = c2ky[5] * c2ky[0] - s2ky[5] * s2ky[0];
1026  s2ky[6] = c2ky[5] * s2ky[0] + c2ky[0] * s2ky[5];
1027  c2ky[7] = 2.0 * c2ky[3] * c2ky[3] - 1.0;
1028  s2ky[7] = 2.0 * c2ky[3] * s2ky[3];
1029 }
1030 
1031 // CLASSIFICATION: UNCLASSIFIED
MIN_SCALE_FACTOR
#define MIN_SCALE_FACTOR
Definition: PolarStereographic.cpp:129
MapProjectionCoordinates.h
MSP::CCS::TransverseMercator::TransverseMercator
TransverseMercator(double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double latitudeOfTrueScale, double falseEasting, double falseNorthing, double scaleFactor, char *ellipsoidCode)
Definition: TransverseMercator.cpp:55
PI_OVER_2
#define PI_OVER_2
Definition: MGRS.cpp:165
MSP::CCS::ErrorMessages::ellipsoidFlattening
static const char * ellipsoidFlattening
Definition: ErrorMessages.h:47
CoordinateConversionException.h
MAX_SCALE_FACTOR
#define MAX_SCALE_FACTOR
Definition: PolarStereographic.cpp:130
MSP::CCS::TransverseMercator::getParameters
MapProjection5Parameters * getParameters() const
Definition: TransverseMercator.cpp:161
MSP::CCS::ErrorMessages::centralMeridian
static const char * centralMeridian
Definition: ErrorMessages.h:51
MSP::CCS::TransverseMercator::convertFromGeodetic
MSP::CCS::MapProjectionCoordinates * convertFromGeodetic(MSP::CCS::GeodeticCoordinates *geodeticCoordinates)
Definition: TransverseMercator.cpp:170
MSP::CCS::ErrorMessages::longitude
static const char * longitude
Definition: ErrorMessages.h:76
MAX_TERMS
#define MAX_TERMS
Definition: TransverseMercator.cpp:45
MapProjection5Parameters.h
MSP::CCS::ErrorMessages::semiMajorAxis
static const char * semiMajorAxis
Definition: ErrorMessages.h:46
MAX_DELTA_LONG
#define MAX_DELTA_LONG
Definition: TransverseMercator.cpp:50
PI
#define PI
Definition: MGRS.cpp:164
MSP::CCS::ErrorMessages::originLatitude
static const char * originLatitude
Definition: ErrorMessages.h:49
MSP::CCS::ErrorMessages::invalidEllipsoidCode
static const char * invalidEllipsoidCode
Definition: ErrorMessages.h:31
ErrorMessages.h
GeodeticCoordinates.h
MSP::CCS::ErrorMessages::northing
static const char * northing
Definition: ErrorMessages.h:78
MSP::CCS::TransverseMercator::operator=
TransverseMercator & operator=(const TransverseMercator &tm)
Definition: TransverseMercator.cpp:129
MSP::CCS::ErrorMessages::scaleFactor
static const char * scaleFactor
Definition: ErrorMessages.h:52
N_TERMS
#define N_TERMS
Definition: TransverseMercator.cpp:44
MSP::CCS::ErrorMessages::easting
static const char * easting
Definition: ErrorMessages.h:77
TransverseMercator.h
MSP::CCS::TransverseMercator::convertToGeodetic
MSP::CCS::GeodeticCoordinates * convertToGeodetic(MSP::CCS::MapProjectionCoordinates *mapProjectionCoordinates)
Definition: TransverseMercator.cpp:284
Generated on Tue Apr 12 2022 10:09:23 for GeographicTranslator by doxygen 1.8.5