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