~ubuntu-branches/ubuntu/raring/geotranz/raring

« back to all changes in this revision

Viewing changes to dt_cc/eckert4/eckert4.c

  • Committer: Bazaar Package Importer
  • Author(s): Roberto Lumbreras
  • Date: 2008-10-17 14:43:09 UTC
  • Revision ID: james.westby@ubuntu.com-20081017144309-jb7uzfi1y1lvez8j
Tags: upstream-2.4.2
ImportĀ upstreamĀ versionĀ 2.4.2

Show diffs side-by-side

added added

removed removed

Lines of Context:
 
1
/***************************************************************************/
 
2
/* RSC IDENTIFIER: ECKERT4
 
3
 *
 
4
 * ABSTRACT
 
5
 *
 
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.
 
11
 *
 
12
 * ERROR HANDLING
 
13
 *
 
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:
 
18
 *
 
19
 *          ECK4_NO_ERROR           : No errors occurred in function
 
20
 *          ECK4_LAT_ERROR          : Latitude outside of valid range
 
21
 *                                      (-90 to 90 degrees)
 
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
 
34
 *                                                                                          (250 to 350)
 
35
 *
 
36
 * REUSE NOTES
 
37
 *
 
38
 *    ECKERT4 is intended for reuse by any application that performs a
 
39
 *    Eckert IV projection or its inverse.
 
40
 *
 
41
 * REFERENCES
 
42
 *
 
43
 *    Further information on ECKERT4 can be found in the Reuse Manual.
 
44
 *
 
45
 *    ECKERT4 originated from :  U.S. Army Topographic Engineering Center
 
46
 *                                Geospatial Information Division
 
47
 *                                7701 Telegraph Road
 
48
 *                                Alexandria, VA  22310-3864
 
49
 *
 
50
 * LICENSES
 
51
 *
 
52
 *    None apply to this component.
 
53
 *
 
54
 * RESTRICTIONS
 
55
 *
 
56
 *    ECKERT4 has no restrictions.
 
57
 *
 
58
 * ENVIRONMENT
 
59
 *
 
60
 *    ECKERT4 was tested and certified in the following environments:
 
61
 *
 
62
 *    1. Solaris 2.5 with GCC 2.8.1
 
63
 *    2. MS Windows 95 with MS Visual C++ 6
 
64
 *
 
65
 * MODIFICATIONS
 
66
 *
 
67
 *    Date              Description
 
68
 *    ----              -----------
 
69
 *    04/16/99          Original Code
 
70
 *
 
71
 */
 
72
 
 
73
 
 
74
/***************************************************************************/
 
75
/*
 
76
 *                               INCLUDES
 
77
 */
 
78
 
 
79
#include <math.h>
 
80
#include "eckert4.h"
 
81
 
 
82
/*
 
83
 *    math.h    - Standard C math library
 
84
 *    eckert4.h - Is for prototype error checking.
 
85
 */
 
86
 
 
87
 
 
88
/***************************************************************************/
 
89
/*
 
90
 *                               DEFINES
 
91
 */
 
92
 
 
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)
 
97
 
 
98
 
 
99
/***************************************************************************/
 
100
/*
 
101
 *                               GLOBALS
 
102
 */
 
103
 
 
104
const double two_PLUS_PI_OVER_2 = (2.0 + PI / 2.0);
 
105
 
 
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  */
 
112
 
 
113
static double Ra0 = 2690082.6043273;                   /* 0.4222382 * Sperical Radius (6371007.1810824) */
 
114
static double Ra1 = 8451143.5741087;                   /* 1.3265004 * Sperical Radius (6371007.1810824) */
 
115
 
 
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;
 
123
/*
 
124
 * These state variables are for optimization purposes.  The only function
 
125
 * that should modify them is Set_Eckert IV_Parameters.
 
126
 */
 
127
 
 
128
 
 
129
/***************************************************************************/
 
130
/*
 
131
 *                              FUNCTIONS
 
132
 */
 
133
 
 
134
 
 
135
long Set_Eckert4_Parameters(double a,
 
136
                            double f,                
 
137
                            double Central_Meridian,
 
138
                            double False_Easting,
 
139
                            double False_Northing)
 
140
 
 
141
{ /* Begin Set_Eckert4_Parameters */
 
142
/*
 
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.
 
147
 *
 
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)
 
151
 *                          the projection
 
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)
 
156
 */
 
157
 
 
158
  double Ra;              /* Spherical radius */
 
159
  double inv_f = 1 / f;
 
160
  long Error_Code = ECK4_NO_ERROR;
 
161
 
 
162
  if (a <= 0.0)
 
163
  { /* Semi-major axis must be greater than zero */
 
164
    Error_Code |= ECK4_A_ERROR;
 
165
  }
 
166
  if ((inv_f < 250) || (inv_f > 350))
 
167
  { /* Inverse flattening must be between 250 and 350 */
 
168
    Error_Code |= ECK4_INV_F_ERROR;
 
169
  }
 
170
  if ((Central_Meridian < -PI) || (Central_Meridian > TWO_PI))
 
171
  { /* origin longitude out of range */
 
172
    Error_Code |= ECK4_CENT_MER_ERROR;
 
173
  }
 
174
  if (!Error_Code)
 
175
  { /* no errors */
 
176
    Eck4_a = a;
 
177
    Eck4_f = f;
 
178
    es2 = 2 * Eck4_f - Eck4_f * Eck4_f;
 
179
    es4 = es2 * es2;
 
180
    es6 = es4 * es2;
 
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)
 
191
    {
 
192
      Eck4_Max_Easting = 16808386.0;
 
193
      Eck4_Min_Easting = -16902288.0;
 
194
    }
 
195
    else if (Eck4_Origin_Long < 0)
 
196
    {
 
197
      Eck4_Max_Easting = 16902288.0;
 
198
      Eck4_Min_Easting = -16808386.0;
 
199
    }
 
200
    else
 
201
    {
 
202
      Eck4_Max_Easting = 16902288.0;
 
203
      Eck4_Min_Easting = -16902288.0;
 
204
    }
 
205
 
 
206
  } /* End if(!Error_Code) */
 
