← Back to branch summary

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

« back to all changes in this revision

Viewing changes to dt_cc/tranmerc/tranmerc.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:
dt_cc/tranmerc/tranmerc.c
 
1
/***************************************************************************/
 
2
/* RSC IDENTIFIER: TRANSVERSE MERCATOR
 
3
 *
 
4
 * ABSTRACT
 
5
 *
 
6
 *    This component provides conversions between Geodetic coordinates 
 
7
 *    (latitude and longitude) and Transverse Mercator projection coordinates
 
8
 *    (easting and northing).
 
9
 *
 
10
 * ERROR HANDLING
 
11
 *
 
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:
 
16
 *
 
17
 *       TRANMERC_NO_ERROR           : No errors occurred in function
 
18
 *       TRANMERC_LAT_ERROR          : Latitude outside of valid range
 
19
 *                                      (-90 to 90 degrees)
 
20
 *       TRANMERC_LON_ERROR          : Longitude outside of valid range
 
21
 *                                      (-180 to 360 degrees, and within
 
22
 *                                        +/-90 of Central Meridian)
 
23
 *       TRANMERC_EASTING_ERROR      : Easting outside of valid range
 
24
 *                                      (depending on ellipsoid and
 
25
 *                                       projection parameters)
 
26
 *       TRANMERC_NORTHING_ERROR     : Northing outside of valid range
 
27
 *                                      (depending on ellipsoid and
 
28
 *                                       projection parameters)
 
29
 *       TRANMERC_ORIGIN_LAT_ERROR   : Origin latitude outside of valid range
 
30
 *                                      (-90 to 90 degrees)
 
31
 *       TRANMERC_CENT_MER_ERROR     : Central meridian outside of valid range
 
32
 *                                      (-180 to 360 degrees)
 
33
 *       TRANMERC_A_ERROR            : Semi-major axis less than or equal to zero
 
34
 *       TRANMERC_INV_F_ERROR        : Inverse flattening outside of valid range
 
35
 *                                                                                        (250 to 350)
 
36
 *       TRANMERC_SCALE_FACTOR_ERROR : Scale factor outside of valid
 
37
 *                                     range (0.3 to 3.0)
 
38
 *               TM_LON_WARNING              : Distortion will result if longitude is more
 
39
 *                                       than 9 degrees from the Central Meridian
 
40
 *
 
41
 * REUSE NOTES
 
42
 *
 
43
 *    TRANSVERSE MERCATOR is intended for reuse by any application that 
 
44
 *    performs a Transverse Mercator projection or its inverse.
 
45
 *    
 
46
 * REFERENCES
 
47
 *
 
48
 *    Further information on TRANSVERSE MERCATOR can be found in the 
 
49
 *    Reuse Manual.
 
50
 *
 
51
 *    TRANSVERSE MERCATOR originated from :  
 
52
 *                      U.S. Army Topographic Engineering Center
 
53
 *                      Geospatial Information Division
 
54
 *                      7701 Telegraph Road
 
55
 *                      Alexandria, VA  22310-3864
 
56
 *
 
57
 * LICENSES
 
58
 *
 
59
 *    None apply to this component.
 
60
 *
 
61
 * RESTRICTIONS
 
62
 *
 
63
 *    TRANSVERSE MERCATOR has no restrictions.
 
64
 *
 
65
 * ENVIRONMENT
 
66
 *
 
67
 *    TRANSVERSE MERCATOR was tested and certified in the following 
 
68
 *    environments:
 
69
 *
 
70
 *    1. Solaris 2.5 with GCC, version 2.8.1
 
71
 *    2. Windows 95 with MS Visual C++, version 6
 
72
 *
 
73
 * MODIFICATIONS
 
74
 *
 
75
 *    Date              Description
 
76
 *    ----              -----------
 
77
 *    10-02-97          Original Code
 
78
 *    03-02-97          Re-engineered Code
 
79
 *
 
80
 */
 
81
 
 
82
 
 
83
/***************************************************************************/
 
84
/*
 
85
 *                               INCLUDES
 
86
 */
 
87
 
 
88
#include <math.h>
 
89
#include "tranmerc.h"
 
90
 
 
91
/*
 
92
 *    math.h      - Standard C math library
 
93
 *    tranmerc.h  - Is for prototype error checking
 
94
 */
 
