UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
TransverseCylindricalEqualArea.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: TRANSVERSE CYLINDRICAL EQUAL AREA
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Transverse Cylindrical Equal Area
10  * projection coordinates (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  * TCEA_NO_ERROR : No errors occurred in function
20  * TCEA_LAT_ERROR : Latitude outside of valid range
21  * (-90 to 90 degrees)
22  * TCEA_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * TCEA_EASTING_ERROR : Easting outside of valid range
25  * (False_Easting +/- ~6,500,000 m,
26  * depending on ellipsoid parameters
27  * and Origin_Latitude)
28  * TCEA_NORTHING_ERROR : Northing outside of valid range
29  * (False_Northing +/- ~20,000,000 m,
30  * depending on ellipsoid parameters
31  * and Origin_Latitude)
32  * TCEA_ORIGIN_LAT_ERROR : Origin latitude outside of valid range
33  * (-90 to 90 degrees)
34  * TCEA_CENT_MER_ERROR : Central meridian outside of valid range
35  * (-180 to 360 degrees)
36  * TCEA_A_ERROR : Semi-major axis less than or equal to zero
37  * TCEA_INV_F_ERROR : Inverse flattening outside of valid range
38  * (250 to 350)
39  * TCEA_SCALE_FACTOR_ERROR : Scale factor outside of valid
40  * range (0.3 to 3.0)
41  * TCEA_LON_WARNING : Distortion will result if longitude is more
42  * than 90 degrees from the Central Meridian
43  *
44  * REUSE NOTES
45  *
46  * TRANSVERSE CYLINDRICAL EQUAL AREA is intended for reuse by any application that
47  * performs a Transverse Cylindrical Equal Area projection or its inverse.
48  *
49  * REFERENCES
50  *
51  * Further information on TRANSVERSE CYLINDRICAL EQUAL AREA can be found in the Reuse Manual.
52  *
53  * TRANSVERSE CYLINDRICAL EQUAL AREA originated from :
54  * U.S. Army Topographic Engineering Center
55  * Geospatial Information Division
56  * 7701 Telegraph Road
57  * Alexandria, VA 22310-3864
58  *
59  * LICENSES
60  *
61  * None apply to this component.
62  *
63  * RESTRICTIONS
64  *
65  * TRANSVERSE CYLINDRICAL EQUAL AREA has no restrictions.
66  *
67  * ENVIRONMENT
68  *
69  * TRANSVERSE CYLINDRICAL EQUAL AREA was tested and certified in the following environments:
70  *
71  * 1. Solaris 2.5 with GCC, version 2.8.1
72  * 2. Windows 95 with MS Visual C++, version 6
73  *
74  * MODIFICATIONS
75  *
76  * Date Description
77  * ---- -----------
78  * 3-1-07 Original C++ Code
79  *
80  */
81 
82 
83 /***************************************************************************/
84 /*
85  * INCLUDES
86  */
87 
88 #include <math.h>
89 #include "TransverseCylindricalEqualArea.h"
90 #include "MapProjection5Parameters.h"
91 #include "MapProjectionCoordinates.h"
92 #include "GeodeticCoordinates.h"
93 #include "CoordinateConversionException.h"
94 #include "ErrorMessages.h"
95 #include "WarningMessages.h"
96 
97 /*
98  * math.h - Standard C math library
99  * TransverseCylindricalEqualArea.h - Is for prototype error checking
100  * MapProjectionCoordinates.h - defines map projection coordinates
101  * GeodeticCoordinates.h - defines geodetic coordinates
102  * CoordinateConversionException.h - Exception handler
103  * ErrorMessages.h - Contains exception messages
104  * WarningMessages.h - Contains warning messages
105  */
106 
107 
108 using namespace MSP::CCS;
109 
110 
111 /***************************************************************************/
112 /* DEFINES
113  *
114  */
115 
116 const double PI = 3.14159265358979323e0; /* PI */
117 const double PI_OVER_2 = ( PI / 2.0);
118 const double TWO_PI = (2.0 * PI);
119 const double MIN_SCALE_FACTOR = 0.3;
120 const double MAX_SCALE_FACTOR = 3.0;
121 
122 
123 /************************************************************************/
124 /* FUNCTIONS
125  *
126  */
127 
128 TransverseCylindricalEqualArea::TransverseCylindricalEqualArea(
129  double ellipsoidSemiMajorAxis,
130  double ellipsoidFlattening,
131  double centralMeridian,
132  double latitudeOfTrueScale,
133  double falseEasting,
134  double falseNorthing,
135  double scaleFactor ) :
136  CoordinateSystem(),
137  es2( 0.0066943799901413800 ),
138  es4( 4.4814723452405e-005 ),
139  es6( 3.0000678794350e-007 ),
140  es( .081819190842622 ),
141  M0( 0.0 ),
142  qp( 1.9955310875028 ),
143  One_MINUS_es2( .99330562000986 ),
144  One_OVER_2es( 6.1110357466348 ),
145  a0( .0022392088624809 ),
146  a1( 2.8830839728915e-006 ),
147  a2( 5.0331826636906e-009 ),
148  b0( .0025188265843907 ),
149  b1( 3.7009490356205e-006 ),
150  b2( 7.4478137675038e-009 ),
151  b3( 1.7035993238596e-011 ),
152  c0( .99832429845280 ),
153  c1( .0025146070605187 ),
154  c2( 2.6390465943377e-006 ),
155  c3( 3.4180460865959e-009 ),
156  Tcea_Origin_Long( 0.0 ),
157  Tcea_Origin_Lat( 0.0 ),
158  Tcea_False_Easting( 0.0 ),
159  Tcea_False_Northing( 0.0 ),
160  Tcea_Scale_Factor( 1.0 ),
161  Tcea_Min_Easting( -6398628.0 ),
162  Tcea_Max_Easting( 6398628.0 ),
163  Tcea_Min_Northing( -20003931.0 ),
164  Tcea_Max_Northing( 20003931.0 )
165 {
166 /*
167  * The constructor receives the ellipsoid parameters and
168  * Transverse Cylindrical Equal Area projection parameters as inputs,
169  * and sets the corresponding state variables.
170  * If any errors occur, an exception is thrown with a description of the error.
171  *
172  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
173  * ellipsoidFlattening : Flattening of ellipsoid (input)
174  * centralMeridian : Longitude in radians at the center of (input)
175  * the projection
176  * latitudeOfTrueScale : Latitude in radians at which the (input)
177  * point scale factor is 1.0
178  * falseEasting : A coordinate value in meters assigned to the
179  * central meridian of the projection. (input)
180  * falseNorthing : A coordinate value in meters assigned to the
181  * origin latitude of the projection (input)
182  * scaleFactor : Multiplier which reduces distances in the
183  * projection to the actual distance on the
184  * ellipsoid (input)
185  */
186 
187  double sin_lat_90 = sin(PI_OVER_2);
188  double x, j, three_es4;
189  double Sqrt_One_MINUS_es2;
190  double e1, e2, e3, e4;
191  double lat, sin2lat, sin4lat, sin6lat;
192  double temp, temp_northing;
193  double inv_f = 1 / ellipsoidFlattening;
194 
195  if (ellipsoidSemiMajorAxis <= 0.0)
196  { /* Semi-major axis must be greater than zero */
197  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
198  }
199  if ((inv_f < 250) || (inv_f > 350))
200  { /* Inverse flattening must be between 250 and 350 */
201  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
202  }
203  if ((latitudeOfTrueScale < -PI_OVER_2) || (latitudeOfTrueScale > PI_OVER_2))
204  { /* origin latitude out of range */
205  throw CoordinateConversionException( ErrorMessages::originLatitude );
206  }
207  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
208  { /* origin longitude out of range */
209  throw CoordinateConversionException( ErrorMessages::centralMeridian );
210  }
211  if ((scaleFactor < MIN_SCALE_FACTOR) || (scaleFactor > MAX_SCALE_FACTOR))
212  {
213  throw CoordinateConversionException( ErrorMessages::scaleFactor );
214  }
215 
216  semiMajorAxis = ellipsoidSemiMajorAxis;
217  flattening = ellipsoidFlattening;
218 
219  Tcea_Origin_Lat = latitudeOfTrueScale;
220  if (centralMeridian > PI)
221  centralMeridian -= TWO_PI;
222  Tcea_Origin_Long = centralMeridian;
223  Tcea_False_Northing = falseNorthing;
224  Tcea_False_Easting = falseEasting;
225  Tcea_Scale_Factor = scaleFactor;
226 
227  es2 = 2 * flattening - flattening * flattening;
228  es = sqrt(es2);
229  One_MINUS_es2 = 1.0 - es2;
230  Sqrt_One_MINUS_es2 = sqrt(One_MINUS_es2);
231  One_OVER_2es = 1.0 / (2.0 * es);
232  es4 = es2 * es2;
233  es6 = es4 * es2;
234  x = es * sin_lat_90;
235  qp = TCEA_Q(sin_lat_90,x);
236 
237  a0 = es2 / 3.0 + 31.0 * es4 / 180.0 + 517.0 * es6 / 5040.0;
238  a1 = 23.0 * es4 / 360.0 + 251.0 * es6 / 3780.0;
239  a2 = 761.0 * es6 / 45360.0;
240 
241  e1 = (1.0 - Sqrt_One_MINUS_es2) / (1.0 + Sqrt_One_MINUS_es2);
242  e2 = e1 * e1;
243  e3 = e2 * e1;
244  e4 = e3 * e1;
245  b0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0;
246  b1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
247  b2 = 151.0 * e3 / 96.0;
248  b3 = 1097.0 * e4 / 512.0;
249 
250  j = 45.0 * es6 / 1024.0;
251  three_es4 = 3.0 * es4;
252  c0 = 1.0 - es2 / 4.0 - three_es4 / 64.0 - 5.0 * es6 / 256.0;
253  c1 = 3.0 * es2 / 8.0 + three_es4 / 32.0 + j;
254  c2 = 15.0 * es4 / 256.0 + j;
255  c3 = 35.0 * es6 / 3072.0;
256  lat = c0 * Tcea_Origin_Lat;
257  sin2lat = TCEA_COEFF_TIMES_SIN(c1, 2.0, Tcea_Origin_Lat);
258  sin4lat = TCEA_COEFF_TIMES_SIN(c2, 4.0, Tcea_Origin_Lat);
259  sin6lat = TCEA_COEFF_TIMES_SIN(c3, 6.0, Tcea_Origin_Lat);
260  M0 = TCEA_M(lat, sin2lat, sin4lat, sin6lat);
261 
262  GeodeticCoordinates gcTemp( CoordinateType::geodetic, PI, PI_OVER_2 );
263  MapProjectionCoordinates* tempCoordinates = convertFromGeodetic( &gcTemp );
264  temp = tempCoordinates->easting();
265  temp_northing = tempCoordinates->northing();
266  delete tempCoordinates;
267 
268  if (temp_northing > 0)
269  {
270  Tcea_Min_Northing = temp_northing - 20003931.458986;
271  Tcea_Max_Northing = temp_northing;
272 
273  if(Tcea_False_Northing)
274  {
275  Tcea_Min_Northing -= Tcea_False_Northing;
276  Tcea_Max_Northing -= Tcea_False_Northing;
277  }
278  }
279  else if (temp_northing < 0)
280  {
281  Tcea_Max_Northing = temp_northing + 20003931.458986;
282  Tcea_Min_Northing = temp_northing;
283 
284  if(Tcea_False_Northing)
285  {
286  Tcea_Min_Northing -= Tcea_False_Northing;
287  Tcea_Max_Northing -= Tcea_False_Northing;
288  }
289  }
290 }
291 
292 
293 TransverseCylindricalEqualArea::TransverseCylindricalEqualArea(
294  const TransverseCylindricalEqualArea &tcea )
295 {
296  semiMajorAxis = tcea.semiMajorAxis;
297  flattening = tcea.flattening;
298  es2 = tcea.es2;
299  es4 = tcea.es4;
300  es6 = tcea.es6;
301  es = tcea.es;
302  M0 = tcea.M0;
303  qp = tcea.qp;
304  One_MINUS_es2 = tcea.One_MINUS_es2;
305  One_OVER_2es = tcea.One_OVER_2es;
306  a0 = tcea.a0;
307  a1 = tcea.a1;
308  a2 = tcea.a2;
309  b0 = tcea.b0;
310  b1 = tcea.b1;
311  b2 = tcea.b2;
312  b3 = tcea.b3;
313  c0 = tcea.c0;
314  c1 = tcea.c1;
315  c2 = tcea.c2;
316  c3 = tcea.c3;
317  Tcea_Origin_Long = tcea.Tcea_Origin_Long;
318  Tcea_Origin_Lat = tcea.Tcea_Origin_Lat;
319  Tcea_False_Easting = tcea.Tcea_False_Easting;
320  Tcea_False_Northing = tcea.Tcea_False_Northing;
321  Tcea_Scale_Factor = tcea.Tcea_Scale_Factor;
322  Tcea_Min_Easting = tcea.Tcea_Min_Easting;
323  Tcea_Max_Easting = tcea.Tcea_Max_Easting;
324  Tcea_Min_Northing = tcea.Tcea_Min_Northing;
325  Tcea_Max_Northing = tcea.Tcea_Max_Northing;
326 }
327 
328 
329 TransverseCylindricalEqualArea::~TransverseCylindricalEqualArea()
330 {
331 }
332 
333 
334 TransverseCylindricalEqualArea& TransverseCylindricalEqualArea::operator=( const TransverseCylindricalEqualArea &tcea )
335 {
336  if( this != &tcea )
337  {
338  semiMajorAxis = tcea.semiMajorAxis;
339  flattening = tcea.flattening;
340  es2 = tcea.es2;
341  es4 = tcea.es4;
342  es6 = tcea.es6;
343  es = tcea.es;
344  M0 = tcea.M0;
345  qp = tcea.qp;
346  One_MINUS_es2 = tcea.One_MINUS_es2;
347  One_OVER_2es = tcea.One_OVER_2es;
348  a0 = tcea.a0;
349  a1 = tcea.a1;
350  a2 = tcea.a2;
351  b0 = tcea.b0;
352  b1 = tcea.b1;
353  b2 = tcea.b2;
354  b3 = tcea.b3;
355  c0 = tcea.c0;
356  c1 = tcea.c1;
357  c2 = tcea.c2;
358  c3 = tcea.c3;
359  Tcea_Origin_Long = tcea.Tcea_Origin_Long;
360  Tcea_Origin_Lat = tcea.Tcea_Origin_Lat;
361  Tcea_False_Easting = tcea.Tcea_False_Easting;
362  Tcea_False_Northing = tcea.Tcea_False_Northing;
363  Tcea_Scale_Factor = tcea.Tcea_Scale_Factor;
364  Tcea_Min_Easting = tcea.Tcea_Min_Easting;
365  Tcea_Max_Easting = tcea.Tcea_Max_Easting;
366  Tcea_Min_Northing = tcea.Tcea_Min_Northing;
367  Tcea_Max_Northing = tcea.Tcea_Max_Northing;
368  }
369 
370  return *this;
371 }
372 
373 
374 MapProjection5Parameters* TransverseCylindricalEqualArea::getParameters() const
375 {
376 /*
377  * The function getParameters returns the current ellipsoid
378  * parameters, Transverse Cylindrical Equal Area projection parameters, and scale factor.
379  *
380  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
381  * ellipsoidFlattening : Flattening of ellipsoid (output)
382  * centralMeridian : Longitude in radians at the center of (output)
383  * the projection
384  * latitudeOfTrueScale : Latitude in radians at which the (output)
385  * point scale factor is 1.0
386  * falseEasting : A coordinate value in meters assigned to the
387  * central meridian of the projection. (output)
388  * falseNorthing : A coordinate value in meters assigned to the
389  * origin latitude of the projection (output)
390  * scaleFactor : Multiplier which reduces distances in the
391  * projection to the actual distance on the
392  * ellipsoid (output)
393  */
394 
395  return new MapProjection5Parameters( CoordinateType::transverseCylindricalEqualArea, Tcea_Origin_Long, Tcea_Origin_Lat, Tcea_Scale_Factor, Tcea_False_Easting, Tcea_False_Northing );
396 }
397 
398 
399 MSP::CCS::MapProjectionCoordinates* TransverseCylindricalEqualArea::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
400 {
401 /*
402  * The function convertFromGeodetic converts geodetic (latitude and
403  * longitude) coordinates to Transverse Cylindrical Equal Area projection (easting and
404  * northing) coordinates, according to the current ellipsoid and Transverse Cylindrical
405  * Equal Area projection parameters. If any errors occur, an exception is thrown with a description
406  * of the error.
407  *
408  * longitude : Longitude (lambda) in radians (input)
409  * latitude : Latitude (phi) in radians (input)
410  * easting : Easting (X) in meters (output)
411  * northing : Northing (Y) in meters (output)
412  */
413 
414  double x;
415  double dlam; /* Longitude - Central Meridan */
416  double qq, qq_OVER_qp;
417  double beta, betac;
418  double sin2betac, sin4betac, sin6betac;
419  double PHIc;
420  double phi, sin2phi, sin4phi, sin6phi;
421  double sinPHIc;
422  double Mc;
423 
424  double longitude = geodeticCoordinates->longitude();
425  double latitude = geodeticCoordinates->latitude();
426  double sin_lat = sin(latitude);
427 
428  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
429  { /* Latitude out of range */
430  throw CoordinateConversionException( ErrorMessages::latitude );
431  }
432  if ((longitude < -PI) || (longitude > TWO_PI))
433  { /* Longitude out of range */
434  throw CoordinateConversionException( ErrorMessages::longitude );
435  }
436 
437  dlam = longitude - Tcea_Origin_Long;
438 
439  char warning[100];
440  warning[0] = '\0';
441  if (fabs(dlam) >= (PI / 2))
442  { /* Distortion resultsif Longitude is > 90 deg from the Central Meridian */
443  strcat( warning, MSP::CCS::WarningMessages::longitude );
444  }
445 
446  if (dlam > PI)
447  {
448  dlam -= TWO_PI;
449  }
450  if (dlam < -PI)
451  {
452  dlam += TWO_PI;
453  }
454  if (latitude == PI_OVER_2)
455  {
456  qq = qp;
457  qq_OVER_qp = 1.0;
458  }
459  else
460  {
461  x = es * sin_lat;
462  qq = TCEA_Q(sin_lat, x);
463  qq_OVER_qp = qq / qp;
464  }
465 
466 
467  if (qq_OVER_qp > 1.0)
468  qq_OVER_qp = 1.0;
469  else if (qq_OVER_qp < -1.0)
470  qq_OVER_qp = -1.0;
471 
472  beta = asin(qq_OVER_qp);
473  betac = atan(tan(beta) / cos(dlam));
474 
475  if ((fabs(betac) - PI_OVER_2) > 1.0e-8)
476  PHIc = betac;
477  else
478  {
479  sin2betac = TCEA_COEFF_TIMES_SIN(a0, 2.0, betac);
480  sin4betac = TCEA_COEFF_TIMES_SIN(a1, 4.0, betac);
481  sin6betac = TCEA_COEFF_TIMES_SIN(a2, 6.0, betac);
482  PHIc = TCEA_L(betac, sin2betac, sin4betac, sin6betac);
483  }
484 
485  sinPHIc = sin(PHIc);
486  double easting = semiMajorAxis * cos(beta) * cos(PHIc) * sin(dlam) /
487  (Tcea_Scale_Factor * cos(betac) * sqrt(1.0 - es2 *
488  sinPHIc * sinPHIc)) + Tcea_False_Easting;
489 
490  phi = c0 * PHIc;
491  sin2phi = TCEA_COEFF_TIMES_SIN(c1, 2.0, PHIc);
492  sin4phi = TCEA_COEFF_TIMES_SIN(c2, 4.0, PHIc);
493  sin6phi = TCEA_COEFF_TIMES_SIN(c3, 6.0, PHIc);
494  Mc = TCEA_M(phi, sin2phi, sin4phi, sin6phi);
495 
496  double northing = Tcea_Scale_Factor * (Mc - M0) + Tcea_False_Northing;
497 
498  return new MapProjectionCoordinates(
499  CoordinateType::transverseCylindricalEqualArea,
500  warning, easting, northing );
501 }
502 
503 
504 MSP::CCS::GeodeticCoordinates* TransverseCylindricalEqualArea::convertToGeodetic(
505  MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
506 {
507 /*
508  * The function convertToGeodetic converts Transverse
509  * Cylindrical Equal Area projection (easting and northing) coordinates
510  * to geodetic (latitude and longitude) coordinates, according to the
511  * current ellipsoid and Transverse Cylindrical Equal Area projection
512  * coordinates. If any errors occur, an exception is thrown with a description
513  * of the error.
514  *
515  * easting : Easting (X) in meters (input)
516  * northing : Northing (Y) in meters (input)
517  * longitude : Longitude (lambda) in radians (output)
518  * latitude : Latitude (phi) in radians (output)
519  */
520 
521  double x;
522  double dx; /* Delta easting - Difference in easting (easting-FE) */
523  double dy; /* Delta northing - Difference in northing (northing-FN) */
524  double Mc;
525  double MUc;
526  double sin2mu, sin4mu, sin6mu, sin8mu;
527  double PHIc;
528  double Qc;
529  double sin_lat;
530  double beta, betac, beta_prime;
531  double sin2beta, sin4beta, sin6beta;
532  double cosbetac;
533  double Qc_OVER_qp;
534  double temp;
535 
536  double easting = mapProjectionCoordinates->easting();
537  double northing = mapProjectionCoordinates->northing();
538 
539  if ((easting < (Tcea_False_Easting + Tcea_Min_Easting))
540  || (easting > (Tcea_False_Easting + Tcea_Max_Easting)))
541  { /* Easting out of range */
542  throw CoordinateConversionException( ErrorMessages::easting );
543  }
544  if( (northing < (Tcea_False_Northing + Tcea_Min_Northing))
545  || (northing > (Tcea_False_Northing + Tcea_Max_Northing)))
546  { /* Northing out of range */
547  throw CoordinateConversionException( ErrorMessages::northing );
548  }
549 
550  dy = northing - Tcea_False_Northing;
551  dx = easting - Tcea_False_Easting;
552  Mc = M0 + dy / Tcea_Scale_Factor;
553  MUc = Mc / (semiMajorAxis * c0);
554 
555  sin2mu = TCEA_COEFF_TIMES_SIN(b0, 2.0, MUc);
556  sin4mu = TCEA_COEFF_TIMES_SIN(b1, 4.0, MUc);
557  sin6mu = TCEA_COEFF_TIMES_SIN(b2, 6.0, MUc);
558  sin8mu = TCEA_COEFF_TIMES_SIN(b3, 8.0, MUc);
559  PHIc = MUc + sin2mu + sin4mu + sin6mu + sin8mu;
560 
561  sin_lat = sin(PHIc);
562  x = es * sin_lat;
563  Qc = TCEA_Q(sin_lat, x);
564  Qc_OVER_qp = Qc / qp;
565 
566  if (Qc_OVER_qp < -1.0)
567  Qc_OVER_qp = -1.0;
568  else if (Qc_OVER_qp > 1.0)
569  Qc_OVER_qp = 1.0;
570 
571  betac = asin(Qc_OVER_qp);
572  cosbetac = cos(betac);
573  temp = Tcea_Scale_Factor * dx * cosbetac * sqrt(1.0 -
574  es2 * sin_lat * sin_lat) / (semiMajorAxis * cos(PHIc));
575  if (temp > 1.0)
576  temp = 1.0;
577  else if (temp < -1.0)
578  temp = -1.0;
579  beta_prime = -asin(temp);
580  beta = asin(cos(beta_prime) * sin(betac));
581 
582  sin2beta = TCEA_COEFF_TIMES_SIN(a0, 2.0, beta);
583  sin4beta = TCEA_COEFF_TIMES_SIN(a1, 4.0, beta);
584  sin6beta = TCEA_COEFF_TIMES_SIN(a2, 6.0, beta);
585  double latitude = TCEA_L(beta, sin2beta, sin4beta, sin6beta);
586 
587  double longitude = Tcea_Origin_Long - atan(tan(beta_prime) / cosbetac);
588 
589  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
590  latitude = PI_OVER_2;
591  else if (latitude < -PI_OVER_2)
592  latitude = -PI_OVER_2;
593 
594  if (longitude > PI)
595  longitude -= TWO_PI;
596  if (longitude < -PI)
597  longitude += TWO_PI;
598 
599  if (longitude > PI) /* force distorted values to 180, -180 degrees */
600  longitude = PI;
601  else if (longitude < -PI)
602  longitude = -PI;
603 
604  return new GeodeticCoordinates(
605  CoordinateType::geodetic, longitude, latitude );
606 }
607 
608 
609 double TransverseCylindricalEqualArea::TCEA_Q( double sinlat, double x )
610 {
611  return One_MINUS_es2*(sinlat/(1.0-es2*sinlat*sinlat)-One_OVER_2es*log((1-x)/(1+x)));
612 }
613 
614 
615 double TransverseCylindricalEqualArea::TCEA_COEFF_TIMES_SIN(
616  double coeff, double x, double latit )
617 {
618  return coeff * sin(x*latit);
619 }
620 
621 
622 double TransverseCylindricalEqualArea::TCEA_M(
623  double c0lat, double c1lat, double c2lat, double c3lat )
624 {
625  return semiMajorAxis * (c0lat - c1lat + c2lat - c3lat);
626 }
627 
628 
629 double TransverseCylindricalEqualArea::TCEA_L(
630  double Beta, double c0lat, double c1lat, double c2lat )
631 {
632  return Beta + c0lat + c1lat + c2lat;
633 }
634 
635 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2