geod_for.c
2.13 KB
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
# include "projects.h"
# include "geodesic.h"
# define MERI_TOL 1e-9
static double
th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1;
static int
merid, signS;
void
geod_pre(void) {
al12 = adjlon(al12); /* reduce to +- 0-PI */
signS = fabs(al12) > HALFPI ? 1 : 0;
th1 = ellipse ? atan(onef * tan(phi1)) : phi1;
costh1 = cos(th1);
sinth1 = sin(th1);
if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {
sina12 = 0.;
cosa12 = fabs(al12) < HALFPI ? 1. : -1.;
M = 0.;
} else {
cosa12 = cos(al12);
M = costh1 * sina12;
}
N = costh1 * cosa12;
if (ellipse) {
if (merid) {
c1 = 0.;
c2 = f4;
D = 1. - c2;
D *= D;
P = c2 / D;
} else {
c1 = geod_f * M;
c2 = f4 * (1. - M * M);
D = (1. - c2)*(1. - c2 - c1 * M);
P = (1. + .5 * c1 * M) * c2 / D;
}
}
if (merid) s1 = HALFPI - th1;
else {
s1 = (fabs(M) >= 1.) ? 0. : acos(M);
s1 = sinth1 / sin(s1);
s1 = (fabs(s1) >= 1.) ? 0. : acos(s1);
}
}
void
geod_for(void) {
double d,sind,u,V,X,ds,cosds,sinds,ss,de;
if (ellipse) {
d = geod_S / (D * geod_a);
if (signS) d = -d;
u = 2. * (s1 - d);
V = cos(u + d);
X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.);
ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind;
ss = s1 + s1 - ds;
} else {
ds = geod_S / geod_a;
if (signS) ds = - ds;
}
cosds = cos(ds);
sinds = sin(ds);
if (signS) sinds = - sinds;
al21 = N * cosds - sinth1 * sinds;
if (merid) {
phi2 = atan( tan(HALFPI + s1 - ds) / onef);
if (al21 > 0.) {
al21 = PI;
if (signS)
de = PI;
else {
phi2 = - phi2;
de = 0.;
}
} else {
al21 = 0.;
if (signS) {
phi2 = - phi2;
de = 0;
} else
de = PI;
}
} else {
al21 = atan(M / al21);
if (al21 > 0)
al21 += PI;
if (al12 < 0.)
al21 -= PI;
al21 = adjlon(al21);
phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) /
(ellipse ? onef * M : M));
de = atan2(sinds * sina12 ,
(costh1 * cosds - sinth1 * sinds * cosa12));
if (ellipse) {
if (signS)
de += c1 * ((1. - c2) * ds +
c2 * sinds * cos(ss));
else
de -= c1 * ((1. - c2) * ds -
c2 * sinds * cos(ss));
}
}
lam2 = adjlon( lam1 + de );
}