1
/***************************************************************************/
2
/* RSC IDENTIFIER: GEOCENTRIC
6
* This component provides conversions between Geodetic coordinates (latitude,
7
* longitude in radians and height in meters) and Geocentric coordinates
12
* This component checks parameters for valid values. If an invalid value
13
* is found, the error code is combined with the current error code using
14
* the bitwise or. This combining allows multiple error codes to be
15
* returned. The possible error codes are:
17
* GEOCENT_NO_ERROR : No errors occurred in function
18
* GEOCENT_LAT_ERROR : Latitude out of valid range
20
* GEOCENT_LON_ERROR : Longitude out of valid range
21
* (-180 to 360 degrees)
22
* GEOCENT_A_ERROR : Semi-major axis less than or equal to zero
23
* GEOCENT_INV_F_ERROR : Inverse flattening outside of valid range
29
* GEOCENTRIC is intended for reuse by any application that performs
30
* coordinate conversions between geodetic coordinates and geocentric
36
* An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
37
* Ralph Toms, February 1996 UCRL-JC-123138.
39
* Further information on GEOCENTRIC can be found in the Reuse Manual.
41
* GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
42
* Geospatial Information Division
44
* Alexandria, VA 22310-3864
48
* None apply to this component.
52
* GEOCENTRIC has no restrictions.
56
* GEOCENTRIC was tested and certified in the following environments:
58
* 1. Solaris 2.5 with GCC version 2.8.1
59
* 2. Windows 95 with MS Visual C++ version 6
65
* 25-02-97 Original Code
70
/***************************************************************************/
77
* math.h - is needed for calls to sin, cos, tan and sqrt.
78
* geocent.h - is needed for Error codes and prototype error checking.
82
/***************************************************************************/
86
#define PI 3.14159265358979323e0
87
#define PI_OVER_2 (PI / 2.0e0)
90
#define COS_67P5 0.38268343236508977 /* cosine of 67.5 degrees */
91
#define AD_C 1.0026000 /* Toms region 1 constant */
94
/***************************************************************************/
98
/* Ellipsoid parameters, default to WGS 84 */
99
double Geocent_a = 6378137.0; /* Semi-major axis of ellipsoid in meters */
100
double Geocent_f = 1 / 298.257223563; /* Flattening of ellipsoid */
102
double Geocent_e2 = 0.0066943799901413800; /* Eccentricity squared */
103
double Geocent_ep2 = 0.00673949675658690300; /* 2nd eccentricity squared */
105
* These state variables are for optimization purposes. The only function
106
* that should modify them is Set_Geocentric_Parameters.
110
/***************************************************************************/
116
long Set_Geocentric_Parameters (double a,
118
{ /* BEGIN Set_Geocentric_Parameters */
120
* The function Set_Geocentric_Parameters receives the ellipsoid parameters
121
* as inputs and sets the corresponding state variables.
123
* a : Semi-major axis of ellipsoid, in meters. (input)
124
* f : Flattening of ellipsoid. (input)
127
double inv_f = 1 / f;
128
long Error_Code = GEOCENT_NO_ERROR;
131
Error_Code |= GEOCENT_A_ERROR;
132
if ((inv_f < 250) || (inv_f > 350))
133
{ /* Inverse flattening must be between 250 and 350 */
134
Error_Code |= GEOCENT_INV_F_ERROR;
140
Geocent_e2 = 2 * Geocent_f - Geocent_f * Geocent_f;
141
Geocent_ep2 = (1 / (1 - Geocent_e2)) - 1;
144
} /* END OF Set_Geocentric_Parameters */
147
void Get_Geocentric_Parameters (double *a,
149
{ /* BEGIN Get_Geocentric_Parameters */
151
* The function Get_Geocentric_Parameters returns the ellipsoid parameters
152
* to be used in geocentric coordinate conversions.
154
* a : Semi-major axis of ellipsoid, in meters. (output)
155
* f : Flattening of ellipsoid. (output)
160
} /* END OF Get_Geocentric_Parameters */
163
long Convert_Geodetic_To_Geocentric (double Latitude,
169
{ /* BEGIN Convert_Geodetic_To_Geocentric */
171
* The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
172
* (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
173
* according to the current ellipsoid parameters.
175
* Latitude : Geodetic latitude in radians (input)
176
* Longitude : Geodetic longitude in radians (input)
177
* Height : Geodetic height, in meters (input)
178
* X : Calculated Geocentric X coordinate, in meters (output)
179
* Y : Calculated Geocentric Y coordinate, in meters (output)
180
* Z : Calculated Geocentric Z coordinate, in meters (output)
183
long Error_Code = GEOCENT_NO_ERROR;
184
double Rn; /* Earth radius at location */
185
double Sin_Lat; /* sin(Latitude) */
186
double Sin2_Lat; /* Square of sin(Latitude) */
187
double Cos_Lat; /* cos(Latitude) */
189
if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
190
{ /* Latitude out of range */
191
Error_Code |= GEOCENT_LAT_ERROR;
193
if ((Longitude < -PI) || (Longitude > (2*PI)))
194
{ /* Longitude out of range */
195
Error_Code |= GEOCENT_LON_ERROR;
201
Sin_Lat = sin(Latitude);
202
Cos_Lat = cos(Latitude);
203
Sin2_Lat = Sin_Lat * Sin_Lat;
204
Rn = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * Sin2_Lat));
205
*X = (Rn + Height) * Cos_Lat * cos(Longitude);
206
*Y = (Rn + Height) * Cos_Lat * sin(Longitude);
207
*Z = ((Rn * (1 - Geocent_e2)) + Height) * Sin_Lat;
211
} /* END OF Convert_Geodetic_To_Geocentric */
214
void Convert_Geocentric_To_Geodetic (double X,
220
{ /* BEGIN Convert_Geocentric_To_Geodetic */
222
* The function Convert_Geocentric_To_Geodetic converts geocentric
223
* coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude,
224
* and height), according to the current ellipsoid parameters.
226
* X : Geocentric X coordinate, in meters. (input)
227
* Y : Geocentric Y coordinate, in meters. (input)
228
* Z : Geocentric Z coordinate, in meters. (input)
229
* Latitude : Calculated latitude value in radians. (output)
230
* Longitude : Calculated longitude value in radians. (output)
231
* Height : Calculated height value, in meters. (output)
233
* The method used here is derived from 'An Improved Algorithm for
234
* Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
237
/* Note: Variable names follow the notation used in Toms, Feb 1996 */
239
double W; /* distance from Z axis */
240
double W2; /* square of distance from Z axis */
241
double T0; /* initial estimate of vertical component */
242
double T1; /* corrected estimate of vertical component */
243
double S0; /* initial estimate of horizontal component */
244
double S1; /* corrected estimate of horizontal component */
245
double Sin_B0; /* sin(B0), B0 is estimate of Bowring aux variable */
246
double Sin3_B0; /* cube of sin(B0) */
247
double Cos_B0; /* cos(B0) */
248
double Sin_p1; /* sin(phi1), phi1 is estimated latitude */
249
double Cos_p1; /* cos(phi1) */
250
double Rn; /* Earth radius at location */
251
double Sum; /* numerator of cos(phi1) */
252
int At_Pole; /* indicates location is in polar region */
253
double Geocent_b = Geocent_a * (1 - Geocent_f); /* Semi-minor axis of ellipsoid, in meters */
258
*Longitude = atan2(Y,X);
264
*Longitude = PI_OVER_2;
268
*Longitude = -PI_OVER_2;
276
*Latitude = PI_OVER_2;
280
*Latitude = -PI_OVER_2;
283
{ /* center of earth */
284
*Latitude = PI_OVER_2;
285
*Height = -Geocent_b;
293
S0 = sqrt(T0 * T0 + W2);
296
Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
297
T1 = Z + Geocent_b * Geocent_ep2 * Sin3_B0;
298
Sum = W - Geocent_a * Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0;
299
S1 = sqrt(T1*T1 + Sum * Sum);
302
Rn = Geocent_a / sqrt(1.0 - Geocent_e2 * Sin_p1 * Sin_p1);
303
if (Cos_p1 >= COS_67P5)
305
*Height = W / Cos_p1 - Rn;
307
else if (Cos_p1 <= -COS_67P5)
309
*Height = W / -Cos_p1 - Rn;
313
*Height = Z / Sin_p1 + Rn * (Geocent_e2 - 1.0);
315
if (At_Pole == FALSE)
317
*Latitude = atan(Sin_p1 / Cos_p1);
319
} /* END OF Convert_Geocentric_To_Geodetic */
added
removed
dt_cc/geocent/geocent.c