95
 
 
96
 
 
97
/***************************************************************************/
 
98
/*                               DEFINES 
 
99
 *
 
100
 */
 
101
 
 
102
#define PI              3.14159265358979323e0   /* PI     */
 
103
#define PI_OVER_2         (PI/2.0e0)            /* PI over 2 */
 
104
#define MAX_LAT         ((PI * 89.99)/180.0)    /* 89.99 degrees in radians */
 
105
#define MAX_DELTA_LONG  ((PI * 90)/180.0)       /* 90 degrees in radians */
 
106
#define MIN_SCALE_FACTOR  0.3
 
107
#define MAX_SCALE_FACTOR  3.0
 
108
 
 
109
#define SPHTMD(Latitude) ((double) (TranMerc_ap * Latitude \
 
110
      - TranMerc_bp * sin(2.e0 * Latitude) + TranMerc_cp * sin(4.e0 * Latitude) \
 
111
      - TranMerc_dp * sin(6.e0 * Latitude) + TranMerc_ep * sin(8.e0 * Latitude) ) )
 
112
 
 
113
#define SPHSN(Latitude) ((double) (TranMerc_a / sqrt( 1.e0 - TranMerc_es * \
 
114
      pow(sin(Latitude), 2))))
 
115
 
 
116
#define SPHSR(Latitude) ((double) (TranMerc_a * (1.e0 - TranMerc_es) / \
 
117
    pow(DENOM(Latitude), 3)))
 
118
 
 
119
#define DENOM(Latitude) ((double) (sqrt(1.e0 - TranMerc_es * pow(sin(Latitude),2))))
 
120
 
 
121
 
 
122
/**************************************************************************/
 
123
/*                               GLOBAL DECLARATIONS
 
124
 *
 
125
 */
 
126
 
 
127
/* Ellipsoid Parameters, default to WGS 84  */
 
128
static double TranMerc_a = 6378137.0;              /* Semi-major axis of ellipsoid in meters */
 
129
static double TranMerc_f = 1 / 298.257223563;      /* Flattening of ellipsoid  */
 
130
static double TranMerc_es = 0.0066943799901413800; /* Eccentricity (0.08181919084262188000) squared */
 
131
static double TranMerc_ebs = 0.0067394967565869;   /* Second Eccentricity squared */
 
132
 
 
133
/* Transverse_Mercator projection Parameters */
 
134
static double TranMerc_Origin_Lat = 0.0;           /* Latitude of origin in radians */
 
135
static double TranMerc_Origin_Long = 0.0;          /* Longitude of origin in radians */
 
136
static double TranMerc_False_Northing = 0.0;       /* False northing in meters */
 
137
static double TranMerc_False_Easting = 0.0;        /* False easting in meters */
 
138
static double TranMerc_Scale_Factor = 1.0;         /* Scale factor  */
 
139
 
 
140
/* Isometeric to geodetic latitude parameters, default to WGS 84 */
 
141
static double TranMerc_ap = 6367449.1458008;
 
142
static double TranMerc_bp = 16038.508696861;
 
143
static double TranMerc_cp = 16.832613334334;
 
144
static double TranMerc_dp = 0.021984404273757;
 
145
static double TranMerc_ep = 3.1148371319283e-005;
 
146
 
 
147
/* Maximum variance for easting and northing values for WGS 84. */
 
148
static double TranMerc_Delta_Easting = 40000000.0;
 
149
static double TranMerc_Delta_Northing = 40000000.0;
 
150
 
 
151
/* These state variables are for optimization purposes. The only function
 
152
 * that should modify them is Set_Tranverse_Mercator_Parameters.         */
 
153
 
 
154
 
 
155
/************************************************************************/
 
156
/*                              FUNCTIONS     
 
157
 *
 
158
 */
 
159
 
 
160
 
 
161
long Set_Transverse_Mercator_Parameters(double a,
 
162
                                        double f,
 
163
                                        double Origin_Latitude,
 
164
                                        double Central_Meridian,
 
165
                                        double False_Easting,
 
166
                                        double False_Northing,
 
167
                                        double Scale_Factor)
 
