Topic Text   Topic Comments (73)   Topic Properties   Topic Information bla@bla...
Topic title: Test review Wednesday June 18, 2008 01:28:34

Download topic text | View in variable-width font | Tab width set to 4 (change to 8)

Files in topic:  
[Jump to] m:\mydoc\eckert4.c   {+384,-0}

[Add General Comment] to topic.

File m:\mydoc\eckert4.c (Revision 1.0) [Add File Comment] [Top]
 
1 /***************************************************************************/
2 /* RSC IDENTIFIER: ECKERT4
3 *
4 * ABSTRACT
5 *
6 * This component provides conversions between Geodetic coordinates
7 * (latitude and longitude in radians) and Eckert IV projection coordinates
8 * (easting and northing in meters). This projection employs a spherical
9 * Earth model. The spherical radius used is the radius of the sphere
10 * having the same area as the ellipsoid.
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 * ECK4_NO_ERROR : No errors occurred in function
20 * ECK4_LAT_ERROR : Latitude outside of valid range
21 * (-90 to 90 degrees)
22 * ECK4_LON_ERROR : Longitude outside of valid range
23 * (-180 to 360 degrees)
24 * ECK4_EASTING_ERROR : Easting outside of valid range
25 * (False_Easting +/- ~17,000,000 m,
26 * depending on ellipsoid parameters)
27 * ECK4_NORTHING_ERROR : Northing outside of valid range
28 * (False_Northing +/- 0 to 8,000,000 m,
29 * depending on ellipsoid parameters)
30 * ECK4_CENT_MER_ERROR : Central_Meridian outside of valid range
31 * (-180 to 360 degrees)
32 * ECK4_A_ERROR : Semi-major axis less than or equal to zero
33 * ECK4_INV_F_ERROR : Inverse flattening outside of valid range
34 * (250 to 350)
35 *
36 * REUSE NOTES
37 *
38 * ECKERT4 is intended for reuse by any application that performs a
39 * Eckert IV projection or its inverse.
40 *
41 * REFERENCES
42 *
43 * Further information on ECKERT4 can be found in the Reuse Manual.
44 *
45 * ECKERT4 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 * ECKERT4 has no restrictions.
57 *
58 * ENVIRONMENT
59 *
60 * ECKERT4 was tested and certified in the following environments:
61 *
62 * 1. Solaris 2.5 with GCC 2.8.1
63 * 2. MS Windows 95 with MS Visual C++ 6
64 *
65 * MODIFICATIONS
66 *
67 * Date Description
68 * ---- -----------
69 * 04/16/99 Original Code
70 *
71 */
72
73
74 /***************************************************************************/
75 /*
76 * INCLUDES
77 */
78
79 #include <math.h>
80 #include "eckert4.h"
81
82 /*
83 * math.h - Standard C math library
84 * eckert4.h - Is for prototype error checking.
85 */
86
87
88 /***************************************************************************/
89 /*
90 * DEFINES
91 */
92
93 #define PI 3.14159265358979323e0 /* PI */
94 #define PI_OVER_2 ( PI / 2.0)
95 #define TWO_PI (2.0 * PI)
96 #define NUM(Theta, SinTheta, CosTheta) (Theta + SinTheta * CosTheta + 2.0 * SinTheta)
97
98
99 /***************************************************************************/
100 /*
101 * GLOBALS
102 */
103
104 const double two_PLUS_PI_OVER_2 = (2.0 + PI / 2.0);
105
106 /* Ellipsoid Parameters, default to WGS 84 */
107 static double Eck4_a = 6378137.0; /* Semi-major axis of ellipsoid in meters */
108 static double Eck4_f = 1 / 298.257223563; /* Flattening of ellipsoid */
109 static double es2 = 0.0066943799901413800; /* Eccentricity (0.08181919084262188000) squared */
110 static double es4 = 4.4814723452405e-005; /* es2 * es2 */
111 static double es6 = 3.0000678794350e-007; /* es4 * es2 */
112
113 static double Ra0 = 2690082.6043273; /* 0.4222382 * Sperical Radius (6371007.1810824) */
114 static double Ra1 = 8451143.5741087; /* 1.3265004 * Sperical Radius (6371007.1810824) */
115
116 /* Eckert4 projection Parameters */
117 static double Eck4_Origin_Long = 0.0; /* Longitude of origin in radians */
118 static double Eck4_False_Easting = 0.0;
119 static double Eck4_False_Northing = 0.0;
120 static double Eck4_Delta_Northing = 8451144.0;
121 static double Eck4_Max_Easting = 16902288.0;
122 static double Eck4_Min_Easting = -16902288.0;
123 /*
124 * These state variables are for optimization purposes. The only function
125 * that should modify them is Set_Eckert IV_Parameters.
126 */
127
128
129 /***************************************************************************/
130 /*
131 * FUNCTIONS
132 */
133
134
135 long Set_Eckert4_Parameters(double a,
136 double f,
137 double Central_Meridian,
138 double False_Easting,
139 double False_Northing)
140
141 { /* Begin Set_Eckert4_Parameters */
142 /*
143 * The function Set_Eckert4_Parameters receives the ellipsoid parameters and
144 * projection parameters as inputs, and sets the corresponding state
145 * variables. If any errors occur, the error code(s) are returned by the
146 * function, otherwise ECK4_NO_ERROR is returned.
147 *
148 * a : Semi-major axis of ellipsoid, in meters (input)
149 * f : Flattening of ellipsoid (input)
150 * Central_Meridian : Longitude in radians at the center of (input)
151 * the projection
152 * False_Easting : A coordinate value in meters assigned to the
153 * central meridian of the projection. (input)
154 * False_Northing : A coordinate value in meters assigned to the
155 * origin latitude of the projection (input)
156 */
157
158 double Ra; /* Spherical radius */
159 double inv_f = 1 / f;
160 long Error_Code = ECK4_NO_ERROR;
161
162 if (a <= 0.0)
163 { /* Semi-major axis must be greater than zero */
164 Error_Code |= ECK4_A_ERROR;
165 }
166 if ((inv_f < 250) || (inv_f > 350))
167 { /* Inverse flattening must be between 250 and 350 */
168 Error_Code |= ECK4_INV_F_ERROR;
169 }
170 if ((Central_Meridian < -PI) || (Central_Meridian > TWO_PI))
171 { /* origin longitude out of range */
172 Error_Code |= ECK4_CENT_MER_ERROR;
173 }
174 if (!Error_Code)
175 { /* no errors */
176 Eck4_a = a;
177 Eck4_f = f;
178 es2 = 2 * Eck4_f - Eck4_f * Eck4_f;
179 es4 = es2 * es2;
180 es6 = es4 * es2;
181 /* spherical radius */
182 Ra = Eck4_a * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 / 3024.0);
183 Ra0 = 0.4222382 * Ra;
184 Ra1 = 1.3265004 * Ra;
185 if (Central_Meridian > PI)
186 Central_Meridian -= TWO_PI;
187 Eck4_Origin_Long = Central_Meridian;
188 Eck4_False_Easting = False_Easting;
189 Eck4_False_Northing = False_Northing;
190 if (Eck4_Origin_Long > 0)
191 {
192 Eck4_Max_Easting = 16808386.0;
193 Eck4_Min_Easting = -16902288.0;
194 }
195 else if (Eck4_Origin_Long < 0)
196 {
197 Eck4_Max_Easting = 16902288.0;
198 Eck4_Min_Easting = -16808386.0;
199 }
200 else
201 {
202 Eck4_Max_Easting = 16902288.0;
203 Eck4_Min_Easting = -16902288.0;
204 }
205
206 } /* End if(!Error_Code) */
207 return (Error_Code);
208 } /* End Set_Eckert4_Parameters */
209
210
211 void Get_Eckert4_Parameters(double *a,
212 double *f,
213 double *Central_Meridian,
214 double *False_Easting,
215 double *False_Northing)
216 { /* Begin Get_Eckert4_Parameters */
217 /*
218 * The function Get_Eckert4_Parameters returns the current ellipsoid
219 * parameters and Eckert IV projection parameters.
220 *
221 * a : Semi-major axis of ellipsoid, in meters (output)
222 * f : Flattening of ellipsoid (output)
223 * Central_Meridian : Longitude in radians at the center of (output)
224 * the projection
225 * False_Easting : A coordinate value in meters assigned to the
226 * central meridian of the projection. (output)
227 * False_Northing : A coordinate value in meters assigned to the
228 * origin latitude of the projection (output)
229 */
230
231 *a = Eck4_a;
232 *f = Eck4_f;
233 *Central_Meridian = Eck4_Origin_Long;
234 *False_Easting = Eck4_False_Easting;
235 *False_Northing = Eck4_False_Northing;
236 return;
237 } /* End Get_Eckert4_Parameters */
238
239
240 long Convert_Geodetic_To_Eckert4 (double Latitude,
241 double Longitude,
242 double *Easting,
243 double *Northing)
244
245 { /* Begin Convert_Geodetic_To_Eckert4 */
246 /*
247 * The function Convert_Geodetic_To_Eckert4 converts geodetic (latitude and
248 * longitude) coordinates to Eckert IV projection (easting and northing)
249 * coordinates, according to the current ellipsoid, spherical radius and
250 * Eckert IV projection parameters.
251 * If any errors occur, the error code(s) are returned by the
252 * function, otherwise ECK4_NO_ERROR is returned.
253 *
254 * Latitude : Latitude (phi) in radians (input)
255 * Longitude : Longitude (lambda) in radians (input)
256 * Easting : Easting (X) in meters (output)
257 * Northing : Northing (Y) in meters (output)
258 */
259
260 double slat = sin(Latitude);
261 double sin_theta, cos_theta;
262 double num;
263 double dlam; /* Longitude - Central Meridan */
264 double theta = Latitude / 2.0;
265 double delta_theta = 1.0;
266 double dt_tolerance = 4.85e-10; /* approximately 1/1000th of
267 an arc second or 1/10th meter */
268 long Error_Code = ECK4_NO_ERROR;
269
270 if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
271 { /* Latitude out of range */
272 Error_Code |= ECK4_LAT_ERROR;
273 }
274 if ((Longitude < -PI) || (Longitude > TWO_PI))
275 { /* Longitude out of range */
276 Error_Code|= ECK4_LON_ERROR;
277 }
278
279 if (!Error_Code)
280 { /* no errors */
281
282 dlam = Longitude - Eck4_Origin_Long;
283 if (dlam > PI)
284 {
285 dlam -= TWO_PI;
286 }
287 if (dlam < -PI)
288 {
289 dlam += TWO_PI;
290 }
291 while (fabs(delta_theta) > dt_tolerance)
292 {
293 sin_theta = sin(theta);
294 cos_theta = cos(theta);
295 num = NUM(theta, sin_theta, cos_theta);
296 delta_theta = -(num - two_PLUS_PI_OVER_2 * slat) /
297 (2.0 * cos_theta * (1.0 + cos_theta));
298 theta += delta_theta;
299 }
300 *Easting = Ra0 * dlam * (1.0 + cos(theta)) + Eck4_False_Easting;
301 *Northing = Ra1 * sin(theta) + Eck4_False_Northing;
302
303 }
304 return (Error_Code);
305
306 } /* End Convert_Geodetic_To_Eckert4 */
307
308
309 long Convert_Eckert4_To_Geodetic(double Easting,
310 double Northing,
311 double *Latitude,
312 double *Longitude)
313 { /* Begin Convert_Eckert4_To_Geodetic */
314 /*
315 * The function Convert_Eckert4_To_Geodetic converts Eckert IV projection
316 * (easting and northing) coordinates to geodetic (latitude and longitude)
317 * coordinates, according to the current ellipsoid, spherical radius and
318 * Eckert IV projection coordinates.
319 * If any errors occur, the error code(s) are returned by the
320 * function, otherwise ECK4_NO_ERROR is returned.
321 *
322 * Easting : Easting (X) in meters (input)
323 * Northing : Northing (Y) in meters (input)
324 * Latitude : Latitude (phi) in radians (output)
325 * Longitude : Longitude (lambda) in radians (output)
326 */
327
328 double theta;
329 double sin_theta, cos_theta;
330 double num;
331 double dx, dy;
332 double i;
333
334 long Error_Code = ECK4_NO_ERROR;
335
336 if ((Easting < (Eck4_False_Easting + Eck4_Min_Easting))
337 || (Easting > (Eck4_False_Easting + Eck4_Max_Easting)))
338 { /* Easting out of range */
339 Error_Code |= ECK4_EASTING_ERROR;
340 }
341 if ((Northing < (Eck4_False_Northing - Eck4_Delta_Northing))
342 || (Northing > (Eck4_False_Northing + Eck4_Delta_Northing)))
343 { /* Northing out of range */
344 Error_Code |= ECK4_NORTHING_ERROR;
345 }
346
347 if (!Error_Code)
348 {
349 dy = Northing - Eck4_False_Northing;
350 dx = Easting - Eck4_False_Easting;
351 i = dy / Ra1;
352 if (i > 1.0)
353 i = 1.0;
354 else if (i < -1.0)
355 i = -1.0;
356
357 theta = asin(i);
358 sin_theta = sin(theta);
359 cos_theta = cos(theta);
360 num = NUM(theta, sin_theta, cos_theta);
361
362 *Latitude = asin(num / two_PLUS_PI_OVER_2);
363 *Longitude = Eck4_Origin_Long + dx / (Ra0 * (1 + cos_theta));
364
365 if (*Latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
366 *Latitude = PI_OVER_2;
367 else if (*Latitude < -PI_OVER_2)
368 *Latitude = -PI_OVER_2;
369
370 if (*Longitude > PI)
371 *Longitude -= TWO_PI;
372 if (*Longitude < -PI)
373 *Longitude += TWO_PI;
374
375 if (*Longitude > PI) /* force distorted values to 180, -180 degrees */
376 *Longitude = PI;
377 else if (*Longitude < -PI)
378 *Longitude = -PI;
379
380 }
381 return (Error_Code);
382
383 } /* End Convert_Eckert4_To_Geodetic */
 
  
Legend:
Removed 
Changed
 Added

[Add General Comment] to topic.