1
/***************************************************************************/
2
/* RSC IDENTIFIER: ECKERT4
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.
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:
19
* ECK4_NO_ERROR : No errors occurred in function
20
* ECK4_LAT_ERROR : Latitude outside of valid range
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
38
* ECKERT4 is intended for reuse by any application that performs a
39
* Eckert IV projection or its inverse.
43
* Further information on ECKERT4 can be found in the Reuse Manual.
45
* ECKERT4 originated from : U.S. Army Topographic Engineering Center
46
* Geospatial Information Division
48
* Alexandria, VA 22310-3864
52
* None apply to this component.
56
* ECKERT4 has no restrictions.
60
* ECKERT4 was tested and certified in the following environments:
62
* 1. Solaris 2.5 with GCC 2.8.1
63
* 2. MS Windows 95 with MS Visual C++ 6
69
* 04/16/99 Original Code
74
/***************************************************************************/
83
* math.h - Standard C math library
84
* eckert4.h - Is for prototype error checking.
88
/***************************************************************************/
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)
99
/***************************************************************************/
104
const double two_PLUS_PI_OVER_2 = (2.0 + PI / 2.0);
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 */
113
static double Ra0 = 2690082.6043273; /* 0.4222382 * Sperical Radius (6371007.1810824) */
114
static double Ra1 = 8451143.5741087; /* 1.3265004 * Sperical Radius (6371007.1810824) */
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;
124
* These state variables are for optimization purposes. The only function
125
* that should modify them is Set_Eckert IV_Parameters.
129
/***************************************************************************/
135
long Set_Eckert4_Parameters(double a,
137
double Central_Meridian,
138
double False_Easting,
139
double False_Northing)
141
{ /* Begin Set_Eckert4_Parameters */
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.
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)
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)
158
double Ra; /* Spherical radius */
159
double inv_f = 1 / f;
160
long Error_Code = ECK4_NO_ERROR;
163
{ /* Semi-major axis must be greater than zero */
164
Error_Code |= ECK4_A_ERROR;
166
if ((inv_f < 250) || (inv_f > 350))
167
{ /* Inverse flattening must be between 250 and 350 */
168
Error_Code |= ECK4_INV_F_ERROR;
170
if ((Central_Meridian < -PI) || (Central_Meridian > TWO_PI))
171
{ /* origin longitude out of range */
172
Error_Code |= ECK4_CENT_MER_ERROR;
178
es2 = 2 * Eck4_f - Eck4_f * Eck4_f;
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)
192
Eck4_Max_Easting = 16808386.0;
193
Eck4_Min_Easting = -16902288.0;
195
else if (Eck4_Origin_Long < 0)
197
Eck4_Max_Easting = 16902288.0;
198
Eck4_Min_Easting = -16808386.0;
202
Eck4_Max_Easting = 16902288.0;
203
Eck4_Min_Easting = -16902288.0;
206
} /* End if(!Error_Code) */
208
} /* End Set_Eckert4_Parameters */
211
void Get_Eckert4_Parameters(double *a,
213
double *Central_Meridian,
214
double *False_Easting,
215
double *False_Northing)
216
{ /* Begin Get_Eckert4_Parameters */
218
* The function Get_Eckert4_Parameters returns the current ellipsoid
219
* parameters and Eckert IV projection parameters.
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)
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)
233
*Central_Meridian = Eck4_Origin_Long;
234
*False_Easting = Eck4_False_Easting;
235
*False_Northing = Eck4_False_Northing;
237
} /* End Get_Eckert4_Parameters */
240
long Convert_Geodetic_To_Eckert4 (double Latitude,
245
{ /* Begin Convert_Geodetic_To_Eckert4 */
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.
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)
260
double slat = sin(Latitude);
261
double sin_theta, cos_theta;
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;
270
if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
271
{ /* Latitude out of range */
272
Error_Code |= ECK4_LAT_ERROR;
274
if ((Longitude < -PI) || (Longitude > TWO_PI))
275
{ /* Longitude out of range */
276
Error_Code|= ECK4_LON_ERROR;
282
dlam = Longitude - Eck4_Origin_Long;
291
while (fabs(delta_theta) > dt_tolerance)
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;
300
*Easting = Ra0 * dlam * (1.0 + cos(theta)) + Eck4_False_Easting;
301
*Northing = Ra1 * sin(theta) + Eck4_False_Northing;
306
} /* End Convert_Geodetic_To_Eckert4 */
309
long Convert_Eckert4_To_Geodetic(double Easting,
313
{ /* Begin Convert_Eckert4_To_Geodetic */
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.
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)
329
double sin_theta, cos_theta;
334
long Error_Code = ECK4_NO_ERROR;
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;
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;
349
dy = Northing - Eck4_False_Northing;
350
dx = Easting - Eck4_False_Easting;
358
sin_theta = sin(theta);
359
cos_theta = cos(theta);
360
num = NUM(theta, sin_theta, cos_theta);
362
*Latitude = asin(num / two_PLUS_PI_OVER_2);
363
*Longitude = Eck4_Origin_Long + dx / (Ra0 * (1 + cos_theta));
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;
371
*Longitude -= TWO_PI;
372
if (*Longitude < -PI)
373
*Longitude += TWO_PI;
375
if (*Longitude > PI) /* force distorted values to 180, -180 degrees */
377
else if (*Longitude < -PI)
383
} /* End Convert_Eckert4_To_Geodetic */