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
Set_Transverse_Mercator_Parameters(double a,double f,double Origin_Latitude,double Central_Meridian,double False_Easting,double False_Northing,double Scale_Factor)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
Get_Transverse_Mercator_Parameters(double * a,double * f,double * Origin_Latitude,double * Central_Meridian,double * False_Easting,double * False_Northing,double * Scale_Factor)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
Convert_Geodetic_To_Transverse_Mercator(double Latitude,double Longitude,double * Easting,double * Northing)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
Convert_Transverse_Mercator_To_Geodetic(double Easting,double Northing,double * Latitude,double * Longitude)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 */
619