UNCLASSIFIED
GeographicTranslator
 All Classes Namespaces Files Functions Variables Enumerations Enumerator Friends Macros
MillerCylindrical.cpp
Go to the documentation of this file.
1 // CLASSIFICATION: UNCLASSIFIED
2 
3 /***************************************************************************/
4 /* RSC IDENTIFIER: MILLER
5  *
6  * ABSTRACT
7  *
8  * This component provides conversions between Geodetic coordinates
9  * (latitude and longitude in radians) and Miller Cylindrical projection
10  * coordinates (easting and northing in meters). The Miller Cylindrical
11  * projection employs a spherical Earth model. The Spherical Radius
12  * used is the the radius of the sphere having the same area as the
13  * ellipsoid.
14  *
15  * ERROR HANDLING
16  *
17  * This component checks parameters for valid values. If an invalid value
18  * is found, the error code is combined with the current error code using
19  * the bitwise or. This combining allows multiple error codes to be
20  * returned. The possible error codes are:
21  *
22  * MILL_NO_ERROR : No errors occurred in function
23  * MILL_LAT_ERROR : Latitude outside of valid range
24  * (-90 to 90 degrees)
25  * MILL_LON_ERROR : Longitude outside of valid range
26  * (-180 to 360 degrees)
27  * MILL_EASTING_ERROR : Easting outside of valid range
28  * (False_Easting +/- ~20,000,000 m,
29  * depending on ellipsoid parameters)
30  * MILL_NORTHING_ERROR : Northing outside of valid range
31  * (False_Northing +/- ~14,000,000 m,
32  * depending on ellipsoid parameters)
33  * MILL_CENT_MER_ERROR : Central meridian outside of valid range
34  * (-180 to 360 degrees)
35  * MILL_A_ERROR : Semi-major axis less than or equal to zero
36  * MILL_INV_F_ERROR : Inverse flattening outside of valid range
37  * (250 to 350)
38  *
39  * REUSE NOTES
40  *
41  * MILLER is intended for reuse by any application that performs a
42  * Miller Cylindrical projection or its inverse.
43  *
44  * REFERENCES
45  *
46  * Further information on MILLER can be found in the Reuse Manual.
47  *
48  * MILLER originated from : U.S. Army Topographic Engineering Center
49  * Geospatial Information Division
50  * 7701 Telegraph Road
51  * Alexandria, VA 22310-3864
52  *
53  * LICENSES
54  *
55  * None apply to this component.
56  *
57  * RESTRICTIONS
58  *
59  * MILLER has no restrictions.
60  *
61  * ENVIRONMENT
62  *
63  * MILLER was tested and certified in the following environments:
64  *
65  * 1. Solaris 2.5 with GCC 2.8.1
66  * 2. Windows 95 with MS Visual C++ 6
67  *
68  * MODIFICATIONS
69  *
70  * Date Description
71  * ---- -----------
72  * 04-16-99 Original Code
73  * 03-05-07 Original C++ Code
74  *
75  */
76 
77 
78 /***************************************************************************/
79 /*
80  * INCLUDES
81  */
82 
83 #include <math.h>
84 #include "MillerCylindrical.h"
85 #include "MapProjection3Parameters.h"
86 #include "MapProjectionCoordinates.h"
87 #include "GeodeticCoordinates.h"
88 #include "CoordinateConversionException.h"
89 #include "ErrorMessages.h"
90 
91 /*
92  * math.h - Standard C math library
93  * MillerCylindrical.h - Is for prototype error checking
94  * MapProjectionCoordinates.h - defines map projection coordinates
95  * GeodeticCoordinates.h - defines geodetic coordinates
96  * CoordinateConversionException.h - Exception handler
97  * ErrorMessages.h - Contains exception messages
98  */
99 
100 
101 using namespace MSP::CCS;
102 
103 
104 /***************************************************************************/
105 /* DEFINES
106  *
107  */
108 
109 const double PI = 3.14159265358979323e0; /* PI */
110 const double PI_OVER_2 = ( PI / 2.0);
111 const double TWO_PI = (2.0 * PI);
112 
113 
114 /************************************************************************/
115 /* FUNCTIONS
116  *
117  */
118 
119 MillerCylindrical::MillerCylindrical( double ellipsoidSemiMajorAxis, double ellipsoidFlattening, double centralMeridian, double falseEasting, double falseNorthing ) :
120  CoordinateSystem(),
121  es2( 0.0066943799901413800 ),
122  es4( 4.4814723452405e-005 ),
123  es6( 3.0000678794350e-007 ),
124  Ra( 6371007.1810824 ),
125  Mill_Origin_Long( 0.0 ),
126  Mill_False_Easting( 0.0 ),
127  Mill_False_Northing( 0.0 ),
128  Mill_Delta_Northing( 14675058.0 ),
129  Mill_Max_Easting( 20015110.0 ),
130  Mill_Min_Easting( -20015110.0 )
131 {
132 /*
133  * The constructor receives the ellipsoid parameters and
134  * Miller Cylindrical projcetion parameters as inputs, and sets the corresponding
135  * state variables. If any errors occur, an exception is thrown with a description
136  * of the error.
137  *
138  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (input)
139  * ellipsoidFlattening : Flattening of ellipsoid (input)
140  * centralMeridian : Longitude in radians at the center of (input)
141  * the projection
142  * falseEasting : A coordinate value in meters assigned to the
143  * central meridian of the projection. (input)
144  * falseNorthing : A coordinate value in meters assigned to the
145  * origin latitude of the projection (input)
146  */
147 
148  double inv_f = 1 / ellipsoidFlattening;
149 
150  if (ellipsoidSemiMajorAxis <= 0.0)
151  { /* Semi-major axis must be greater than zero */
152  throw CoordinateConversionException( ErrorMessages::semiMajorAxis );
153  }
154  if ((inv_f < 250) || (inv_f > 350))
155  { /* Inverse flattening must be between 250 and 350 */
156  throw CoordinateConversionException( ErrorMessages::ellipsoidFlattening );
157  }
158  if ((centralMeridian < -PI) || (centralMeridian > TWO_PI))
159  { /* origin longitude out of range */
160  throw CoordinateConversionException( ErrorMessages::centralMeridian );
161  }
162 
163  semiMajorAxis = ellipsoidSemiMajorAxis;
164  flattening = ellipsoidFlattening;
165 
166  es2 = 2 * flattening - flattening * flattening;
167  es4 = es2 * es2;
168  es6 = es4 * es2;
169  /* spherical radius */
170  Ra = semiMajorAxis * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
171  if (centralMeridian > PI)
172  centralMeridian -= TWO_PI;
173  Mill_Origin_Long = centralMeridian;
174  Mill_False_Easting = falseEasting;
175  Mill_False_Northing = falseNorthing;
176  if (Mill_Origin_Long > 0)
177  {
178  Mill_Max_Easting = 19903915.0;
179  Mill_Min_Easting = -20015110.0;
180  }
181  else if (Mill_Origin_Long < 0)
182  {
183  Mill_Max_Easting = 20015110.0;
184  Mill_Min_Easting = -19903915.0;
185  }
186  else
187  {
188  Mill_Max_Easting = 20015110.0;
189  Mill_Min_Easting = -20015110.0;
190  }
191 }
192 
193 
194 MillerCylindrical::MillerCylindrical( const MillerCylindrical &mc )
195 {
196  semiMajorAxis = mc.semiMajorAxis;
197  flattening = mc.flattening;
198  es2 = mc.es2;
199  es4 = mc.es4;
200  es6 = mc.es6;
201  Ra = mc.Ra;
202  Mill_Origin_Long = mc.Mill_Origin_Long;
203  Mill_False_Easting = mc.Mill_False_Easting;
204  Mill_False_Northing = mc.Mill_False_Northing;
205  Mill_Delta_Northing = mc.Mill_Delta_Northing;
206  Mill_Max_Easting = mc.Mill_Max_Easting;
207  Mill_Min_Easting = mc.Mill_Min_Easting;
208 }
209 
210 
211 MillerCylindrical::~MillerCylindrical()
212 {
213 }
214 
215 
216 MillerCylindrical& MillerCylindrical::operator=( const MillerCylindrical &mc )
217 {
218  if( this != &mc )
219  {
220  semiMajorAxis = mc.semiMajorAxis;
221  flattening = mc.flattening;
222  es2 = mc.es2;
223  es4 = mc.es4;
224  es6 = mc.es6;
225  Ra = mc.Ra;
226  Mill_Origin_Long = mc.Mill_Origin_Long;
227  Mill_False_Easting = mc.Mill_False_Easting;
228  Mill_False_Northing = mc.Mill_False_Northing;
229  Mill_Delta_Northing = mc.Mill_Delta_Northing;
230  Mill_Max_Easting = mc.Mill_Max_Easting;
231  Mill_Min_Easting = mc.Mill_Min_Easting;
232  }
233 
234  return *this;
235 }
236 
237 
238 MapProjection3Parameters* MillerCylindrical::getParameters() const
239 {
240 /*
241  * The function getParameters returns the current ellipsoid
242  * parameters and Miller Cylindrical projection parameters.
243  *
244  * ellipsoidSemiMajorAxis : Semi-major axis of ellipsoid, in meters (output)
245  * ellipsoidFlattening : Flattening of ellipsoid (output)
246  * centralMeridian : Longitude in radians at the center of (output)
247  * the projection
248  * falseEasting : A coordinate value in meters assigned to the
249  * central meridian of the projection. (output)
250  * falseNorthing : A coordinate value in meters assigned to the
251  * origin latitude of the projection (output)
252  */
253 
254  return new MapProjection3Parameters( CoordinateType::millerCylindrical, Mill_Origin_Long, Mill_False_Easting, Mill_False_Northing );
255 }
256 
257 
258 MSP::CCS::MapProjectionCoordinates* MillerCylindrical::convertFromGeodetic( MSP::CCS::GeodeticCoordinates* geodeticCoordinates )
259 {
260 /*
261  * The function convertFromGeodetic converts geodetic (latitude and
262  * longitude) coordinates to Miller Cylindrical projection (easting and northing)
263  * coordinates, according to the current ellipsoid and Miller Cylindrical projection
264  * parameters. If any errors occur, an exception is thrown with a description
265  * of the error.
266  *
267  * longitude : Longitude (lambda) in radians (input)
268  * latitude : Latitude (phi) in radians (input)
269  * easting : Easting (X) in meters (output)
270  * northing : Northing (Y) in meters (output)
271  */
272 
273  double dlam; /* Longitude - Central Meridan */
274  double longitude = geodeticCoordinates->longitude();
275  double latitude = geodeticCoordinates->latitude();
276  double slat = sin(0.8 * latitude);
277 
278  if ((latitude < -PI_OVER_2) || (latitude > PI_OVER_2))
279  { /* Latitude out of range */
280  throw CoordinateConversionException( ErrorMessages::latitude );
281  }
282  if ((longitude < -PI) || (longitude > TWO_PI))
283  { /* Longitude out of range */
284  throw CoordinateConversionException( ErrorMessages::longitude );
285  }
286 
287  dlam = longitude - Mill_Origin_Long;
288  if (dlam > PI)
289  {
290  dlam -= TWO_PI;
291  }
292  if (dlam < -PI)
293  {
294  dlam += TWO_PI;
295  }
296  double easting = Ra * dlam + Mill_False_Easting;
297  double northing = (Ra / 1.6) * log((1.0 + slat) /
298  (1.0 - slat)) + Mill_False_Northing;
299 
300  return new MapProjectionCoordinates( CoordinateType::millerCylindrical, easting, northing );
301 }
302 
303 
304 MSP::CCS::GeodeticCoordinates* MillerCylindrical::convertToGeodetic(
305  MSP::CCS::MapProjectionCoordinates* mapProjectionCoordinates )
306 {
307 /*
308  * The function convertToGeodetic converts Miller Cylindrical projection
309  * (easting and northing) coordinates to geodetic (latitude and longitude)
310  * coordinates, according to the current ellipsoid and Miller Cylindrical projection
311  * coordinates. If any errors occur, an exception is thrown with a description
312  * of the error.
313  *
314  * easting : Easting (X) in meters (input)
315  * northing : Northing (Y) in meters (input)
316  * longitude : Longitude (lambda) in radians (output)
317  * latitude : Latitude (phi) in radians (output)
318  */
319 
320  double dx, dy;
321 
322  double easting = mapProjectionCoordinates->easting();
323  double northing = mapProjectionCoordinates->northing();
324 
325  if ((easting < (Mill_False_Easting + Mill_Min_Easting))
326  || (easting > (Mill_False_Easting + Mill_Max_Easting)))
327  { /* Easting out of range */
328  throw CoordinateConversionException( ErrorMessages::easting );
329  }
330  if ((northing < (Mill_False_Northing - Mill_Delta_Northing)) ||
331  (northing > (Mill_False_Northing + Mill_Delta_Northing) ))
332  { /* Northing out of range */
333  throw CoordinateConversionException( ErrorMessages::northing );
334  }
335 
336  dy = northing - Mill_False_Northing;
337  dx = easting - Mill_False_Easting;
338  double latitude = atan(sinh(0.8 * dy / Ra)) / 0.8;
339  double longitude = Mill_Origin_Long + dx / Ra;
340 
341  if (latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
342  latitude = PI_OVER_2;
343  else if (latitude < -PI_OVER_2)
344  latitude = -PI_OVER_2;
345 
346  if (longitude > PI)
347  longitude -= TWO_PI;
348  if (longitude < -PI)
349  longitude += TWO_PI;
350 
351  if (longitude > PI) /* force distorted values to 180, -180 degrees */
352  longitude = PI;
353  else if (longitude < -PI)
354  longitude = -PI;
355 
356  return new GeodeticCoordinates( CoordinateType::geodetic, longitude, latitude );
357 }
358 
359 
360 
361 // CLASSIFICATION: UNCLASSIFIED
Generated on Tue Feb 16 2016 14:54:02 for GeographicTranslator by doxygen 1.8.2