1 /**********************************************************
2 * Version $Id: grinten.c 911 2011-02-14 16:38:15Z reklov_w $
3 *********************************************************/
4 /***************************************************************************/
5 /* RSC IDENTIFIER: VAN DER GRINTEN
6 *
7 * ABSTRACT
8 *
9 * This component provides conversions between Geodetic coordinates
10 * (latitude and longitude in radians) and Van Der Grinten projection
11 * coordinates (easting and northing in meters). The Van Der Grinten
12 * projection employs a spherical Earth model. The Spherical Radius
13 * used is the the radius of the sphere having the same area as the
14 * ellipsoid.
15 *
16 * ERROR HANDLING
17 *
18 * This component checks parameters for valid values. If an invalid value
19 * is found, the error code is combined with the current error code using
20 * the bitwise or. This combining allows multiple error codes to be
21 * returned. The possible error codes are:
22 *
23 * GRIN_NO_ERROR : No errors occurred in function
24 * GRIN_LAT_ERROR : Latitude outside of valid range
25 * (-90 to 90 degrees)
26 * GRIN_LON_ERROR : Longitude outside of valid range
27 * (-180 to 360 degrees)
28 * GRIN_EASTING_ERROR : Easting outside of valid range
29 * (False_Easting +/- ~20,000,000 m,
30 * depending on ellipsoid parameters)
31 * GRIN_NORTHING_ERROR : Northing outside of valid range
32 * (False_Northing +/- ~20,000,000 m,
33 * depending on ellipsoid parameters)
34 * GRIN_RADIUS_ERROR : Coordinates too far from pole,
35 * depending on ellipsoid and
36 * projection parameters
37 * GRIN_CENT_MER_ERROR : Central meridian outside of valid range
38 * (-180 to 360 degrees)
39 * GRIN_A_ERROR : Semi-major axis less than or equal to zero
40 * GRIN_INV_F_ERROR : Inverse flattening outside of valid range
41 * (250 to 350)
42 *
43 * REUSE NOTES
44 *
45 * VAN DER GRINTEN is intended for reuse by any application that performs a
46 * Van Der Grinten projection or its inverse.
47 *
48 * REFERENCES
49 *
50 * Further information on VAN DER GRINTEN can be found in the Reuse Manual.
51 *
52 * VAN DER GRINTEN originated from : U.S. Army Topographic Engineering Center
53 * Geospatial Information Division
54 * 7701 Telegraph Road
55 * Alexandria, VA 22310-3864
56 *
57 * LICENSES
58 *
59 * None apply to this component.
60 *
61 * RESTRICTIONS
62 *
63 * VAN DER GRINTEN has no restrictions.
64 *
65 * ENVIRONMENT
66 *
67 * VAN DER GRINTEN was tested and certified in the following environments:
68 *
69 * 1. Solaris 2.5 with GCC, version 2.8.1
70 * 2. Windows 95 with MS Visual C++, version 6
71 *
72 * MODIFICATIONS
73 *
74 * Date Description
75 * ---- -----------
76 * 10-02-97 Original Code
77 *
78 */
79
80
81 /***************************************************************************/
82 /*
83 * INCLUDES
84 */
85
86 #include <math.h>
87 #include "grinten.h"
88
89 /*
90 * math.h - Standard C math library
91 * grinten.h - Is for prototype error checking
92 */
93
94
95 /***************************************************************************/
96 /*
97 * DEFINES
98 */
99
100 #define PI 3.14159265358979323e0 /* PI */
101 #define PI_OVER_2 ( PI / 2.0)
102 #define MAX_LAT ( 90 * (PI / 180.0) ) /* 90 degrees in radians */
103 #define TWO_PI (2.0 * PI)
104 #define FLOAT_EQ(x,v,epsilon) (((v - epsilon) < x) && (x < (v + epsilon)))
105
106
107 /***************************************************************************/
108 /*
109 * GLOBALS
110 */
111
112 const double TWO_OVER_PI = (2.0 / PI);
113 const double PI_OVER_3 = (PI / 3.0);
114 const double ONE_THIRD = (1.0 / 3.0);
115
116 /* Ellipsoid Parameters, default to WGS 84 */
117 static double Grin_a = 6378137.0; /* Semi-major axis of ellipsoid in meters */
118 static double Grin_f = 1 / 298.257223563; /* Flattening of ellipsoid */
119 static double es2 = 0.0066943799901413800; /* Eccentricity (0.08181919084262188000) squared */
120 static double es4 = 4.4814723452405e-005; /* es2 * es2 */
121 static double es6 = 3.0000678794350e-007; /* es4 * es2 */
122 static double Ra = 6371007.1810824; /* Spherical Radius */
123 static double PI_Ra = 20015109.356056;
124
125 /* Van Der Grinten projection Parameters */
126 static double Grin_Origin_Long = 0.0; /* Longitude of origin in radians */
127 static double Grin_False_Easting = 0.0;
128 static double Grin_False_Northing = 0.0;
129 /*
130 * These state variables are for optimization purposes. The only function
131 * that should modify them is Set_Van_der_Grinten_Parameters.
132 */
133
134
135 /***************************************************************************/
136 /*
137 * FUNCTIONS
138 */
139
140
Set_Van_der_Grinten_Parameters(double a,double f,double Central_Meridian,double False_Easting,double False_Northing)141 long Set_Van_der_Grinten_Parameters(double a,
142 double f,
143 double Central_Meridian,
144 double False_Easting,
145 double False_Northing)
146
147 { /* BEGIN Set_Van_der_Grinten_Parameters */
148 /*
149 * The function Set_Van_der_Grinten_Parameters receives the ellipsoid parameters and
150 * projection parameters as inputs, and sets the corresponding state
151 * variables. If any errors occur, the error code(s) are returned by the function,
152 * otherwise Grin_NO_ERROR is returned.
153 *
154 * a : Semi-major axis of ellipsoid, in meters (input)
155 * f : Flattening of ellipsoid (input)
156 * Central_Meridian : Longitude in radians at the center of (input)
157 * the projection
158 * False_Easting : A coordinate value in meters assigned to the
159 * central meridian of the projection. (input)
160 * False_Northing : A coordinate value in meters assigned to the
161 * origin latitude of the projection (input)
162 */
163
164 double inv_f = 1 / f;
165 long Error_Code = GRIN_NO_ERROR;
166
167 if (a <= 0.0)
168 { /* Semi-major axis must be greater than zero */
169 Error_Code |= GRIN_A_ERROR;
170 }
171 if ((inv_f < 250) || (inv_f > 350))
172 { /* Inverse flattening must be between 250 and 350 */
173 Error_Code |= GRIN_INV_F_ERROR;
174 }
175 if ((Central_Meridian < -PI) || (Central_Meridian > TWO_PI))
176 { /* origin longitude out of range */
177 Error_Code |= GRIN_CENT_MER_ERROR;
178 }
179 if (!Error_Code)
180 { /* no errors */
181 Grin_a = a;
182 Grin_f = f;
183 es2 = 2 * Grin_f - Grin_f * Grin_f;
184 es4 = es2 * es2;
185 es6 = es4 * es2;
186 /* spherical radius */
187 Ra = Grin_a * (1.0 - es2 / 6.0 - 17.0 * es4 / 360.0 - 67.0 * es6 /3024.0);
188 PI_Ra = PI * Ra;
189 if (Central_Meridian > PI)
190 Central_Meridian -= TWO_PI;
191 Grin_Origin_Long = Central_Meridian;
192 Grin_False_Easting = False_Easting;
193 Grin_False_Northing = False_Northing;
194
195 } /* END OF if(!Error_Code) */
196 return (Error_Code);
197 } /* END OF Set_Van_der_Grinten_Parameters */
198
199
Get_Van_der_Grinten_Parameters(double * a,double * f,double * Central_Meridian,double * False_Easting,double * False_Northing)200 void Get_Van_der_Grinten_Parameters(double *a,
201 double *f,
202 double *Central_Meridian,
203 double *False_Easting,
204 double *False_Northing)
205
206 { /* BEGIN Get_Van_der_Grinten_Parameters */
207 /*
208 * The function Get_Van_der_Grinten_Parameters returns the current ellipsoid
209 * parameters, and Van Der Grinten projection parameters.
210 *
211 * a : Semi-major axis of ellipsoid, in meters (output)
212 * f : Flattening of ellipsoid (output)
213 * Central_Meridian : Longitude in radians at the center of (output)
214 * the projection
215 * False_Easting : A coordinate value in meters assigned to the
216 * central meridian of the projection. (output)
217 * False_Northing : A coordinate value in meters assigned to the
218 * origin latitude of the projection (output)
219 */
220
221 *a = Grin_a;
222 *f = Grin_f;
223 *Central_Meridian = Grin_Origin_Long;
224 *False_Easting = Grin_False_Easting;
225 *False_Northing = Grin_False_Northing;
226 return;
227 } /* END OF Get_Van_der_Grinten_Parameters */
228
229
Convert_Geodetic_To_Van_der_Grinten(double Latitude,double Longitude,double * Easting,double * Northing)230 long Convert_Geodetic_To_Van_der_Grinten (double Latitude,
231 double Longitude,
232 double *Easting,
233 double *Northing)
234
235 { /* BEGIN Convert_Geodetic_To_Van_der_Grinten */
236 /*
237 * The function Convert_Geodetic_To_Van_der_Grinten converts geodetic (latitude and
238 * longitude) coordinates to Van Der Grinten projection (easting and northing)
239 * coordinates, according to the current ellipsoid and Van Der Grinten projection
240 * parameters. If any errors occur, the error code(s) are returned by the
241 * function, otherwise GRIN_NO_ERROR is returned.
242 *
243 * Latitude : Latitude (phi) in radians (input)
244 * Longitude : Longitude (lambda) in radians (input)
245 * Easting : Easting (X) in meters (output)
246 * Northing : Northing (Y) in meters (output)
247 */
248
249 double dlam; /* Longitude - Central Meridan */
250 double aa, aasqr;
251 double gg;
252 double pp, ppsqr;
253 double gg_MINUS_ppsqr, ppsqr_PLUS_aasqr;
254 double in_theta;
255 double theta;
256 double sin_theta, cos_theta;
257 double qq;
258 long Error_Code = GRIN_NO_ERROR;
259
260 if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
261 { /* Latitude out of range */
262 Error_Code |= GRIN_LAT_ERROR;
263 }
264 if ((Longitude < -PI) || (Longitude > TWO_PI))
265 { /* Longitude out of range */
266 Error_Code|= GRIN_LON_ERROR;
267 }
268
269 if (!Error_Code)
270 { /* no errors */
271
272 dlam = Longitude - Grin_Origin_Long;
273 if (dlam > PI)
274 {
275 dlam -= TWO_PI;
276 }
277 if (dlam < -PI)
278 {
279 dlam += TWO_PI;
280 }
281
282 if (Latitude == 0.0)
283 {
284 *Easting = Ra * dlam + Grin_False_Easting;
285 *Northing = 0.0;
286 }
287 else if (dlam == 0.0 || FLOAT_EQ(Latitude,MAX_LAT,.00001) || FLOAT_EQ(Latitude,-MAX_LAT,.00001))
288 {
289 in_theta = fabs(TWO_OVER_PI * Latitude);
290
291 if (in_theta > 1.0)
292 in_theta = 1.0;
293 else if (in_theta < -1.0)
294 in_theta = -1.0;
295
296 theta = asin(in_theta);
297 *Easting = 0.0;
298 *Northing = PI_Ra * tan(theta / 2) + Grin_False_Northing;
299 if (Latitude < 0.0)
300 *Northing *= -1.0;
301 }
302 else
303 {
304 aa = 0.5 * fabs(PI / dlam - dlam / PI);
305 in_theta = fabs(TWO_OVER_PI * Latitude);
306
307 if (in_theta > 1.0)
308 in_theta = 1.0;
309 else if (in_theta < -1.0)
310 in_theta = -1.0;
311
312 theta = asin(in_theta);
313 sin_theta = sin(theta);
314 cos_theta = cos(theta);
315 gg = cos_theta / (sin_theta + cos_theta - 1);
316 pp = gg * (2 / sin_theta - 1);
317 aasqr = aa * aa;
318 ppsqr = pp * pp;
319 gg_MINUS_ppsqr = gg - ppsqr;
320 ppsqr_PLUS_aasqr = ppsqr + aasqr;
321 qq = aasqr + gg;
322 *Easting = PI_Ra * (aa * (gg_MINUS_ppsqr) +
323 sqrt(aasqr * (gg_MINUS_ppsqr) * (gg_MINUS_ppsqr) -
324 (ppsqr_PLUS_aasqr) * (gg * gg - ppsqr))) /
325 (ppsqr_PLUS_aasqr) + Grin_False_Easting;
326 if (dlam < 0.0)
327 *Easting *= -1.0;
328 *Northing = PI_Ra * (pp * qq - aa * sqrt ((aasqr + 1) * (ppsqr_PLUS_aasqr) - qq * qq)) /
329 (ppsqr_PLUS_aasqr) + Grin_False_Northing;
330 if (Latitude < 0.0)
331 *Northing *= -1.0;
332 }
333 }
334 return (Error_Code);
335
336 } /* END OF Convert_Geodetic_To_Van_der_Grinten */
337
338
Convert_Van_der_Grinten_To_Geodetic(double Easting,double Northing,double * Latitude,double * Longitude)339 long Convert_Van_der_Grinten_To_Geodetic(double Easting,
340 double Northing,
341 double *Latitude,
342 double *Longitude)
343 { /* BEGIN Convert_Van_der_Grinten_To_Geodetic */
344 /*
345 * The function Convert_Van_der_Grinten_To_Geodetic converts Grinten projection
346 * (easting and northing) coordinates to geodetic (latitude and longitude)
347 * coordinates, according to the current ellipsoid and Grinten projection
348 * coordinates. If any errors occur, the error code(s) are returned by the
349 * function, otherwise GRIN_NO_ERROR is returned.
350 *
351 * Easting : Easting (X) in meters (input)
352 * Northing : Northing (Y) in meters (input)
353 * Latitude : Latitude (phi) in radians (output)
354 * Longitude : Longitude (lambda) in radians (output)
355 */
356
357 double dx, dy;
358 double xx, xxsqr;
359 double yy, yysqr, two_yysqr;
360 double xxsqr_PLUS_yysqr;
361 double c1;
362 double c2;
363 double c3, c3sqr;
364 double c2_OVER_3c3;
365 double dd;
366 double a1;
367 double m1;
368 double i;
369 double theta1;
370 double temp;
371 const double epsilon = 1.0e-2;
372
373 long Error_Code = GRIN_NO_ERROR;
374
375 if ((Easting > (Grin_False_Easting + PI_Ra + epsilon)) ||
376 (Easting < (Grin_False_Easting - PI_Ra - epsilon)))
377 { /* Easting out of range */
378 Error_Code |= GRIN_EASTING_ERROR;
379 }
380 if ((Northing > (Grin_False_Northing + PI_Ra + epsilon)) ||
381 (Northing < (Grin_False_Northing - PI_Ra - epsilon)))
382 { /* Northing out of range */
383 Error_Code |= GRIN_NORTHING_ERROR;
384 }
385 if (!Error_Code)
386 {
387 temp = sqrt(Easting * Easting + Northing * Northing);
388
389 if ((temp > (Grin_False_Easting + PI_Ra + epsilon)) ||
390 (temp > (Grin_False_Northing + PI_Ra + epsilon)) ||
391 (temp < (Grin_False_Easting - PI_Ra - epsilon)) ||
392 (temp < (Grin_False_Northing - PI_Ra - epsilon)))
393 { /* Point is outside of projection area */
394 Error_Code |= GRIN_RADIUS_ERROR;
395 }
396 }
397
398 if (!Error_Code)
399 {
400 dy = Northing - Grin_False_Northing;
401 dx = Easting - Grin_False_Easting;
402 xx = dx / PI_Ra;
403 yy = dy / PI_Ra;
404 xxsqr = xx * xx;
405 yysqr = yy * yy;
406 xxsqr_PLUS_yysqr = xxsqr + yysqr;
407 two_yysqr = 2 * yysqr;
408
409 if (Northing == 0.0)
410 *Latitude = 0.0;
411
412 else
413 {
414 c1 = - fabs(yy) * (1 + xxsqr_PLUS_yysqr);
415 c2 = c1 - two_yysqr + xxsqr;
416 c3 = - 2 * c1 + 1 + two_yysqr + (xxsqr_PLUS_yysqr) * (xxsqr_PLUS_yysqr);
417 c2_OVER_3c3 = c2 / (3.0 * c3);
418 c3sqr = c3 * c3;
419 dd = yysqr / c3 + ((2 * c2 * c2 * c2) / (c3sqr * c3) - (9 * c1 * c2) / (c3sqr)) / 27;
420 a1 = (c1 - c2 * c2_OVER_3c3) /c3;
421 m1 = 2 * sqrt(-ONE_THIRD * a1);
422 i = 3 * dd/ (a1 * m1);
423 if ((i > 1.0)||(i < -1.0))
424 *Latitude = MAX_LAT;
425 else
426 {
427 theta1 = ONE_THIRD * acos(3 * dd / (a1 * m1));
428 *Latitude = PI * (-m1 * cos(theta1 + PI_OVER_3) - c2_OVER_3c3);
429 }
430 }
431 if (Northing < 0.0)
432 *Latitude *= -1.0;
433
434 if (xx == 0.0)
435 *Longitude = Grin_Origin_Long;
436 else
437 {
438 *Longitude = PI * (xxsqr_PLUS_yysqr - 1 +
439 sqrt(1 + (2 * xxsqr - two_yysqr) + (xxsqr_PLUS_yysqr) * (xxsqr_PLUS_yysqr))) /
440 (2 * xx) + Grin_Origin_Long;
441 }
442 if (*Latitude > PI_OVER_2) /* force distorted values to 90, -90 degrees */
443 *Latitude = PI_OVER_2;
444 else if (*Latitude < -PI_OVER_2)
445 *Latitude = -PI_OVER_2;
446
447 if (*Longitude > PI)
448 *Longitude -= TWO_PI;
449 if (*Longitude < -PI)
450 *Longitude += TWO_PI;
451
452 if (*Longitude > PI) /* force distorted values to 180, -180 degrees */
453 *Longitude = PI;
454 else if (*Longitude < -PI)
455 *Longitude = -PI;
456
457 }
458 return (Error_Code);
459
460 } /* END OF Convert_Van_der_Grinten_To_Geodetic */
461
462