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 lessthan or equal to zero
23
* GEOCENT_B_ERROR : Semi-minor axis lessthan or equal to zero
24
* GEOCENT_A_LESS_B_ERROR : Semi-major axis less than semi-minor axis
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
68
* Revision 1.7 2007/09/11 20:19:36 fwarmerdam
69
* avoid use of static variables to make reentrant
71
* Revision 1.6 2006/01/12 22:29:01 fwarmerdam
72
* make geocent.c globals static to avoid conflicts
74
* Revision 1.5 2004/10/25 15:34:36 fwarmerdam
75
* make names of geodetic funcs from geotrans unique
77
* Revision 1.4 2004/05/03 16:28:01 warmerda
78
* Apply iterative solution to geocentric_to_geodetic as suggestion from
80
* http://bugzilla.remotesensing.org/show_bug.cgi?id=563
82
* Revision 1.3 2002/01/08 15:04:08 warmerda
83
* The latitude clamping fix from September in Convert_Geodetic_To_Geocentric
84
* was botched. Fixed up now.
89
/***************************************************************************/
96
* math.h - is needed for calls to sin, cos, tan and sqrt.
97
* geocent.h - is needed for Error codes and prototype error checking.
101
/***************************************************************************/
105
#define PI 3.14159265358979323e0
106
#define PI_OVER_2 (PI / 2.0e0)
109
#define COS_67P5 0.38268343236508977 /* cosine of 67.5 degrees */
110
#define AD_C 1.0026000 /* Toms region 1 constant */
113
/***************************************************************************/
119
long pj_Set_Geocentric_Parameters (GeocentricInfo *gi, double a, double b)
121
{ /* BEGIN Set_Geocentric_Parameters */
123
* The function Set_Geocentric_Parameters receives the ellipsoid parameters
124
* as inputs and sets the corresponding state variables.
126
* a : Semi-major axis, in meters. (input)
127
* b : Semi-minor axis, in meters. (input)
129
long Error_Code = GEOCENT_NO_ERROR;
132
Error_Code |= GEOCENT_A_ERROR;
134
Error_Code |= GEOCENT_B_ERROR;
136
Error_Code |= GEOCENT_A_LESS_B_ERROR;
141
gi->Geocent_a2 = a * a;
142
gi->Geocent_b2 = b * b;
143
gi->Geocent_e2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_a2;
144
gi->Geocent_ep2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_b2;
147
} /* END OF Set_Geocentric_Parameters */
150
void pj_Get_Geocentric_Parameters (GeocentricInfo *gi,
153
{ /* BEGIN Get_Geocentric_Parameters */
155
* The function Get_Geocentric_Parameters returns the ellipsoid parameters
156
* to be used in geocentric coordinate conversions.
158
* a : Semi-major axis, in meters. (output)
159
* b : Semi-minor axis, in meters. (output)
164
} /* END OF Get_Geocentric_Parameters */
167
long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo *gi,
174
{ /* BEGIN Convert_Geodetic_To_Geocentric */
176
* The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
177
* (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
178
* according to the current ellipsoid parameters.
180
* Latitude : Geodetic latitude in radians (input)
181
* Longitude : Geodetic longitude in radians (input)
182
* Height : Geodetic height, in meters (input)
183
* X : Calculated Geocentric X coordinate, in meters (output)
184
* Y : Calculated Geocentric Y coordinate, in meters (output)
185
* Z : Calculated Geocentric Z coordinate, in meters (output)
188
long Error_Code = GEOCENT_NO_ERROR;
189
double Rn; /* Earth radius at location */
190
double Sin_Lat; /* sin(Latitude) */
191
double Sin2_Lat; /* Square of sin(Latitude) */
192
double Cos_Lat; /* cos(Latitude) */
195
** Don't blow up if Latitude is just a little out of the value
196
** range as it may just be a rounding issue. Also removed longitude
197
** test, it should be wrapped by cos() and sin(). NFW for PROJ.4, Sep/2001.
199
if( Latitude < -PI_OVER_2 && Latitude > -1.001 * PI_OVER_2 )
200
Latitude = -PI_OVER_2;
201
else if( Latitude > PI_OVER_2 && Latitude < 1.001 * PI_OVER_2 )
202
Latitude = PI_OVER_2;
203
else if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
204
{ /* Latitude out of range */
205
Error_Code |= GEOCENT_LAT_ERROR;
212
Sin_Lat = sin(Latitude);
213
Cos_Lat = cos(Latitude);
214
Sin2_Lat = Sin_Lat * Sin_Lat;
215
Rn = gi->Geocent_a / (sqrt(1.0e0 - gi->Geocent_e2 * Sin2_Lat));
216
*X = (Rn + Height) * Cos_Lat * cos(Longitude);
217
*Y = (Rn + Height) * Cos_Lat * sin(Longitude);
218
*Z = ((Rn * (1 - gi->Geocent_e2)) + Height) * Sin_Lat;
222
} /* END OF Convert_Geodetic_To_Geocentric */
225
* The function Convert_Geocentric_To_Geodetic converts geocentric
226
* coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude,
227
* and height), according to the current ellipsoid parameters.
229
* X : Geocentric X coordinate, in meters. (input)
230
* Y : Geocentric Y coordinate, in meters. (input)
231
* Z : Geocentric Z coordinate, in meters. (input)
232
* Latitude : Calculated latitude value in radians. (output)
233
* Longitude : Calculated longitude value in radians. (output)
234
* Height : Calculated height value, in meters. (output)
237
#define USE_ITERATIVE_METHOD
239
void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
246
{ /* BEGIN Convert_Geocentric_To_Geodetic */
247
#if !defined(USE_ITERATIVE_METHOD)
249
* The method used here is derived from 'An Improved Algorithm for
250
* Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
253
/* Note: Variable names follow the notation used in Toms, Feb 1996 */
255
double W; /* distance from Z axis */
256
double W2; /* square of distance from Z axis */
257
double T0; /* initial estimate of vertical component */
258
double T1; /* corrected estimate of vertical component */
259
double S0; /* initial estimate of horizontal component */
260
double S1; /* corrected estimate of horizontal component */
261
double Sin_B0; /* sin(B0), B0 is estimate of Bowring aux variable */
262
double Sin3_B0; /* cube of sin(B0) */
263
double Cos_B0; /* cos(B0) */
264
double Sin_p1; /* sin(phi1), phi1 is estimated latitude */
265
double Cos_p1; /* cos(phi1) */
266
double Rn; /* Earth radius at location */
267
double Sum; /* numerator of cos(phi1) */
268
int At_Pole; /* indicates location is in polar region */
273
*Longitude = atan2(Y,X);
279
*Longitude = PI_OVER_2;
283
*Longitude = -PI_OVER_2;
291
*Latitude = PI_OVER_2;
295
*Latitude = -PI_OVER_2;
298
{ /* center of earth */
299
*Latitude = PI_OVER_2;
300
*Height = -Geocent_b;
308
S0 = sqrt(T0 * T0 + W2);
311
Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
312
T1 = Z + gi->Geocent_b * gi->Geocent_ep2 * Sin3_B0;
313
Sum = W - gi->Geocent_a * gi->Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0;
314
S1 = sqrt(T1*T1 + Sum * Sum);
317
Rn = gi->Geocent_a / sqrt(1.0 - gi->Geocent_e2 * Sin_p1 * Sin_p1);
318
if (Cos_p1 >= COS_67P5)
320
*Height = W / Cos_p1 - Rn;
322
else if (Cos_p1 <= -COS_67P5)
324
*Height = W / -Cos_p1 - Rn;
328
*Height = Z / Sin_p1 + Rn * (gi->Geocent_e2 - 1.0);
330
if (At_Pole == FALSE)
332
*Latitude = atan(Sin_p1 / Cos_p1);
334
#else /* defined(USE_ITERATIVE_METHOD) */
338
* Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für
339
* das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover
340
* Nr. 137, p. 130-131.
342
* Programmed by GGA- Leibniz-Institue of Applied Geophysics
345
* Federal Republic of Germany
346
* Internet: www.gga-hannover.de
348
* Hannover, March 1999, April 2004.
349
* see also: comments in statements
351
* Mathematically exact and because of symmetry of rotation-ellipsoid,
352
* each point (X,Y,Z) has at least two solutions (Latitude1,Longitude1,Height1) and
353
* (Latitude2,Longitude2,Height2). Is point=(0.,0.,Z) (P=0.), so you get even
354
* four solutions, every two symmetrical to the semi-minor axis.
355
* Here Height1 and Height2 have at least a difference in order of
356
* radius of curvature (e.g. (0,0,b)=> (90.,0.,0.) or (-90.,0.,-2b);
357
* (a+100.)*(sqrt(2.)/2.,sqrt(2.)/2.,0.) => (0.,45.,100.) or
358
* (0.,225.,-(2a+100.))).
359
* The algorithm always computes (Latitude,Longitude) with smallest |Height|.
360
* For normal computations, that means |Height|<10000.m, algorithm normally
361
* converges after to 2-3 steps!!!
362
* But if |Height| has the amount of length of ellipsoid's axis
363
* (e.g. -6300000.m), algorithm needs about 15 steps.
366
/* local defintions and variables */
367
/* end-criterium of loop, accuracy of sin(Latitude) */
369
#define genau2 (genau*genau)
372
double P; /* distance between semi-minor axis and location */
373
double RR; /* distance between center and location */
374
double CT; /* sin of geocentric latitude */
375
double ST; /* cos of geocentric latitude */
378
double RN; /* Earth radius at location */
379
double CPHI0; /* cos of start or old geodetic latitude in iterations */
380
double SPHI0; /* sin of start or old geodetic latitude in iterations */
381
double CPHI; /* cos of searched geodetic latitude */
382
double SPHI; /* sin of searched geodetic latitude */
383
double SDPHI; /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
384
int At_Pole; /* indicates location is in polar region */
385
int iter; /* # of continous iteration, max. 30 is always enough (s.a.) */
389
RR = sqrt(X*X+Y*Y+Z*Z);
391
/* special cases for latitude and longitude */
392
if (P/gi->Geocent_a < genau) {
394
/* special case, if P=0. (X=0., Y=0.) */
398
/* if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
399
* of ellipsoid (=center of mass), Latitude becomes PI/2 */
400
if (RR/gi->Geocent_a < genau) {
401
*Latitude = PI_OVER_2;
402
*Height = -gi->Geocent_b;
408
/* ellipsoidal (geodetic) longitude
409
* interval: -PI < Longitude <= +PI */
410
*Longitude=atan2(Y,X);
413
/* --------------------------------------------------------------
414
* Following iterative algorithm was developped by
415
* "Institut für Erdmessung", University of Hannover, July 1988.
416
* Internet: www.ife.uni-hannover.de
417
* Iterative computation of CPHI,SPHI and Height.
418
* Iteration of CPHI and SPHI to 10**-12 radian resp.
420
* --------------------------------------------------------------
424
RX = 1.0/sqrt(1.0-gi->Geocent_e2*(2.0-gi->Geocent_e2)*ST*ST);
425
CPHI0 = ST*(1.0-gi->Geocent_e2)*RX;
429
/* loop to find sin(Latitude) resp. Latitude
430
* until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
434
RN = gi->Geocent_a/sqrt(1.0-gi->Geocent_e2*SPHI0*SPHI0);
436
/* ellipsoidal (geodetic) height */
437
*Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi->Geocent_e2*SPHI0*SPHI0);
439
RK = gi->Geocent_e2*RN/(RN+*Height);
440
RX = 1.0/sqrt(1.0-RK*(2.0-RK)*ST*ST);
441
CPHI = ST*(1.0-RK)*RX;
443
SDPHI = SPHI*CPHI0-CPHI*SPHI0;
447
while (SDPHI*SDPHI > genau2 && iter < maxiter);
449
/* ellipsoidal (geodetic) latitude */
450
*Latitude=atan(SPHI/fabs(CPHI));
453
#endif /* defined(USE_ITERATIVE_METHOD) */
454
} /* END OF Convert_Geocentric_To_Geodetic */
added
removed
include/ggl/projections/impl/geocent.c