168
 
 
169
{ /* BEGIN Set_Tranverse_Mercator_Parameters */
 
170
  /*
 
171
   * The function Set_Tranverse_Mercator_Parameters receives the ellipsoid
 
172
   * parameters and Tranverse Mercator projection parameters as inputs, and
 
173
   * sets the corresponding state variables. If any errors occur, the error
 
174
   * code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
 
175
   * returned.
 
176
   *
 
177
   *    a                 : Semi-major axis of ellipsoid, in meters    (input)
 
178
   *    f                 : Flattening of ellipsoid                                                      (input)
 
179
   *    Origin_Latitude   : Latitude in radians at the origin of the   (input)
 
180
   *                         projection
 
181
   *    Central_Meridian  : Longitude in radians at the center of the  (input)
 
182
   *                         projection
 
183
   *    False_Easting     : Easting/X at the center of the projection  (input)
 
184
   *    False_Northing    : Northing/Y at the center of the projection (input)
 
185
   *    Scale_Factor      : Projection scale factor                    (input) 
 
186
   */
 
187
 
 
188
  double tn;        /* True Meridianal distance constant  */
 
189
  double tn2;
 
190
  double tn3;
 
191
  double tn4;
 
192
  double tn5;
 
193
  double dummy_northing;
 
194
  double TranMerc_b; /* Semi-minor axis of ellipsoid, in meters */
 
195
  double inv_f = 1 / f;
 
196
  long Error_Code = TRANMERC_NO_ERROR;
 
197
 
 
198
  if (a <= 0.0)
 
199
  { /* Semi-major axis must be greater than zero */
 
200
    Error_Code |= TRANMERC_A_ERROR;
 
201
  }
 
202
  if ((inv_f < 250) || (inv_f > 350))
 
203
  { /* Inverse flattening must be between 250 and 350 */
 
204
    Error_Code |= TRANMERC_INV_F_ERROR;
 
205
  }
 
206
  if ((Origin_Latitude < -PI_OVER_2) || (Origin_Latitude > PI_OVER_2))
 
207
  { /* origin latitude out of range */
 
208
    Error_Code |= TRANMERC_ORIGIN_LAT_ERROR;
 
209
  }
 
210
  if ((Central_Meridian < -PI) || (Central_Meridian > (2*PI)))
 
211
  { /* origin longitude out of range */
 
212
    Error_Code |= TRANMERC_CENT_MER_ERROR;
 
213
  }
 
214
  if ((Scale_Factor < MIN_SCALE_FACTOR) || (Scale_Factor > MAX_SCALE_FACTOR))
 
215
  {
 
216
    Error_Code |= TRANMERC_SCALE_FACTOR_ERROR;
 
217
  }
 
218
  if (!Error_Code)
 
219
  { /* no errors */
 
220
    TranMerc_a = a;
 
221
    TranMerc_f = f;
 
222
    TranMerc_Origin_Lat = Origin_Latitude;
 
223
    if (Central_Meridian > PI)
 
224
      Central_Meridian -= (2*PI);
 
225
    TranMerc_Origin_Long = Central_Meridian;
 
226
    TranMerc_False_Northing = False_Northing;
 
227
    TranMerc_False_Easting = False_Easting; 
 
228
    TranMerc_Scale_Factor = Scale_Factor;
 
229
 
 
230
    /* Eccentricity Squared */
 
231
    TranMerc_es = 2 * TranMerc_f - TranMerc_f * TranMerc_f;
 
232
    /* Second Eccentricity Squared */
 
233
    TranMerc_ebs = (1 / (1 - TranMerc_es)) - 1;
 
234
 
 
235
    TranMerc_b = TranMerc_a * (1 - TranMerc_f);    
 
236
    /*True meridianal constants  */
 
237
    tn = (TranMerc_a - TranMerc_b) / (TranMerc_a + TranMerc_b);
 
238
    tn2 = tn * tn;
 
239
    tn3 = tn2 * tn;
 
240
    tn4 = tn3 * tn;
 
241
    tn5 = tn4 * tn;
 
242
 
 
243
    TranMerc_ap = TranMerc_a * (1.e0 - tn + 5.e0 * (tn2 - tn3)/4.e0
 
244
                                + 81.e0 * (tn4 - tn5)/64.e0 );
 
245
    TranMerc_bp = 3.e0 * TranMerc_a * (tn - tn2 + 7.e0 * (tn3 - tn4)
 
246
                                       /8.e0 + 55.e0 * tn5/64.e0 )/2.e0;
 
247
    TranMerc_cp = 15.e0 * TranMerc_a * (tn2 - tn3 + 3.e0 * (tn4 - tn5 )/4.e0) /16.0;
 
248
    TranMerc_dp = 35.e0 * TranMerc_a * (tn3 - tn4 + 11.e0 * tn5 / 16.e0) / 48.e0;
 
249
    TranMerc_ep = 315.e0 * TranMerc_a * (tn4 - tn5) / 512.e0;
 
250
    Convert_Geodetic_To_Transverse_Mercator(MAX_LAT,
 
251
                                            MAX_DELTA_LONG + Central_Meridian,
 
252
                                            &TranMerc_Delta_Easting,
 
253
                                            &TranMerc_Delta_Northing);
 
254
    Convert_Geodetic_To_Transverse_Mercator(0,
 
255
                                            MAX_DELTA_LONG + Central_Meridian,
 
256
                                            &TranMerc_Delta_Easting,
 
257
                                            &dummy_northing);
 
258
    TranMerc_Delta_Northing++;
 
259
    TranMerc_Delta_Easting++;
 
260
 
 
261
  } /* END OF if(!Error_Code) */
 
262
  return (Error_Code);
 
263
}  /* END of Set_Transverse_Mercator_Parameters  */
 
