← Back to branch summary

~ubuntu-branches/ubuntu/natty/proj/natty

« back to all changes in this revision

Viewing changes to src/geocent.c

  • Committer: Bazaar Package Importer
  • Author(s): Peter S Galbraith
  • Date: 2007-12-21 21:16:56 UTC
  • mfrom: (1.2.4 upstream)
  • mto: This revision was merged to the branch mainline in revision 5.
  • Revision ID: james.westby@ubuntu.com-20071221211656-mtez2521i313j9mn
Tags: 4.6.0-1
New upstream release.

Show diffs side-by-side

added added

removed removed

Lines of Context:
src/geocent.c
65
65
 *    25-02-97          Original Code
66
66
 *
67
67
 * $Log: geocent.c,v $
 
68
 * Revision 1.7  2007/09/11 20:19:36  fwarmerdam
 
69
 * avoid use of static variables to make reentrant
 
70
 *
68
71
 * Revision 1.6  2006/01/12 22:29:01  fwarmerdam
69
72
 * make geocent.c globals static to avoid conflicts
70
73
 *
109
112
 
110
113
/***************************************************************************/
111
114
/*
112
 
 *                              GLOBAL DECLARATIONS
113
 
 */
114
 
/* Ellipsoid parameters, default to WGS 84 */
115
 
static double Geocent_a = 6378137.0;     /* Semi-major axis of ellipsoid in meters */
116
 
static double Geocent_b = 6356752.3142;  /* Semi-minor axis of ellipsoid           */
117
 
 
118
 
static double Geocent_a2 = 40680631590769.0;        /* Square of semi-major axis */
119
 
static double Geocent_b2 = 40408299984087.05;       /* Square of semi-minor axis */
120
 
static double Geocent_e2 = 0.0066943799901413800;   /* Eccentricity squared  */
121
 
static double Geocent_ep2 = 0.00673949675658690300; /* 2nd eccentricity squared */
122
 
/*
123
 
 * These state variables are for optimization purposes.  The only function
124
 
 * that should modify them is Set_Geocentric_Parameters.
125
 
 */
126
 
 
127
 
 
128
 
/***************************************************************************/
129
 
/*
130
115
 *                              FUNCTIONS     
131
116
 */
132
117
 
133
118
 
134
 
long pj_Set_Geocentric_Parameters (double a, double b) 
 
119
long pj_Set_Geocentric_Parameters (GeocentricInfo *gi, double a, double b) 
135
120
 
136
121
{ /* BEGIN Set_Geocentric_Parameters */
137
122
/*
141
126
 *    a  : Semi-major axis, in meters.          (input)
142
127
 *    b  : Semi-minor axis, in meters.          (input)
143
128
 */
144
 
  long Error_Code = GEOCENT_NO_ERROR;
 
129
    long Error_Code = GEOCENT_NO_ERROR;
145
130
 
146
 
  if (a <= 0.0)
147
 
    Error_Code |= GEOCENT_A_ERROR;
148
 
  if (b <= 0.0)
149
 
    Error_Code |= GEOCENT_B_ERROR;
150
 
  if (a < b)
151
 
    Error_Code |= GEOCENT_A_LESS_B_ERROR;
152
 
  if (!Error_Code)
153
 
  {
154
 
    Geocent_a = a;
155
 
    Geocent_b = b;
156
 
    Geocent_a2 = a * a;
157
 
    Geocent_b2 = b * b;
158
 
    Geocent_e2 = (Geocent_a2 - Geocent_b2) / Geocent_a2;
159
 
    Geocent_ep2 = (Geocent_a2 - Geocent_b2) / Geocent_b2;
160
 
  }
161
 
  return (Error_Code);
 
131
    if (a <= 0.0)
 
132
        Error_Code |= GEOCENT_A_ERROR;
 
133
    if (b <= 0.0)
 
134
        Error_Code |= GEOCENT_B_ERROR;
 
135
    if (a < b)
 
136
        Error_Code |= GEOCENT_A_LESS_B_ERROR;
 
137
    if (!Error_Code)
 
138
    {
 
139
        gi->Geocent_a = a;
 
140
        gi->Geocent_b = b;
 
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;
 
145
    }
 
146
    return (Error_Code);
162
147
} /* END OF Set_Geocentric_Parameters */
163
148
 
164
149
 
165
 
void pj_Get_Geocentric_Parameters (double *a, 
166
 
                                double *b)
 
150
void pj_Get_Geocentric_Parameters (GeocentricInfo *gi,
 
151
                                   double *a, 
 
152
                                   double *b)
167
153
{ /* BEGIN Get_Geocentric_Parameters */
168
154
/*
169
155
 * The function Get_Geocentric_Parameters returns the ellipsoid parameters
173
159
 *    b  : Semi-minor axis, in meters.          (output)
174
160
 */
175
161
 
176
 
  *a = Geocent_a;
177
 
  *b = Geocent_b;
 
162
    *a = gi->Geocent_a;
 
163
    *b = gi->Geocent_b;
178
164
} /* END OF Get_Geocentric_Parameters */
179
165
 
180
166
 
181
 
long pj_Convert_Geodetic_To_Geocentric (double Latitude,
182
 
                                     double Longitude,
183
 
                                     double Height,
184
 
                                     double *X,
185
 
                                     double *Y,
186
 
                                     double *Z) 
 
167
long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo *gi,
 
