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
|
/****************************************************************************/
/* */
/* ./grid/ups.c - Convert to and from UPS Grid Format */
/* */
/* 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>
/* prototype functions */
static void calcPhi(double *phi, double e, double t);
/****************************************************************************/
/* Convert degree to UPS Grid. */
/****************************************************************************/
void toUPS(double lat, double lon, double *x, double *y)
{
extern struct PREFS gPrefs;
double a, t, e, es, rho;
const double k0 = 0.994;
double lambda = lon * Degree;
double phi = fabs(lat * Degree);
datumParams(gPrefs.datum, &a, &es);
e = sqrt(es);
t = tan(Pi/4.0 - phi/2.0) / pow( (1.0 - e * sin(phi)) /
(1.0 + e * sin(phi)), (e/2.0) );
rho = 2.0 * a * k0 * t / sqrt(pow(1.0+e, 1.0+e) * pow(1.0-e, 1.0-e));
*x = rho * sin(lambda);
*y = rho * cos(lambda);
if (lat > 0.0) *y = -(*y); /* Northern hemisphere */
*x += 2.0e6; /* Add in false easting and northing */
*y += 2.0e6;
}
/****************************************************************************/
/* Convert UPS Grid to degree. */
/****************************************************************************/
void fromUPS(short southernHemisphere, double x, double y, double *lat,
double *lon)
{
extern struct PREFS gFilePrefs;
double a, es, e, t, rho;
const double k0 = 0.994;
datumParams(gFilePrefs.datum, &a, &es);
e = sqrt(es);
/* Remove false easting and northing */
x -= 2.0e6;
y -= 2.0e6;
rho = sqrt(x*x + y*y);
t = rho * sqrt(pow(1.0+e, 1.0+e) * pow(1.0-e, 1.0-e)) / (2.0 * a * k0);
calcPhi(lat, e, t);
*lat /= Degree;
if (y != 0.0)
t = atan(fabs(x/y));
else {
t = Pi / 2.0;
if (x < 0.0) t = -t;
}
if (!southernHemisphere)
y = -y;
if (y < 0.0) t = Pi - t;
if (x < 0.0) t = -t;
*lon = t / Degree;
}
/****************************************************************************/
/* Calculate Phi. */
/****************************************************************************/
static void calcPhi(double *phi, double e, double t)
{
double old = Pi/2.0 - 2.0 * atan(t);
short maxIterations = 20;
while ((fabs((*phi - old) / *phi) > 1.0e-8) && maxIterations-- ) {
old = *phi;
*phi = Pi/ 2.0 - 2.0 * atan( t * pow((1.0 - e * sin(*phi)) /
((1.0 + e * sin(*phi))), (e / 2.0)) );
}
}
|
