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
|
/**
* Navit, a modular navigation system.
* Copyright (C) 2005-2008 Navit Team
*
* This program is free software; you can redistribute it and/or
* modify it under the terms of the GNU General Public License
* version 2 as published by the Free Software Foundation.
*
* This program is distributed in the hope that it will be useful,
* but WITHOUT ANY WARRANTY; without even the implied warranty of
* MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE. See the
* GNU General Public License for more details.
*
* You should have received a copy of the GNU General Public License
* along with this program; if not, write to the
* Free Software Foundation, Inc., 51 Franklin Street, Fifth Floor,
* Boston, MA 02110-1301, USA.
*/
#ifndef lint
static const char SCCSID[] =
"@(#)geod_inv.c 4.5 95/09/23 GIE REL";
#endif
# include "projects.h"
# include "geodesic.h"
# define DTOL 1e-12
void geod_inv(void)
{
double th1,
th2,
thm,
dthm,
dlamm,
dlam,
sindlamm,
costhm,
sinthm,
cosdthm,
sindthm, L, E, cosd, d, X, Y, T, sind, tandlammp, u, v, D, A, B;
if (ellipse) {
th1 = atan(onef * tan(phi1));
th2 = atan(onef * tan(phi2));
} else {
th1 = phi1;
th2 = phi2;
}
thm = .5 * (th1 + th2);
dthm = .5 * (th2 - th1);
dlamm = .5 * (dlam = adjlon(lam2 - lam1));
if (fabs(dlam) < DTOL && fabs(dthm) < DTOL) {
al12 = al21 = geod_S = 0.;
return;
}
sindlamm = sin(dlamm);
costhm = cos(thm);
sinthm = sin(thm);
cosdthm = cos(dthm);
sindthm = sin(dthm);
L = sindthm * sindthm + (cosdthm * cosdthm - sinthm * sinthm)
* sindlamm * sindlamm;
d = acos(cosd = 1 - L - L);
if (ellipse) {
E = cosd + cosd;
sind = sin(d);
Y = sinthm * cosdthm;
Y *= (Y + Y) / (1. - L);
T = sindthm * costhm;
T *= (T + T) / L;
X = Y + T;
Y -= T;
T = d / sind;
D = 4. * T * T;
A = D * E;
B = D + D;
geod_S = geod_a * sind * (T - f4 * (T * X - Y) +
f64 * (X * (A + (T - .5 * (A - E)) * X) -
Y * (B + E * Y) + D * X * Y));
tandlammp = tan(.5 * (dlam - .25 * (Y + Y - E * (4. - X)) *
(f2 * T + f64 * (32. * T - (20. * T - A)
* X - (B +
4.) * Y)) *
tan(dlam)));
} else {
geod_S = geod_a * d;
tandlammp = tan(dlamm);
}
u = atan2(sindthm, (tandlammp * costhm));
v = atan2(cosdthm, (tandlammp * sinthm));
al12 = adjlon(TWOPI + v - u);
al21 = adjlon(TWOPI - v - u);
}
|
