1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
|
/****************************************************************************/
/* */
/* ./gps/tm.c - Convert to selected datum */
/* */
/* This file is part of gpstrans - a program to communicate with garmin gps */
/* Parts are taken from John F. Waers (jfwaers@csn.net) program MacGPS. */
/* */
/* */
/* Copyright (c) 1995 by Carsten Tschach (tschach@zedat.fu-berlin.de) */
/* */
/* Permission to use, copy, modify, and distribute this software and its */
/* documentation for non-commercial purpose, is hereby granted without fee, */
/* providing that the copyright notice appear in all copies and that both */
/* the copyright notice and this permission notice appear in supporting */
/* documentation. I make no representations about the suitability of this */
/* software for any purpose. It is provides "as is" without express or */
/* implid warranty. */
/* */
/****************************************************************************/
#include "defs.h"
#include "Garmin.h"
#include "Prefs.h"
#include <math.h>
/****************************************************************************/
/* Equations from "Map Projections -- A Working Manual", USGS Professional */
/* Paper 1395 */
/****************************************************************************/
/* prototype function */
static double M(double phi, double a, double es);
/****************************************************************************/
/* Convert latitude and longitude to converted position. */
/****************************************************************************/
void toTM(double lat, double lon, double lat0, double lon0, double k0,
double *x, double *y)
{
extern struct PREFS gPrefs;
double m, et2, n, t, c, A, a, m0, es, lambda, phi, lambda0, phi0;
datumParams(gPrefs.datum, &a, &es);
lambda = lon * Degree;
phi = lat * Degree;
phi0 = lat0 * Degree;
lambda0 = lon0 * Degree;
m0 = M(phi0, a, es);
m = M(phi, a, es);
et2 = es / (1 - es);
n = a / sqrt(1 - es * pow(sin(phi), 2.0));
t = pow(tan(phi), 2.0);
c = et2 * pow(cos(phi), 2.0);
A = (lambda - lambda0) * cos(phi);
*x = k0*n*(A + (1.0 - t + c)*A*A*A/6.0
+ (5.0 - 18.0*t + t*t + 72.0*c - 58.0*et2)*pow(A, 5.0) / 120.0);
*y = k0*(m - m0 + n*tan(phi)*(A*A/2.0
+ (5.0 - t + 9.0*c + 4*c*c)*pow(A, 4.0)/24.0
+ (61.0 - 58.0*t + t*t + 600.0*c - 330.0*et2)*
pow(A, 6.0)/720.0) );
}
/****************************************************************************/
/* Convert converted position to latitude and longitude. */
/****************************************************************************/
void fromTM(double x, double y, double lat0, double lon0, double k0,
double *lat, double *lon)
{
extern struct PREFS gFilePrefs;
double a, m0, es, et2, m, e1, mu, phi1, c1, t1, n1, r1, d, phi0, lambda0;
phi0 = lat0 * Degree;
lambda0 = lon0 * Degree;
datumParams(gFilePrefs.datum, &a, &es);
m0 = M(phi0, a, es);
et2 = es / (1.0 - es);
m = m0 + y / k0;
e1 = (1.0 - sqrt(1.0 - es)) / (1.0 + sqrt(1.0 - es));
mu = m / (a * (1.0 - es/4.0 - 3.0 * es*es/64.0 - 5.0 * es*es*es/256.0));
phi1 = mu + (3.0 * e1/2.0 - 27.0 * pow(e1, 3.0)/32.0) * sin(2.0 * mu)
+ (21.0 * e1*e1/16.0 - 55.0 * pow(e1, 4.0)/32.0)
* sin(4.0 * mu) + 151.0 * pow(e1, 3.0)/96.0 * sin(6.0 * mu)
+ 1097.0 * pow(e1, 4.0)/512.0 * sin(8.0 * mu);
c1 = et2 * pow(cos(phi1), 2.0);
t1 = pow(tan(phi1), 2.0);
n1 = a / sqrt(1 - es * pow(sin(phi1), 2.0));
r1 = a * (1.0 - es) / pow(1.0 - es * pow(sin(phi1), 2.0), 1.5);
d = x / (n1 * k0);
*lat = (phi1 - n1 * tan(phi1) / r1
* (d*d / 2.0 - (5.0 + 3.0 * t1 + 10.0 * c1 - 4.0 * c1*c1 - 9.0 * et2)
* pow(d, 4.0) / 24.0 + (61.0 + 90.0 * t1 + 298.0 * c1 + 45.0 * t1*t1
- 252.0 * et2 - 3.0 * c1*c1) * pow(d, 6.0) / 720.0 )) / Degree;
*lon = (lambda0 + (d - (1.0 + 2.0 * t1 + c1) * pow(d, 3.0)/6.0
+ (5.0 -2.0 * c1 + 28.0 * t1 - 3.0 * c1*c1 + 8.0 * et2 + 24.0 * t1*t1)
* pow(d, 5.0)/120.0) / cos(phi1)) / Degree;
}
/****************************************************************************/
/* Do some calculations - don't really know whats happen here. */
/****************************************************************************/
static double M(double phi, double a, double es)
{
double fix;
if (phi == 0.0)
return 0.0;
else {
fix = a * (
(1.0 - es/4.0 - 3.0*es*es/64.0 - 5.0*es*es*es/256.0 ) * phi -
(3.0 * es/8.0 + 3.0*es*es/32.0 + 45.0*es*es*es/1024.0 ) *
sin(2.0 * phi) +
(15.0 * es*es/256.0 + 45.0*es*es*es/1024.0 ) * sin(4.0 * phi) -
(35.0*es*es*es/3072.0 ) * sin(6.0 * phi) );
return fix;
}
}
|