264
 
 
265
 
 
266
void Get_Transverse_Mercator_Parameters(double *a,
 
267
                                        double *f,
 
268
                                        double *Origin_Latitude,
 
269
                                        double *Central_Meridian,
 
270
                                        double *False_Easting,
 
271
                                        double *False_Northing,
 
272
                                        double *Scale_Factor)
 
273
 
 
274
{ /* BEGIN Get_Tranverse_Mercator_Parameters  */
 
275
  /*
 
276
   * The function Get_Transverse_Mercator_Parameters returns the current
 
277
   * ellipsoid and Transverse Mercator projection parameters.
 
278
   *
 
279
   *    a                 : Semi-major axis of ellipsoid, in meters    (output)
 
280
   *    f                 : Flattening of ellipsoid                                                      (output)
 
281
   *    Origin_Latitude   : Latitude in radians at the origin of the   (output)
 
282
   *                         projection
 
283
   *    Central_Meridian  : Longitude in radians at the center of the  (output)
 
284
   *                         projection
 
285
   *    False_Easting     : Easting/X at the center of the projection  (output)
 
286
   *    False_Northing    : Northing/Y at the center of the projection (output)
 
287
   *    Scale_Factor      : Projection scale factor                    (output) 
 
288
   */
 
289
 
 
290
  *a = TranMerc_a;
 
291
  *f = TranMerc_f;
 
292
  *Origin_Latitude = TranMerc_Origin_Lat;
 
293
  *Central_Meridian = TranMerc_Origin_Long;
 
294
  *False_Easting = TranMerc_False_Easting;
 
295
  *False_Northing = TranMerc_False_Northing;
 
296
  *Scale_Factor = TranMerc_Scale_Factor;
 
297
  return;
 
298
} /* END OF Get_Tranverse_Mercator_Parameters */
 
299
 
 
300
 
 
301
 
 
302
long Convert_Geodetic_To_Transverse_Mercator (double Latitude,
 
303
                                              double Longitude,
 
304
                                              double *Easting,
 
305
                                              double *Northing)
 
