UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
Sinusoidal.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: SINUSOIDAL
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Sinusoid 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  * SINU_NO_ERROR : No errors occurred in function
20  * SINU_LAT_ERROR : Latitude outside of valid range
21  * (-90 to 90 degrees)
22  * SINU_LON_ERROR : Longitude outside of valid range
23  * (-180 to 360 degrees)
24  * SINU_EASTING_ERROR : Easting outside of valid range
25  * (False_Easting +/- ~20,000,000 m,
26  * depending on ellipsoid parameters)
27  * SINU_NORTHING_ERROR : Northing outside of valid range
28  * (False_Northing +/- ~10,000,000 m,
29  * depending on ellipsoid parameters)
30  * SINU_CENT_MER_ERROR : Origin longitude outside of valid range
31  * (-180 to 360 degrees)
32  * SINU_A_ERROR : Semi-major axis less than or equal to zero
33  * SINU_INV_F_ERROR : Inverse flattening outside of valid range
34  * (250 to 350)
35  *
36  * REUSE NOTES
37  *
38  * SINUSOIDAL is intended for reuse by any application that performs a
39  * Sinusoid projection or its inverse.
40  *
41  * REFERENCES
42  *
43  * Further information on SINUSOIDAL can be found in the Reuse Manual.
44  *
45  * SINUSOIDAL originated from : U.S. Army Topographic Engineering Center
46  * Geospatial Information Division
47  * 7701 Telegraph Road
48  * Alexandria, VA 22310-3864
49  *
50  * LICENSES
51  *
52  * None apply to this component.
53  *
54  * RESTRICTIONS
55  *
56  * SINUSOIDAL has no restrictions.
57  *
58  * ENVIRONMENT
59  *
60  * SINUSOIDAL was tested and certified in the following environments:
61  *
62  * 1. Solaris 2.5 with GCC, version 2.8.1
63  * 2. Windows 95 with MS Visual C++, version 6
64  *
65  * MODIFICATIONS
66  *
67  * Date Description
68  * ---- -----------
69  * 07-15-99 Original Code
70  * 03-05-07 Original C++ Code
71  *
72  */
73 
74 
75 /***************************************************************************/
76 /*
77  * INCLUDES
78  */
79 
80 #include <math.h>
81 #include "Sinusoidal.h"
82 #include "MapProjection3Parameters.h"
83 #include "MapProjectionCoordinates.h"
84 #include "GeodeticCoordinates.h"
85 #include "CoordinateConversionException.h"
86 #include "ErrorMessages.h"
87 
88 /*
89  * math.h - Standard C math library
90  * Sinusoidal.h - Is for prototype error checking
91  * MapProjectionCoordinates.h - defines map projection coordinates
92  * GeodeticCoordinates.h - defines geodetic coordinates
93  * CoordinateConversionException.h - Exception handler
94  * ErrorMessages.h - Contains exception messages
95  */
96 
97 
98 using namespace MSP::CCS;
99 
100 
101 /***************************************************************************/
102 /* DEFINES
103  *
104  */
105 
106 const double PI = 3.14159265358979323e0; /* PI */
107 const double PI_OVER_2 = (PI / 2.0);
108 const double TWO_PI = (2.0 * PI);
109 
110 
111 double sinuCoeffTimesSine( double coeff, double x, double latit )
112 {
113  return coeff * sin(x * latit);
114 }
115 
116 
117 double floatEq( double x, double v, double epsilon )
118 {
119  return ((v - epsilon) < x) && (x < (v + epsilon));
120 }
121 
122 
123 /************************************************************************/
124 /* FUNCTIONS
125  *
126  */
127 
128 Sinusoidal::Sinusoidal( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
129  CoordinateSystem(),
130  es2( 0.0066943799901413800 ),
131  es4( 4.4814723452405e-005 ),
132  es6( 3.0000678794350e-007 ),
133  c0( .99832429845280 ),
134  c1( .0025146070605187 ),
135  c2( 2.6390465943377e-006 ),
136  c3( 3.4180460865959e-009 ),
137  a0( .0025188265843907 ),
138  a1( 3.7009490356205e-006 ),
139  a2( 7.4478137675038e-009 ),
140  a3( 1.7035993238596e-011 ),
141  Sinu_Origin_Long( 0.0 ),
142  Sinu_False_Easting( 0.0 ),
143  Sinu_False_Northing( 0.0 ),
144  Sinu_Max_Easting( 20037509.0 ),
145  Sinu_Min_Easting( -20037509.0 ),
146  Sinu_Delta_Northing( 10001966.0 )
147 {
148 /*
149  * The constructor receives the ellipsoid parameters and
150  * Sinusoidal projection parameters as inputs, and sets the corresponding state
151  * variables. If any errors occur, an exception is thrown with a description
152  * of the error.
153  *
154  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
155  * ellipsoidFlattening : Flattening of ellipsoid (input)
156  * centralMeridian : Longitude in radians at the center of (input)
157  * the projection
158  * falseEasting : A coordinate value in meters assigned to the
159  * central meridian of the projection. (input)
160  * falseNorthing : A coordinate value in meters assigned to the
161  * origin latitude of the projection (input)
162  */
163 
164  double j;
165  double One_MINUS_es2, Sqrt_One_MINUS_es2, e1, e2, e3, e4;
166  double inv_f = 1 / ellipsoidFlattening;
167 
168  if (ellipsoidSemiMajorAxis <= 0.0)
169  { /* Semi-major axis must be greater than zero */
170  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
171  }
172  if ((inv_f < 250) || (inv_f > 350))
173  { /* Inverse flattening must be between 250 and 350 */
174  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
175  }
176  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
177  { /* origin longitude out of range */
178  throw CoordinateConversionException( ErrorMessages::centralMeridian );
179  }
180 
181  semiMajorAxis = ellipsoidSemiMajorAxis;
182  flattening = ellipsoidFlattening;
183 
184  es2 = 2 * flattening - flattening * flattening;
185  es4 = es2 * es2;
186  es6 = es4 * es2;
187  j = 45.0 * es6 / 1024.0;
188  c0 = 1.0 - es2 / 4.0 - 3.0 * es4 / 64.0 - 5.0 * es6 / 256.0;
189  c1 = 3.0 * es2 / 8.0 + 3.0 * es4 / 32.0 + j;
190  c2 = 15.0 * es4 / 256.0 + j;
191  c3 = 35.0 * es6 / 3072.0;
192  One_MINUS_es2 = 1.0 - es2;
193  Sqrt_One_MINUS_es2 = sqrt(One_MINUS_es2);
194  e1 = (1.0 - Sqrt_One_MINUS_es2) / (1.0 + Sqrt_One_MINUS_es2);
195  e2 = e1 * e1;
196  e3 = e2 * e1;
197  e4 = e3 * e1;
198  a0 = 3.0 * e1 / 2.0 - 27.0 * e3 / 32.0 ;
199  a1 = 21.0 * e2 / 16.0 - 55.0 * e4 / 32.0;
200  a2 = 151.0 * e3 / 96.0;
201  a3 = 1097.0 * e4 / 512.0;
202  if (centralMeridian > PI)
203  centralMeridian -= TWO_PI;
204  Sinu_Origin_Long = centralMeridian;
205  Sinu_False_Northing = falseNorthing;
206  Sinu_False_Easting = falseEasting;
207 
208  if (Sinu_Origin_Long > 0)
209  {
210  Sinu_Max_Easting = 19926189.0;
211  Sinu_Min_Easting = -20037509.0;
212  }
213  else if (Sinu_Origin_Long < 0)
214  {
215  Sinu_Max_Easting = 20037509.0;
216  Sinu_Min_Easting = -19926189.0;
217  }
218  else
219  {
220  Sinu_Max_Easting = 20037509.0;
221  Sinu_Min_Easting = -20037509.0;
222  }
223 }
224 
225 
226 Sinusoidal::Sinusoidal( const Sinusoidal &s )
227 {
228  semiMajorAxis = s.semiMajorAxis;
229  flattening = s.flattening;
230  es2 = s.es2;
231  es4 = s.es4;
232  es6 = s.es6;
233  c0 = s.c0;
234  c1 = s.c1;
235  c2 = s.c2;
236  c3 = s.c3;
237  a0 = s.a0;
238  a1 = s.a1;
239  a2 = s.a2;
240  a3 = s.a3;
241  Sinu_Origin_Long = s.Sinu_Origin_Long;
242  Sinu_False_Easting = s.Sinu_False_Easting;
243  Sinu_False_Northing = s.Sinu_False_Northing;
244  Sinu_Max_Easting = s.Sinu_Max_Easting;
245  Sinu_Min_Easting = s.Sinu_Min_Easting;
246  Sinu_Delta_Northing = s.Sinu_Delta_Northing;
247 }
248 
249 
250 Sinusoidal::~Sinusoidal()
251 {
252 }
253 
254 
255 Sinusoidal& Sinusoidal::operator=( const Sinusoidal &s )
256 {
257  if( this != &s )
258  {
259  semiMajorAxis = s.semiMajorAxis;
260  flattening = s.flattening;
261  es2 = s.es2;
262  es4 = s.es4;
263  es6 = s.es6;
264  c0 = s.c0;
265  c1 = s.c1;
266  c2 = s.c2;
267  c3 = s.c3;
268  a0 = s.a0;
269  a1 = s.a1;
270  a2 = s.a2;
271  a3 = s.a3;
272  Sinu_Origin_Long = s.Sinu_Origin_Long;
273  Sinu_False_Easting = s.Sinu_False_Easting;
274  Sinu_False_Northing = s.Sinu_False_Northing;
275  Sinu_Max_Easting = s.Sinu_Max_Easting;
276  Sinu_Min_Easting = s.Sinu_Min_Easting;
277  Sinu_Delta_Northing = s.Sinu_Delta_Northing;
278  }
279 
280  return *this;
281 }
282 
283 
284 MapProjection3Parameters* Sinusoidal::getParameters() const
285 {
286 /*
287  * The function getParameters returns the current ellipsoid
288  * parameters, and Sinusoidal projection parameters.
289  *
290  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
291  * ellipsoidFlattening : Flattening of ellipsoid (output)
292  * centralMeridian : Longitude in radians at the center of (output)
293  * the projection
294  * falseEasting : A coordinate value in meters assigned to the
295  * central meridian of the projection. (output)
296  * falseNorthing : A coordinate value in meters assigned to the
297  * origin latitude of the projection (output)
298  */
299 
300  return new MapProjection3Parameters(
301  CoordinateType::sinusoidal,
302  Sinu_Origin_Long,
303  Sinu_False_Easting,
304  Sinu_False_Northing );
305 }
306 
307 
308 MSP::CCS::MapProjectionCoordinates* Sinusoidal::convertFromGeodetic(
309  MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
310 {
311 /*
312  * The function convertFromGeodetic converts geodetic (latitude and
313  * longitude) coordinates to Sinusoidal projection (easting and northing)
314  * coordinates, according to the current ellipsoid and Sinusoidal projection
315  * parameters. If any errors occur, an exception is thrown with a description
316  * of the error.
317  *
318  * longitude : Longitude (lambda) in radians (input)
319  * latitude : Latitude (phi) in radians (input)
320  * easting : Easting (X) in meters (output)
321  * northing : Northing (Y) in meters (output)
322  */
323 
324  double sin2lat, sin4lat, sin6lat;
325  double dlam; /* Longitude - Central Meridan */
326  double mm;
327  double MM;
328 
329  double longitude = geodeticCoordinates->longitude();
330  double latitude = geodeticCoordinates->latitude();
331  double slat = sin(latitude);
332 
333  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
334  { /* Latitude out of range */
335  throw CoordinateConversionException( ErrorMessages::latitude );
336  }
337  if ((longitude < -PI) || (longitude > TWO_PI))
338  { /* Longitude out of range */
339  throw CoordinateConversionException( ErrorMessages::longitude );
340  }
341 
342  dlam = longitude - Sinu_Origin_Long;
343  if (dlam > PI)
344  {
345  dlam -= TWO_PI;
346  }
347  if (dlam < -PI)
348  {
349  dlam += TWO_PI;
350  }
351  mm = sqrt(1.0 - es2 * slat * slat);
352 
353  sin2lat = sinuCoeffTimesSine(c1, 2.0, latitude);
354  sin4lat = sinuCoeffTimesSine(c2, 4.0, latitude);
355  sin6lat = sinuCoeffTimesSine(c3, 6.0, latitude);
356  MM = semiMajorAxis * (c0 * latitude - sin2lat + sin4lat - sin6lat);
357 
358  double easting = semiMajorAxis * dlam * cos(latitude) / mm + Sinu_False_Easting;
359  double northing = MM + Sinu_False_Northing;
360 
361  return new MapProjectionCoordinates( CoordinateType::sinusoidal, easting, northing );
362 }
363 
364 
365 MSP::CCS::GeodeticCoordinates* Sinusoidal::convertToGeodetic( MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
366 {
367 /*
368  * The function convertToGeodetic converts Sinusoidal projection
369  * (easting and northing) coordinates to geodetic (latitude and longitude)
370  * coordinates, according to the current ellipsoid and Sinusoidal projection
371  * coordinates. If any errors occur, an exception is thrown with a description
372  * of the error.
373  *
374  * easting : Easting (X) in meters (input)
375  * northing : Northing (Y) in meters (input)
376  * longitude : Longitude (lambda) in radians (output)
377  * latitude : Latitude (phi) in radians (output)
378  */
379 
380  double dx; /* Delta easting - Difference in easting (easting-FE) */
381  double dy; /* Delta northing - Difference in northing (northing-FN) */
382  double mu;
383  double sin2mu, sin4mu, sin6mu, sin8mu;
384  double sin_lat;
385  double longitude, latitude;
386 
387  double easting = mapProjectionCoordinates->easting();
388  double northing = mapProjectionCoordinates->northing();
389 
390  if ((easting < (Sinu_False_Easting + Sinu_Min_Easting))
391  || (easting > (Sinu_False_Easting + Sinu_Max_Easting)))
392  { /* Easting out of range */
393  throw CoordinateConversionException( ErrorMessages::easting );
394  }
395  if ((northing < (Sinu_False_Northing - Sinu_Delta_Northing))
396  || (northing > (Sinu_False_Northing + Sinu_Delta_Northing)))
397  { /* Northing out of range */
398  throw CoordinateConversionException( ErrorMessages::northing );
399  }
400 
401  dy = northing - Sinu_False_Northing;
402  dx = easting - Sinu_False_Easting;
403 
404  mu = dy / (semiMajorAxis * c0);
405  sin2mu = sinuCoeffTimesSine(a0, 2.0, mu);
406  sin4mu = sinuCoeffTimesSine(a1, 4.0, mu);
407  sin6mu = sinuCoeffTimesSine(a2, 6.0, mu);
408  sin8mu = sinuCoeffTimesSine(a3, 8.0, mu);
409  latitude = mu + sin2mu + sin4mu + sin6mu + sin8mu;
410 
411  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
412  latitude = PI_OVER_2;
413  else if (latitude < -PI_OVER_2)
414  latitude = -PI_OVER_2;
415 
416  if (floatEq(fabs(latitude),PI_OVER_2,1.0e-8))
417  longitude = Sinu_Origin_Long;
418  else
419  {
420  sin_lat = sin(latitude);
421  longitude = Sinu_Origin_Long + dx * sqrt(1.0 - es2 *
422  sin_lat * sin_lat) / (semiMajorAxis * cos(latitude));
423 
424 
425  if (longitude > PI)
426  longitude -= TWO_PI;
427  if (longitude < -PI)
428  longitude += TWO_PI;
429 
430  if (longitude > PI) /* force distorted values to 180, -180 degrees */
431  longitude = PI;
432  else if (longitude < -PI)
433  longitude = -PI;
434  }
435 
436  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
437 }
438 
439 
440 
441 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2