UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
Mercator.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: MERCATOR
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Mercator projection coordinates
10  * (easting and northing in meters).
11  *
12  * ERROR HANDLING
13  *
14  * This component checks parameters for valid values. If an invalid value
15  * is found, the error code is combined with the current error code using
16  * the bitwise or. This combining allows multiple error codes to be
17  * returned. The possible error codes are:
18  *
19  * MERC_NO_ERROR : No errors occurred in function
20  * MERC_LAT_ERROR : Latitude outside of valid range
21  * (-89.5 to 89.5 degrees)
22  * MERC_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * MERC_EASTING_ERROR : Easting outside of valid range
25  * (False_Easting +/- ~20,500,000 m,
26  * depending on ellipsoid parameters
27  * and Origin_Latitude)
28  * MERC_NORTHING_ERROR : Northing outside of valid range
29  * (False_Northing +/- ~23,500,000 m,
30  * depending on ellipsoid parameters
31  * and Origin_Latitude)
32  * MERC_LAT_OF_TRUE_SCALE_ERROR : Latitude of true scale outside of valid range
33  * (-89.5 to 89.5 degrees)
34  * MERC_CENT_MER_ERROR : Central meridian outside of valid range
35  * (-180 to 360 degrees)
36  * MERC_A_ERROR : Semi-major axis less than or equal to zero
37  * MERC_INV_F_ERROR : Inverse flattening outside of valid range
38  * (250 to 350)
39  *
40  * REUSE NOTES
41  *
42  * MERCATOR is intended for reuse by any application that performs a
43  * Mercator projection or its inverse.
44  *
45  * REFERENCES
46  *
47  * Further information on MERCATOR can be found in the Reuse Manual.
48  *
49  * MERCATOR originated from : U.S. Army Topographic Engineering Center
50  * Geospatial Information Division
51  * 7701 Telegraph Road
52  * Alexandria, VA 22310-3864
53  *
54  * LICENSES
55  *
56  * None apply to this component.
57  *
58  * RESTRICTIONS
59  *
60  * MERCATOR has no restrictions.
61  *
62  * ENVIRONMENT
63  *
64  * MERCATOR was tested and certified in the following environments:
65  *
66  * 1. Solaris 2.5 with GCC, version 2.8.1
67  * 2. Windows 95 with MS Visual C++, version 6
68  *
69  * MODIFICATIONS
70  *
71  * Date Description
72  * ---- -----------
73  * 10-02-97 Original Code
74  * 03-06-07 Original C++ Code
75  *
76  */
77 
78 
79 /***************************************************************************/
80 /*
81  * INCLUDES
82  */
83 
84 #include <stdio.h>
85 #include <math.h>
86 #include "Mercator.h"
87 #include "MercatorStandardParallelParameters.h"
88 #include "MercatorScaleFactorParameters.h"
89 #include "MapProjectionCoordinates.h"
90 #include "GeodeticCoordinates.h"
91 #include "CoordinateConversionException.h"
92 #include "ErrorMessages.h"
93 
94 /*
95  * math.h - Standard C math library
96  * Mercator.h - Is for prototype error checking
97  * MapProjectionCoordinates.h - defines map projection coordinates
98  * GeodeticCoordinates.h - defines geodetic coordinates
99  * CoordinateConversionException.h - Exception handler
100  * ErrorMessages.h - Contains exception messages
101  */
102 
103 
104 using namespace MSP::CCS;
105 
106 
107 /***************************************************************************/
108 /* DEFINES
109  *
110  */
111 
112 const double PI = 3.14159265358979323e0; /* PI */
113 const double PI_OVER_2 = ( PI / 2.0e0);
114 const double TWO_PI = (2.0 * PI);
115 const double MAX_LAT = ( (PI * 89.5) / 180.0 ); /* 89.5 degrees in radians */
116 const double MIN_SCALE_FACTOR = 0.3;
117 const double MAX_SCALE_FACTOR = 3.0;
118 
119 
120 /************************************************************************/
121 /* FUNCTIONS
122  *
123  */
124 
125 Mercator::Mercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double standardParallel, double falseEasting, double falseNorthing, double* scaleFactor ) :
126  CoordinateSystem(),
127  coordinateType( CoordinateType::mercatorStandardParallel ),
128  Merc_e( 0.08181919084262188000 ),
129  Merc_es( 0.0066943799901413800 ),
130  Merc_Cent_Mer( 0.0 ),
131  Merc_Standard_Parallel( 0.0 ),
132  Merc_False_Easting( 0.0 ),
133  Merc_False_Northing( 0.0 ),
134  Merc_Scale_Factor( 1.0 ),
135  Merc_ab( 0.00335655146887969400 ),
136  Merc_bb( 0.00000657187271079536 ),
137  Merc_cb( 0.00000001764564338702 ),
138  Merc_db( 0.00000000005328478445 ),
139  Merc_Delta_Easting( 20237883.0 ),
140  Merc_Delta_Northing( 23421740.0 )
141 {
142 /*
143  * The constructor receives the ellipsoid parameters and
144  * Mercator (Standard Parallel) projection parameters as inputs, and sets the corresponding state
145  * variables. It calculates and returns the scale factor. If any errors occur,
146  * an exception is thrown with a description of the error.
147  *
148  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
149  * ellipsoidFlattening : Flattening of ellipsoid (input)
150  * centralMeridian : Longitude in radians at the center of (input)
151  * the projection
152  * standardParallel : Latitude in radians at which the (input)
153  * point scale factor is 1.0
154  * falseEasting : A coordinate value in meters assigned to the
155  * central meridian of the projection. (input)
156  * falseNorthing : A coordinate value in meters assigned to the
157  * origin latitude of the projection (input)
158  * scaleFactor : Multiplier which reduces distances in the
159  * projection to the actual distance on the
160  * ellipsoid (output)
161  */
162 
163  double es2; /* Eccentricity squared of ellipsoid to the second power */
164  double es3; /* Eccentricity squared of ellipsoid to the third power */
165  double es4; /* Eccentricity squared of ellipsoid to the fourth power */
166  double sin_olat; /* sin(Origin_Latitude), temp variable */
167  double inv_f = 1 / ellipsoidFlattening;
168 
169  if (ellipsoidSemiMajorAxis <= 0.0)
170  { /* Semi-major axis must be greater than zero */
171  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
172  }
173  if ((inv_f < 250) || (inv_f > 350))
174  { /* Inverse flattening must be between 250 and 350 */
175  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
176  }
177  if ((standardParallel < -MAX_LAT) || (standardParallel > MAX_LAT))
178  { /* latitude of true scale out of range */
179  throw CoordinateConversionException( ErrorMessages::originLatitude );
180  }
181  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
182  { /* central meridian out of range */
183  throw CoordinateConversionException( ErrorMessages::centralMeridian );
184  }
185 
186  semiMajorAxis = ellipsoidSemiMajorAxis;
187  flattening = ellipsoidFlattening;
188 
189  Merc_Standard_Parallel = standardParallel;
190  if (centralMeridian > PI)
191  centralMeridian -= TWO_PI;
192  Merc_Cent_Mer = centralMeridian;
193  Merc_False_Northing = falseNorthing;
194  Merc_False_Easting = falseEasting;
195 
196  Merc_es = 2 * flattening - flattening * flattening;
197  Merc_e = sqrt(Merc_es);
198  sin_olat = sin(Merc_Standard_Parallel);
199  Merc_Scale_Factor = cos(Merc_Standard_Parallel) / sqrt(1.e0 - Merc_es * sin_olat * sin_olat );
200 // Merc_Scale_Factor = 1.0 / ( sqrt(1.e0 - Merc_es * sin_olat * sin_olat)
201 // / cos(Merc_Standard_Parallel) );
202  es2 = Merc_es * Merc_es;
203  es3 = es2 * Merc_es;
204  es4 = es3 * Merc_es;
205  Merc_ab = Merc_es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0
206  + 13.e0 * es4 / 360.e0;
207  Merc_bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0
208  + 811.e0 * es4 / 11520.e0;
209  Merc_cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
210  Merc_db = 4279.e0 * es4 / 161280.e0;
211  *scaleFactor = Merc_Scale_Factor;
212 
213  /* Calculate the width of the bounding box */
214  /* Note: The width of the bounding box needs to be relative */
215  /* to a false origin of 0,0, so subtract the false easting */
216  /* and false northing values from the delta easting and */
217  /* delta northing values */
218  GeodeticCoordinates gcTemp( CoordinateType::geodetic, ( Merc_Cent_Mer + PI ), MAX_LAT );
219  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
220  Merc_Delta_Easting = tempCoordinates->easting();
221  Merc_Delta_Northing = tempCoordinates->northing();
222  delete tempCoordinates;
223 
224  if(Merc_False_Easting)
225  Merc_Delta_Easting -= Merc_False_Easting;
226  if (Merc_Delta_Easting < 0)
227  Merc_Delta_Easting = -Merc_Delta_Easting;
228  Merc_Delta_Easting *= 1.01;
229 
230  if(Merc_False_Northing)
231  Merc_Delta_Northing -= Merc_False_Northing;
232  if (Merc_Delta_Northing < 0)
233  Merc_Delta_Northing = -Merc_Delta_Northing;
234  Merc_Delta_Northing *= 1.01;
235 }
236 
237 
238 Mercator::Mercator( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing, double scaleFactor ) :
239  CoordinateSystem(),
240  coordinateType( CoordinateType::mercatorScaleFactor ),
241  Merc_e( 0.08181919084262188000 ),
242  Merc_es( 0.0066943799901413800 ),
243  Merc_Cent_Mer( 0.0 ),
244  Merc_Standard_Parallel( 0.0 ),
245  Merc_False_Easting( 0.0 ),
246  Merc_False_Northing( 0.0 ),
247  Merc_Scale_Factor( 1.0 ),
248  Merc_ab( 0.00335655146887969400 ),
249  Merc_bb( 0.00000657187271079536 ),
250  Merc_cb( 0.00000001764564338702 ),
251  Merc_db( 0.00000000005328478445 ),
252  Merc_Delta_Easting( 20237883.0 ),
253  Merc_Delta_Northing( 23421740.0 )
254 {
255 /*
256  * The constructor receives the ellipsoid parameters and
257  * Mercator (Scale Factor) projection parameters as inputs, and sets the corresponding state
258  * variables. It receives the scale factor as input. If any errors occur,
259  * an exception is thrown with a description of the error.
260  *
261  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
262  * ellipsoidFlattening : Flattening of ellipsoid (input)
263  * centralMeridian : Longitude in radians at the center of (input)
264  * the projection
265  * falseEasting : A coordinate value in meters assigned to the
266  * central meridian of the projection. (input)
267  * falseNorthing : A coordinate value in meters assigned to the
268  * origin latitude of the projection (input)
269  * scaleFactor : Multiplier which reduces distances in the
270  * projection to the actual distance on the
271  * ellipsoid (input)
272  */
273 
274  double es2; /* Eccentricity squared of ellipsoid to the second power */
275  double es3; /* Eccentricity squared of ellipsoid to the third power */
276  double es4; /* Eccentricity squared of ellipsoid to the fourth power */
277  double Merc_Scale_Factor2;
278  double inv_f = 1 / ellipsoidFlattening;
279 
280  if (ellipsoidSemiMajorAxis <= 0.0)
281  { /* Semi-major axis must be greater than zero */
282  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
283  }
284  if ((inv_f < 250) || (inv_f > 350))
285  { /* Inverse flattening must be between 250 and 350 */
286  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
287  }
288  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
289  { /* central meridian out of range */
290  throw CoordinateConversionException( ErrorMessages::centralMeridian );
291  }
292  if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
293  {
294  throw CoordinateConversionException( ErrorMessages::scaleFactor );
295  }
296 
297  semiMajorAxis = ellipsoidSemiMajorAxis;
298  flattening = ellipsoidFlattening;
299 
301  if (centralMeridian > PI)
302  centralMeridian -= TWO_PI;
303  Merc_Cent_Mer = centralMeridian;
304  Merc_False_Northing = falseNorthing;
305  Merc_False_Easting = falseEasting;
306  Merc_Scale_Factor = scaleFactor;
307 
308  Merc_es = 2 * flattening - flattening * flattening;
309  Merc_e = sqrt(Merc_es);
310 // sin_olat = sin(Merc_Standard_Parallel);
311 // Merc_Scale_Factor = cos(Merc_Standard_Parallel) / sqrt(1.e0 - Merc_es * sin_olat * sin_olat );
312 // Merc_Scale_Factor = 1.0 / ( sqrt(1.e0 - Merc_es * sin_olat * sin_olat)
313 // / cos(Merc_Standard_Parallel) );
314  es2 = Merc_es * Merc_es;
315  es3 = es2 * Merc_es;
316  es4 = es3 * Merc_es;
317  Merc_ab = Merc_es / 2.e0 + 5.e0 * es2 / 24.e0 + es3 / 12.e0
318  + 13.e0 * es4 / 360.e0;
319  Merc_bb = 7.e0 * es2 / 48.e0 + 29.e0 * es3 / 240.e0
320  + 811.e0 * es4 / 11520.e0;
321  Merc_cb = 7.e0 * es3 / 120.e0 + 81.e0 * es4 / 1120.e0;
322  Merc_db = 4279.e0 * es4 / 161280.e0;
323 // *scaleFactor = Merc_Scale_Factor;
324 
325  Merc_Scale_Factor2 = Merc_Scale_Factor * Merc_Scale_Factor;
326  Merc_Standard_Parallel = asin(sqrt((1 - Merc_Scale_Factor2)/(1 - Merc_Scale_Factor2 * Merc_es)));
327 
328  /* Calculate the width of the bounding box */
329  /* Note: The width of the bounding box needs to be relative */
330  /* to a false origin of 0,0, so subtract the false easting */
331  /* and false northing values from the delta easting and */
332  /* delta northing values */
333  GeodeticCoordinates gcTemp( CoordinateType::geodetic, ( Merc_Cent_Mer + PI ), MAX_LAT );
334  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
335  Merc_Delta_Easting = tempCoordinates->easting();
336  Merc_Delta_Northing = tempCoordinates->northing();
337  delete tempCoordinates;
338 
339  if(Merc_False_Easting)
340  Merc_Delta_Easting -= Merc_False_Easting;
341  if (Merc_Delta_Easting < 0)
342  Merc_Delta_Easting = -Merc_Delta_Easting;
343  Merc_Delta_Easting *= 1.01;
344 
345  if(Merc_False_Northing)
346  Merc_Delta_Northing -= Merc_False_Northing;
347  if (Merc_Delta_Northing < 0)
348  Merc_Delta_Northing = -Merc_Delta_Northing;
349  Merc_Delta_Northing *= 1.01;
350 }
351 
352 
353 Mercator::Mercator( const Mercator &m )
354 {
355  coordinateType = m.coordinateType;
356  semiMajorAxis = m.semiMajorAxis;
357  flattening = m.flattening;
358  Merc_e = m.Merc_e;
359  Merc_es = m.Merc_es;
360  Merc_Cent_Mer = m.Merc_Cent_Mer;
361  Merc_Standard_Parallel = m.Merc_Standard_Parallel;
362  Merc_False_Easting = m.Merc_False_Easting;
363  Merc_False_Northing = m.Merc_False_Northing;
364  Merc_Scale_Factor = m.Merc_Scale_Factor;
365  Merc_ab = m.Merc_ab;
366  Merc_bb = m.Merc_bb;
367  Merc_cb = m.Merc_cb;
368  Merc_db = m.Merc_db;
369  Merc_Delta_Easting = m.Merc_Delta_Easting;
370  Merc_Delta_Northing = m.Merc_Delta_Northing;
371 }
372 
373 
374 Mercator::~Mercator()
375 {
376 }
377 
378 
379 Mercator& Mercator::operator=( const Mercator &m )
380 {
381  if( this != &m )
382  {
383  coordinateType = m.coordinateType;
384  semiMajorAxis = m.semiMajorAxis;
385  flattening = m.flattening;
386  Merc_e = m.Merc_e;
387  Merc_es = m.Merc_es;
388  Merc_Cent_Mer = m.Merc_Cent_Mer;
389  Merc_Standard_Parallel = m.Merc_Standard_Parallel;
390  Merc_False_Easting = m.Merc_False_Easting;
391  Merc_False_Northing = m.Merc_False_Northing;
392  Merc_Scale_Factor = m.Merc_Scale_Factor;
393  Merc_ab = m.Merc_ab;
394  Merc_bb = m.Merc_bb;
395  Merc_cb = m.Merc_cb;
396  Merc_db = m.Merc_db;
397  Merc_Delta_Easting = m.Merc_Delta_Easting;
398  Merc_Delta_Northing = m.Merc_Delta_Northing;
399  }
400 
401  return *this;
402 }
403 
404 
405 MercatorStandardParallelParameters* Mercator::getStandardParallelParameters() const
406 {
407 /*
408  * The function getStandardParallelParameters returns the current ellipsoid
409  * parameters and Mercator (Standard Parallel) projection parameters.
410  *
411  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
412  * ellipsoidFlattening : Flattening of ellipsoid (output)
413  * centralMeridian : Longitude in radians at the center of (output)
414  * the projection
415  * standardParallel : Latitude in radians at which the (output)
416  * point scale factor is 1.0
417  * falseEasting : A coordinate value in meters assigned to the
418  * central meridian of the projection. (output)
419  * falseNorthing : A coordinate value in meters assigned to the
420  * origin latitude of the projection (output)
421  */
422 
423  return new MercatorStandardParallelParameters(CoordinateType::mercatorStandardParallel, Merc_Cent_Mer, Merc_Standard_Parallel, Merc_Scale_Factor, Merc_False_Easting, Merc_False_Northing);
424 }
425 
426 
427 MercatorScaleFactorParameters* Mercator::getScaleFactorParameters() const
428 {
429 /*
430  * The function getScaleFactorParameters returns the current ellipsoid
431  * parameters and Mercator (Scale Factor) projection parameters.
432  *
433  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
434  * ellipsoidFlattening : Flattening of ellipsoid (output)
435  * centralMeridian : Longitude in radians at the center of (output)
436  * the projection
437  * falseEasting : A coordinate value in meters assigned to the
438  * central meridian of the projection. (output)
439  * falseNorthing : A coordinate value in meters assigned to the
440  * origin latitude of the projection (output)
441  * scaleFactor : Multiplier which reduces distances in the
442  * projection to the actual distance on the
443  * ellipsoid (output)
444  */
445 
446  return new MercatorScaleFactorParameters( CoordinateType::mercatorScaleFactor, Merc_Cent_Mer, Merc_Scale_Factor, Merc_False_Easting, Merc_False_Northing );
447 }
448 
449 
450 MSP::CCS::MapProjectionCoordinates* Mercator::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
451 {
452 /*
453  * The function convertFromGeodetic converts geodetic (latitude and
454  * longitude) coordinates to Mercator projection (easting and northing)
455  * coordinates, according to the current ellipsoid and Mercator projection
456  * parameters. If any errors occur, an exception is thrown with a description
457  * of the error.
458  *
459  * longitude : Longitude (lambda) in radians (input)
460  * latitude : Latitude (phi) in radians (input)
461  * easting : Easting (X) in meters (output)
462  * northing : Northing (Y) in meters (output)
463  */
464 
465  double ctanz2; /* Cotangent of z/2 - z - Isometric colatitude */
466  double e_x_sinlat; /* e * sin(latitude) */
467  double Delta_Long; /* Difference in origin longitude and longitude */
468  double tan_temp;
469  double pow_temp;
470 
471  double longitude = geodeticCoordinates->longitude();
472  double latitude = geodeticCoordinates->latitude();
473 
474  if ((latitude < -MAX_LAT) || (latitude > MAX_LAT))
475  { /* Latitude out of range */
476  throw CoordinateConversionException( ErrorMessages::latitude );
477  }
478  if ((longitude < -PI) || (longitude > TWO_PI))
479  { /* Longitude out of range */
480  throw CoordinateConversionException( ErrorMessages::longitude );
481  }
482 
483  if (longitude > PI)
484  longitude -= TWO_PI;
485  e_x_sinlat = Merc_e * sin(latitude);
486  tan_temp = tan(PI / 4.e0 + latitude / 2.e0);
487  pow_temp = pow( ((1.e0 - e_x_sinlat) / (1.e0 + e_x_sinlat)),
488  (Merc_e / 2.e0) );
489  ctanz2 = tan_temp * pow_temp;
490  double northing = Merc_Scale_Factor * semiMajorAxis * log(ctanz2) + Merc_False_Northing;
491  Delta_Long = longitude - Merc_Cent_Mer;
492  if (Delta_Long > PI)
493  Delta_Long -= TWO_PI;
494  if (Delta_Long < -PI)
495  Delta_Long += TWO_PI;
496  double easting = Merc_Scale_Factor * semiMajorAxis * Delta_Long
497  + Merc_False_Easting;
498 
499  return new MapProjectionCoordinates( coordinateType, easting, northing );
500 }
501 
502 
503 MSP::CCS::GeodeticCoordinates* Mercator::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
504 {
505 /*
506  * The function convertToGeodetic converts Mercator projection
507  * (easting and northing) coordinates to geodetic (latitude and longitude)
508  * coordinates, according to the current ellipsoid and Mercator projection
509  * coordinates. If any errors occur, an exception is thrown with a description
510  * of the error.
511  *
512  * easting : Easting (X) in meters (input)
513  * northing : Northing (Y) in meters (input)
514  * longitude : Longitude (lambda) in radians (output)
515  * latitude : Latitude (phi) in radians (output)
516  */
517 
518  double dx; /* Delta easting - Difference in easting (easting-FE) */
519  double dy; /* Delta northing - Difference in northing (northing-FN) */
520  double xphi; /* Isometric latitude */
521 
522  double easting = mapProjectionCoordinates->easting();
523  double northing = mapProjectionCoordinates->northing();
524 
525  if ((easting < (Merc_False_Easting - Merc_Delta_Easting))
526  || (easting > (Merc_False_Easting + Merc_Delta_Easting)))
527  { /* Easting out of range */
528  throw CoordinateConversionException( ErrorMessages::easting );
529  }
530  if ((northing < (Merc_False_Northing - Merc_Delta_Northing))
531  || (northing > (Merc_False_Northing + Merc_Delta_Northing)))
532  { /* Northing out of range */
533  throw CoordinateConversionException( ErrorMessages::northing );
534  }
535 
536  dy = northing - Merc_False_Northing;
537  dx = easting - Merc_False_Easting;
538  double longitude = Merc_Cent_Mer + dx / (Merc_Scale_Factor * semiMajorAxis);
539  xphi = PI_OVER_2
540  - 2.e0 * atan(1.e0 / exp(dy / (Merc_Scale_Factor * semiMajorAxis)));
541  double latitude = xphi + Merc_ab * sin(2.e0 * xphi) + Merc_bb * sin(4.e0 * xphi)
542  + Merc_cb * sin(6.e0 * xphi) + Merc_db * sin(8.e0 * xphi);
543 
544  if (longitude > PI)
545  longitude -= TWO_PI;
546  if (longitude < -PI)
547  longitude += TWO_PI;
548 
549  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
550 }
551 
552 
553 
554 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2