168
                                        double Latitude,
 
169
                                        double Longitude,
 
170
                                        double Height,
 
171
                                        double *X,
 
172
                                        double *Y,
 
173
                                        double *Z) 
187
174
{ /* BEGIN Convert_Geodetic_To_Geocentric */
188
175
/*
189
176
 * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
225
212
    Sin_Lat = sin(Latitude);
226
213
    Cos_Lat = cos(Latitude);
227
214
    Sin2_Lat = Sin_Lat * Sin_Lat;
228
 
    Rn = Geocent_a / (sqrt(1.0e0 - Geocent_e2 * Sin2_Lat));
 
215
    Rn = gi->Geocent_a / (sqrt(1.0e0 - gi->Geocent_e2 * Sin2_Lat));
229
216
    *X = (Rn + Height) * Cos_Lat * cos(Longitude);
230
217
    *Y = (Rn + Height) * Cos_Lat * sin(Longitude);
231
 
    *Z = ((Rn * (1 - Geocent_e2)) + Height) * Sin_Lat;
 
218
    *Z = ((Rn * (1 - gi->Geocent_e2)) + Height) * Sin_Lat;
232
219
 
233
220
  }
234
221
  return (Error_Code);
249
236
 
250
237
#define USE_ITERATIVE_METHOD
251
238
 
252
 
void pj_Convert_Geocentric_To_Geodetic (double X,
253
 
                                     double Y, 
254
 
                                     double Z,
255
 
                                     double *Latitude,
256
 
                                     double *Longitude,
257
 
                                     double *Height)
 
239
void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
 
240
                                        double X,
 
241
                                        double Y, 
 
242
                                        double Z,
 
243
                                        double *Latitude,
 
244
                                        double *Longitude,
 
245
                                        double *Height)
258
246
{ /* BEGIN Convert_Geocentric_To_Geodetic */
259
247
#if !defined(USE_ITERATIVE_METHOD)
260
248
/*
321
309
    Sin_B0 = T0 / S0;
322
310
    Cos_B0 = W / S0;
323
311
    Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
324
 
    T1 = Z + Geocent_b * Geocent_ep2 * Sin3_B0;
325
 
    Sum = W - Geocent_a * Geocent_e2 * Cos_B0 * Cos_B0 * Cos_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;
326
314
    S1 = sqrt(T1*T1 + Sum * Sum);
327
315
    Sin_p1 = T1 / S1;
328
316
    Cos_p1 = Sum / S1;
329
 
    Rn = Geocent_a / sqrt(1.0 - Geocent_e2 * Sin_p1 * Sin_p1);
 
317
    Rn = gi->Geocent_a / sqrt(1.0 - gi->Geocent_e2 * Sin_p1 * Sin_p1);
330
318
    if (Cos_p1 >= COS_67P5)
331
319
    {
332
320
        *Height = W / Cos_p1 - Rn;
337
325
    }
338
326
    else
339
327
    {
340
 
        *Height = Z / Sin_p1 + Rn * (Geocent_e2 - 1.0);
 
328
        *Height = Z / Sin_p1 + Rn * (gi->Geocent_e2 - 1.0);
341
329
    }
342
330
    if (At_Pole == FALSE)
343
331
    {
401
389
    RR = sqrt(X*X+Y*Y+Z*Z);
402
390
 
403
391
/*      special cases for latitude and longitude */
404
 
    if (P/Geocent_a < genau) {
 
392
    if (P/gi->Geocent_a < genau) {
405
393
 
406
394
/*  special case, if P=0. (X=0., Y=0.) */
407
395
        At_Pole = TRUE;
409
397
 
410
398
/*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
411
399
 *  of ellipsoid (=center of mass), Latitude becomes PI/2 */
412
 
        if (RR/Geocent_a < genau) {
 
400
        if (RR/gi->Geocent_a < genau) {
413
401
            *Latitude = PI_OVER_2;
414
 
            *Height   = -Geocent_b;
 
402
            *Height   = -gi->Geocent_b;
415
403
            return ;
416
404
 
417
405
        }
433
421
 */
434
422
    CT = Z/RR;
435
423
    ST = P/RR;
436
 
    RX = 1.0/sqrt(1.0-Geocent_e2*(2.0-Geocent_e2)*ST*ST);
437
 
    CPHI0 = ST*(1.0-Geocent_e2)*RX;
 
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;
438
426
    SPHI0 = CT*RX;
439
427
    iter = 0;
440
428
 
443
431
    do
444
432
    {
445
433
        iter++;
446
 
        RN = Geocent_a/sqrt(1.0-Geocent_e2*SPHI0*SPHI0);
 
434
        RN = gi->Geocent_a/sqrt(1.0-gi->Geocent_e2*SPHI0*SPHI0);
447
435
 
448
436
/*  ellipsoidal (geodetic) height */
449
 
        *Height = P*CPHI0+Z*SPHI0-RN*(1.0-Geocent_e2*SPHI0*SPHI0);
 
437
        *Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi->Geocent_e2*SPHI0*SPHI0);
450
438
 
451
 
        RK = Geocent_e2*RN/(RN+*Height);
 
439
        RK = gi->Geocent_e2*RN/(RN+*Height);
452
440
        RX = 1.0/sqrt(1.0-RK*(2.0-RK)*ST*ST);
453
441
        CPHI = ST*(1.0-RK)*RX;
454
442
        SPHI = CT*RX;