207
  return (Error_Code);
 
208
} /* End Set_Eckert4_Parameters */
 
209
 
 
210
 
 
211
void Get_Eckert4_Parameters(double *a,
 
212
                            double *f,                           
 
213
                            double *Central_Meridian,
 
214
                            double *False_Easting,
 
215
                            double *False_Northing)
 
216
{ /* Begin Get_Eckert4_Parameters */
 
217
/*
 
218
 * The function Get_Eckert4_Parameters returns the current ellipsoid
 
219
 * parameters and Eckert IV projection parameters.
 
220
 *
 
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)
 
224
 *                          the projection
 
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)
 
229
 */
 
230
 
 
231
  *a = Eck4_a;
 
232
  *f = Eck4_f;
 
233
  *Central_Meridian = Eck4_Origin_Long;
 
234
  *False_Easting = Eck4_False_Easting;
 
235
  *False_Northing = Eck4_False_Northing;
 
236
  return;
 
237
} /* End Get_Eckert4_Parameters */
 
238
 
 
239
 
 
240
long Convert_Geodetic_To_Eckert4 (double Latitude,
 
241
                                  double Longitude,
 
242
                                  double *Easting,
 
243
                                  double *Northing)
 
244
 
 
245
{ /* Begin Convert_Geodetic_To_Eckert4 */
 
246
/*
 
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.
 
253
 *
 
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)
 
258
 */
 
259
 
 
260
  double slat = sin(Latitude);
 
261
  double sin_theta, cos_theta;
 
262
  double num;
 
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;
 
269
 
 
270
  if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
 
271
  {  /* Latitude out of range */
 
272
    Error_Code |= ECK4_LAT_ERROR;
 
273
  }
 
274
  if ((Longitude < -PI) || (Longitude > TWO_PI))
 
275
  {  /* Longitude out of range */
 
276
    Error_Code|= ECK4_LON_ERROR;
 
277
  }
 
278
 
 
279
  if (!Error_Code)
 
280
  { /* no errors */
 
281
 
 
282
    dlam = Longitude - Eck4_Origin_Long;
 
283
    if (dlam > PI)
 
284
    {
 
285
      dlam -= TWO_PI;
 
286
    }
 
287
    if (dlam < -PI)
 
288
    {
 
289
      dlam += TWO_PI;
 
290
    }
 
291
    while (fabs(delta_theta) > dt_tolerance)
 
292
    {
 
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;
 
299
    }
 
300
    *Easting = Ra0 * dlam * (1.0 + cos(theta)) + Eck4_False_Easting;
 
301
    *Northing = Ra1 * sin(theta) + Eck4_False_Northing;
 
302
 
 
303
  }
 
304
  return (Error_Code);
 
305
 
 
306
} /* End Convert_Geodetic_To_Eckert4 */
 
307
 
 
308
 
 
309
long Convert_Eckert4_To_Geodetic(double Easting,
 
310
                                 double Northing,
 
311
                                 double *Latitude,
 
312
                                 double *Longitude)
 
313
{ /* Begin Convert_Eckert4_To_Geodetic */
 
314
/*
 
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.
 
321
 *
 
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)
 
326
 */
 
327
 
 
328
  double theta;
 
329
  double sin_theta, cos_theta;
 
330
  double num;
 
331
  double dx, dy;
 
332
  double i;
 
333
 
 
334
  long Error_Code = ECK4_NO_ERROR;
 
335
 
 
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;
 
340
  }
 
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;
 
345
  }
 
346
 
 
347
  if (!Error_Code)
 
348
  {
 
349
    dy = Northing - Eck4_False_Northing;
 
350
    dx = Easting - Eck4_False_Easting;
 
351
    i = dy / Ra1;
 
352
    if (i > 1.0)
 
353
      i = 1.0;
 
354
    else if (i < -1.0)
 
355
      i = -1.0;
 
356
 
 
357
    theta = asin(i);
 
358
    sin_theta = sin(theta);
 
359
    cos_theta = cos(theta);
 
360
    num = NUM(theta, sin_theta, cos_theta);
 
361
 
 
362
    *Latitude = asin(num / two_PLUS_PI_OVER_2);
 
363
    *Longitude = Eck4_Origin_Long + dx / (Ra0 * (1 + cos_theta));
 
364
 
 
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;
 
369
 
 
370
    if (*Longitude > PI)
 
371
      *Longitude -= TWO_PI;
 
372
    if (*Longitude < -PI)
 
373
      *Longitude += TWO_PI;
 
374
 
 
375
    if (*Longitude > PI)  /* force distorted values to 180, -180 degrees */
 
376
      *Longitude = PI;
 
377
    else if (*Longitude < -PI)
 
378
      *Longitude = -PI;
 
379
 
 
380
  }
 
381
  return (Error_Code);
 
382
 
 
383
} /* End Convert_Eckert4_To_Geodetic */
 
384