306
 
 
307
{      /* BEGIN Convert_Geodetic_To_Transverse_Mercator */
 
308
 
 
309
  /*
 
310
   * The function Convert_Geodetic_To_Transverse_Mercator converts geodetic
 
311
   * (latitude and longitude) coordinates to Transverse Mercator projection
 
312
   * (easting and northing) coordinates, according to the current ellipsoid
 
313
   * and Transverse Mercator projection coordinates.  If any errors occur, the
 
314
   * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
 
315
   * returned.
 
316
   *
 
317
   *    Latitude      : Latitude in radians                         (input)
 
318
   *    Longitude     : Longitude in radians                        (input)
 
319
   *    Easting       : Easting/X in meters                         (output)
 
320
   *    Northing      : Northing/Y in meters                        (output)
 
321
   */
 
322
 
 
323
  double c;       /* Cosine of latitude                          */
 
324
  double c2;
 
325
  double c3;
 
326
  double c5;
 
327
  double c7;
 
328
  double dlam;    /* Delta longitude - Difference in Longitude       */
 
329
  double eta;     /* constant - TranMerc_ebs *c *c                   */
 
330
  double eta2;
 
331
  double eta3;
 
332
  double eta4;
 
333
  double s;       /* Sine of latitude                        */
 
334
  double sn;      /* Radius of curvature in the prime vertical       */
 
335
  double t;       /* Tangent of latitude                             */
 
336
  double tan2;
 
337
  double tan3;
 
338
  double tan4;
 
339
  double tan5;
 
340
  double tan6;
 
341
  double t1;      /* Term in coordinate conversion formula - GP to Y */
 
342
  double t2;      /* Term in coordinate conversion formula - GP to Y */
 
343
  double t3;      /* Term in coordinate conversion formula - GP to Y */
 
344
  double t4;      /* Term in coordinate conversion formula - GP to Y */
 
345
  double t5;      /* Term in coordinate conversion formula - GP to Y */
 
346
  double t6;      /* Term in coordinate conversion formula - GP to Y */
 
347
  double t7;      /* Term in coordinate conversion formula - GP to Y */
 
348
  double t8;      /* Term in coordinate conversion formula - GP to Y */
 
349
  double t9;      /* Term in coordinate conversion formula - GP to Y */
 
350
  double tmd;     /* True Meridional distance                        */
 
351
  double tmdo;    /* True Meridional distance for latitude of origin */
 
352
  long    Error_Code = TRANMERC_NO_ERROR;
 
353
  double temp_Origin;
 
354
  double temp_Long;
 
355
 
 
356
  if ((Latitude < -MAX_LAT) || (Latitude > MAX_LAT))
 
357
  {  /* Latitude out of range */
 
358
    Error_Code|= TRANMERC_LAT_ERROR;
 
359
  }
 
360
  if (Longitude > PI)
 
361
    Longitude -= (2 * PI);
 
362
  if ((Longitude < (TranMerc_Origin_Long - MAX_DELTA_LONG))
 
363
      || (Longitude > (TranMerc_Origin_Long + MAX_DELTA_LONG)))
 
364
  {
 
365
    if (Longitude < 0)
 
366
      temp_Long = Longitude + 2 * PI;
 
367
    else
 
368
      temp_Long = Longitude;
 
369
    if (TranMerc_Origin_Long < 0)
 
370
      temp_Origin = TranMerc_Origin_Long + 2 * PI;
 
371
    else
 
372
      temp_Origin = TranMerc_Origin_Long;
 
373
    if ((temp_Long < (temp_Origin - MAX_DELTA_LONG))
 
374
        || (temp_Long > (temp_Origin + MAX_DELTA_LONG)))
 
375
      Error_Code|= TRANMERC_LON_ERROR;
 
376
  }
 
377
  if (!Error_Code)
 
378
  { /* no errors */
 
379
 
 
380
    /* 
 
381
     *  Delta Longitude
 
382
     */
 
383
    dlam = Longitude - TranMerc_Origin_Long;
 
384
 
 
385
    if (fabs(dlam) > (9.0 * PI / 180))
 
386
    { /* Distortion will result if Longitude is more than 9 degrees from the Central Meridian */
 
387
      Error_Code |= TRANMERC_LON_WARNING;
 
388
    }
 
389
 
 
390
    if (dlam > PI)
 
391
      dlam -= (2 * PI);
 
392
    if (dlam < -PI)
 
393
      dlam += (2 * PI);
 
394
    if (fabs(dlam) < 2.e-10)
 
395
      dlam = 0.0;
 
396
 
 
397
    s = sin(Latitude);
 
398
    c = cos(Latitude);
 
399
    c2 = c * c;
 
400
    c3 = c2 * c;
 
401
    c5 = c3 * c2;
 
402
    c7 = c5 * c2;
 
403
    t = tan (Latitude);
 
404
    tan2 = t * t;
 
405
    tan3 = tan2 * t;
 
406
    tan4 = tan3 * t;
 
407
    tan5 = tan4 * t;
 
408
    tan6 = tan5 * t;
 
409
    eta = TranMerc_ebs * c2;
 
410
    eta2 = eta * eta;
 
411
    eta3 = eta2 * eta;
 
412
    eta4 = eta3 * eta;
 
413
 
 
414
    /* radius of curvature in prime vertical */
 
415
    sn = SPHSN(Latitude);
 
416
 
 
417
    /* True Meridianal Distances */
 
418
    tmd = SPHTMD(Latitude);
 
419
 
 
420
    /*  Origin  */
 
421
    tmdo = SPHTMD (TranMerc_Origin_Lat);
 
422
 
 
423
    /* northing */
 
424
    t1 = (tmd - tmdo) * TranMerc_Scale_Factor;
 
425
    t2 = sn * s * c * TranMerc_Scale_Factor/ 2.e0;
 
426
    t3 = sn * s * c3 * TranMerc_Scale_Factor * (5.e0 - tan2 + 9.e0 * eta 
 
427
                                                + 4.e0 * eta2) /24.e0; 
 
428
 
 
429
    t4 = sn * s * c5 * TranMerc_Scale_Factor * (61.e0 - 58.e0 * tan2
 
430
                                                + tan4 + 270.e0 * eta - 330.e0 * tan2 * eta + 445.e0 * eta2
 
431
                                                + 324.e0 * eta3 -680.e0 * tan2 * eta2 + 88.e0 * eta4 
 
432
                                                -600.e0 * tan2 * eta3 - 192.e0 * tan2 * eta4) / 720.e0;
 
433
 
 
434
    t5 = sn * s * c7 * TranMerc_Scale_Factor * (1385.e0 - 3111.e0 * 
 
435
                                                tan2 + 543.e0 * tan4 - tan6) / 40320.e0;
 
436
 
 
437
    *Northing = TranMerc_False_Northing + t1 + pow(dlam,2.e0) * t2
 
438
                + pow(dlam,4.e0) * t3 + pow(dlam,6.e0) * t4
 
439
                + pow(dlam,8.e0) * t5; 
 
440
 
 
441
    /* Easting */
 
442
    t6 = sn * c * TranMerc_Scale_Factor;
 
443
    t7 = sn * c3 * TranMerc_Scale_Factor * (1.e0 - tan2 + eta ) /6.e0;
 
444
    t8 = sn * c5 * TranMerc_Scale_Factor * (5.e0 - 18.e0 * tan2 + tan4
 
445
                                            + 14.e0 * eta - 58.e0 * tan2 * eta + 13.e0 * eta2 + 4.e0 * eta3 
 
446
                                            - 64.e0 * tan2 * eta2 - 24.e0 * tan2 * eta3 )/ 120.e0;
 
447
    t9 = sn * c7 * TranMerc_Scale_Factor * ( 61.e0 - 479.e0 * tan2
 
448
                                             + 179.e0 * tan4 - tan6 ) /5040.e0;
 
449
 
 
450
    *Easting = TranMerc_False_Easting + dlam * t6 + pow(dlam,3.e0) * t7 
 
451
               + pow(dlam,5.e0) * t8 + pow(dlam,7.e0) * t9;
 
452
  }
 
453
  return (Error_Code);
 
454
} /* END OF Convert_Geodetic_To_Transverse_Mercator */
 
455
 
 
456
 
 
457
long Convert_Transverse_Mercator_To_Geodetic (
 
458
                                             double Easting,
 
459
                                             double Northing,
 
460
                                             double *Latitude,
 
461
                                             double *Longitude)
 
462
{      /* BEGIN Convert_Transverse_Mercator_To_Geodetic */
 
463
 
 
464
  /*
 
465
   * The function Convert_Transverse_Mercator_To_Geodetic converts Transverse
 
466
   * Mercator projection (easting and northing) coordinates to geodetic
 
467
   * (latitude and longitude) coordinates, according to the current ellipsoid
 
468
   * and Transverse Mercator projection parameters.  If any errors occur, the
 
469
   * error code(s) are returned by the function, otherwise TRANMERC_NO_ERROR is
 
470
   * returned.
 
471
   *
 
472
   *    Easting       : Easting/X in meters                         (input)
 
473
   *    Northing      : Northing/Y in meters                        (input)
 
474
   *    Latitude      : Latitude in radians                         (output)
 
475
   *    Longitude     : Longitude in radians                        (output)
 
476
   */
 
477
 
 
478
  double c;       /* Cosine of latitude                          */
 
479
  double de;      /* Delta easting - Difference in Easting (Easting-Fe)    */
 
480
  double dlam;    /* Delta longitude - Difference in Longitude       */
 
481
  double eta;     /* constant - TranMerc_ebs *c *c                   */
 
482
  double eta2;
 
483
  double eta3;
 
484
  double eta4;
 
485
  double ftphi;   /* Footpoint latitude                              */
 
486
  int    i;       /* Loop iterator                   */
 
487
  double s;       /* Sine of latitude                        */
 
488
  double sn;      /* Radius of curvature in the prime vertical       */
 
489
  double sr;      /* Radius of curvature in the meridian             */
 
490
  double t;       /* Tangent of latitude                             */
 
491
  double tan2;
 
492
  double tan4;
 
493
  double t10;     /* Term in coordinate conversion formula - GP to Y */
 
494
  double t11;     /* Term in coordinate conversion formula - GP to Y */
 
495
  double t12;     /* Term in coordinate conversion formula - GP to Y */
 
496
  double t13;     /* Term in coordinate conversion formula - GP to Y */
 
497
  double t14;     /* Term in coordinate conversion formula - GP to Y */
 
498
  double t15;     /* Term in coordinate conversion formula - GP to Y */
 
499
  double t16;     /* Term in coordinate conversion formula - GP to Y */
 
500
  double t17;     /* Term in coordinate conversion formula - GP to Y */
 
501
  double tmd;     /* True Meridional distance                        */
 
502
  double tmdo;    /* True Meridional distance for latitude of origin */
 
503
  long Error_Code = TRANMERC_NO_ERROR;
 
504
 
 
505
  if ((Easting < (TranMerc_False_Easting - TranMerc_Delta_Easting))
 
506
      ||(Easting > (TranMerc_False_Easting + TranMerc_Delta_Easting)))
 
507
  { /* Easting out of range  */
 
508
    Error_Code |= TRANMERC_EASTING_ERROR;
 
509
  }
 
510
  if ((Northing < (TranMerc_False_Northing - TranMerc_Delta_Northing))
 
511
      || (Northing > (TranMerc_False_Northing + TranMerc_Delta_Northing)))
 
512
  { /* Northing out of range */
 
513
    Error_Code |= TRANMERC_NORTHING_ERROR;
 
514
  }
 
515
 
 
516
  if (!Error_Code)
 
517
  {
 
518
    /* True Meridional Distances for latitude of origin */
 
519
    tmdo = SPHTMD(TranMerc_Origin_Lat);
 
520
 
 
521
    /*  Origin  */
 
522
    tmd = tmdo +  (Northing - TranMerc_False_Northing) / TranMerc_Scale_Factor; 
 
523
 
 
524
    /* First Estimate */
 
525
    sr = SPHSR(0.e0);
 
526
    ftphi = tmd/sr;
 
527
 
 
528
    for (i = 0; i < 5 ; i++)
 
529
    {
 
530
      t10 = SPHTMD (ftphi);
 
531
      sr = SPHSR(ftphi);
 
532
      ftphi = ftphi + (tmd - t10) / sr;
 
533
    }
 
534
 
 
535
    /* Radius of Curvature in the meridian */
 
536
    sr = SPHSR(ftphi);
 
537
 
 
538
    /* Radius of Curvature in the meridian */
 
539
    sn = SPHSN(ftphi);
 
540
 
 
541
    /* Sine Cosine terms */
 
542
    s = sin(ftphi);
 
543
    c = cos(ftphi);
 
544
 
 
545
    /* Tangent Value  */
 
546
    t = tan(ftphi);
 
547
    tan2 = t * t;
 
548
    tan4 = tan2 * tan2;
 
549
    eta = TranMerc_ebs * pow(c,2);
 
550
    eta2 = eta * eta;
 
551
    eta3 = eta2 * eta;
 
552
    eta4 = eta3 * eta;
 
553
    de = Easting - TranMerc_False_Easting;
 
554
    if (fabs(de) < 0.0001)
 
555
      de = 0.0;
 
556
 
 
557
    /* Latitude */
 
558
    t10 = t / (2.e0 * sr * sn * pow(TranMerc_Scale_Factor, 2));
 
559
    t11 = t * (5.e0  + 3.e0 * tan2 + eta - 4.e0 * pow(eta,2)
 
560
               - 9.e0 * tan2 * eta) / (24.e0 * sr * pow(sn,3) 
 
561
                                       * pow(TranMerc_Scale_Factor,4));
 
562
    t12 = t * (61.e0 + 90.e0 * tan2 + 46.e0 * eta + 45.E0 * tan4
 
563
               - 252.e0 * tan2 * eta  - 3.e0 * eta2 + 100.e0 
 
564
               * eta3 - 66.e0 * tan2 * eta2 - 90.e0 * tan4
 
565
               * eta + 88.e0 * eta4 + 225.e0 * tan4 * eta2
 
566
               + 84.e0 * tan2* eta3 - 192.e0 * tan2 * eta4)
 
567
          / ( 720.e0 * sr * pow(sn,5) * pow(TranMerc_Scale_Factor, 6) );
 
568
    t13 = t * ( 1385.e0 + 3633.e0 * tan2 + 4095.e0 * tan4 + 1575.e0 
 
569
                * pow(t,6))/ (40320.e0 * sr * pow(sn,7) * pow(TranMerc_Scale_Factor,8));
 
570
    *Latitude = ftphi - pow(de,2) * t10 + pow(de,4) * t11 - pow(de,6) * t12 
 
571
                + pow(de,8) * t13;
 
572
 
 
573
    t14 = 1.e0 / (sn * c * TranMerc_Scale_Factor);
 
574
 
 
575
    t15 = (1.e0 + 2.e0 * tan2 + eta) / (6.e0 * pow(sn,3) * c * 
 
576
                                        pow(TranMerc_Scale_Factor,3));
 
577
 
 
578
    t16 = (5.e0 + 6.e0 * eta + 28.e0 * tan2 - 3.e0 * eta2
 
579
           + 8.e0 * tan2 * eta + 24.e0 * tan4 - 4.e0 
 
580
           * eta3 + 4.e0 * tan2 * eta2 + 24.e0 
 
581
           * tan2 * eta3) / (120.e0 * pow(sn,5) * c  
 
582
                             * pow(TranMerc_Scale_Factor,5));
 
583
 
 
584
    t17 = (61.e0 +  662.e0 * tan2 + 1320.e0 * tan4 + 720.e0 
 
585
           * pow(t,6)) / (5040.e0 * pow(sn,7) * c 
 
586
                          * pow(TranMerc_Scale_Factor,7));
 
587
 
 
588
    /* Difference in Longitude */
 
589
    dlam = de * t14 - pow(de,3) * t15 + pow(de,5) * t16 - pow(de,7) * t17;
 
590
 
 
591
    /* Longitude */
 
592
    (*Longitude) = TranMerc_Origin_Long + dlam;
 
593
 
 
594
    if((fabs)(*Latitude) > (90.0 * PI / 180.0))
 
595
      Error_Code |= TRANMERC_NORTHING_ERROR;
 
596
 
 
597
    if((*Longitude) > (PI))
 
598
    {
 
599
      *Longitude -= (2 * PI);
 
600
      if((fabs)(*Longitude) > PI)
 
601
        Error_Code |= TRANMERC_EASTING_ERROR;
 
602
    }
 
603
    else if((*Longitude) < (-PI))
 
604
    {
 
605
      *Longitude += (2 * PI);
 
606
      if((fabs)(*Longitude) > PI)
 
607
        Error_Code |= TRANMERC_EASTING_ERROR;
 
608
    }
 
609
 
 
610
    if (fabs(dlam) > (9.0 * PI / 180) * cos(*Latitude))
 
611
    { /* Distortion will result if Longitude is more than 9 degrees from the Central Meridian at the equator */
 
612
      /* and decreases to 0 degrees at the poles */
 
613
      /* As you move towards the poles, distortion will become more significant */
 
614
      Error_Code |= TRANMERC_LON_WARNING;
 
615
    }
 
616
  }
 
617
  return (Error_Code);
 
618
} /* END OF Convert_Transverse_Mercator_To_Geodetic */