svn-gvsig-desktop / tags / v1_1_2_Build_1045 / libraries / libjni-proj4 / src / geod_for.c @ 38629
History | View | Annotate | Download (2.2 KB)
| 1 |
#ifndef lint
|
|---|---|
| 2 |
static const char SCCSID[]="@(#)geod_for.c 4.6 95/09/23 GIE REL"; |
| 3 |
#endif
|
| 4 |
# include "projects.h" |
| 5 |
# include "geodesic.h" |
| 6 |
# define MERI_TOL 1e-9 |
| 7 |
static double |
| 8 |
th1,costh1,sinth1,sina12,cosa12,M,N,c1,c2,D,P,s1; |
| 9 |
static int |
| 10 |
merid, signS; |
| 11 |
void
|
| 12 |
geod_pre(void) {
|
| 13 |
al12 = adjlon(al12); /* reduce to +- 0-PI */
|
| 14 |
signS = fabs(al12) > HALFPI ? 1 : 0; |
| 15 |
th1 = ellipse ? atan(onef * tan(phi1)) : phi1; |
| 16 |
costh1 = cos(th1); |
| 17 |
sinth1 = sin(th1); |
| 18 |
if ((merid = fabs(sina12 = sin(al12)) < MERI_TOL)) {
|
| 19 |
sina12 = 0.;
|
| 20 |
cosa12 = fabs(al12) < HALFPI ? 1. : -1.; |
| 21 |
M = 0.;
|
| 22 |
} else {
|
| 23 |
cosa12 = cos(al12); |
| 24 |
M = costh1 * sina12; |
| 25 |
} |
| 26 |
N = costh1 * cosa12; |
| 27 |
if (ellipse) {
|
| 28 |
if (merid) {
|
| 29 |
c1 = 0.;
|
| 30 |
c2 = f4; |
| 31 |
D = 1. - c2;
|
| 32 |
D *= D; |
| 33 |
P = c2 / D; |
| 34 |
} else {
|
| 35 |
c1 = geod_f * M; |
| 36 |
c2 = f4 * (1. - M * M);
|
| 37 |
D = (1. - c2)*(1. - c2 - c1 * M); |
| 38 |
P = (1. + .5 * c1 * M) * c2 / D; |
| 39 |
} |
| 40 |
} |
| 41 |
if (merid) s1 = HALFPI - th1;
|
| 42 |
else {
|
| 43 |
s1 = (fabs(M) >= 1.) ? 0. : acos(M); |
| 44 |
s1 = sinth1 / sin(s1); |
| 45 |
s1 = (fabs(s1) >= 1.) ? 0. : acos(s1); |
| 46 |
} |
| 47 |
} |
| 48 |
void
|
| 49 |
geod_for(void) {
|
| 50 |
double d,sind,u,V,X,ds,cosds,sinds,ss,de;
|
| 51 |
|
| 52 |
if (ellipse) {
|
| 53 |
d = geod_S / (D * geod_a); |
| 54 |
if (signS) d = -d;
|
| 55 |
u = 2. * (s1 - d);
|
| 56 |
V = cos(u + d); |
| 57 |
X = c2 * c2 * (sind = sin(d)) * cos(d) * (2. * V * V - 1.); |
| 58 |
ds = d + X - 2. * P * V * (1. - 2. * P * cos(u)) * sind; |
| 59 |
ss = s1 + s1 - ds; |
| 60 |
} else {
|
| 61 |
ds = geod_S / geod_a; |
| 62 |
if (signS) ds = - ds;
|
| 63 |
} |
| 64 |
cosds = cos(ds); |
| 65 |
sinds = sin(ds); |
| 66 |
if (signS) sinds = - sinds;
|
| 67 |
al21 = N * cosds - sinth1 * sinds; |
| 68 |
if (merid) {
|
| 69 |
phi2 = atan( tan(HALFPI + s1 - ds) / onef); |
| 70 |
if (al21 > 0.) { |
| 71 |
al21 = PI; |
| 72 |
if (signS)
|
| 73 |
de = PI; |
| 74 |
else {
|
| 75 |
phi2 = - phi2; |
| 76 |
de = 0.;
|
| 77 |
} |
| 78 |
} else {
|
| 79 |
al21 = 0.;
|
| 80 |
if (signS) {
|
| 81 |
phi2 = - phi2; |
| 82 |
de = 0;
|
| 83 |
} else
|
| 84 |
de = PI; |
| 85 |
} |
| 86 |
} else {
|
| 87 |
al21 = atan(M / al21); |
| 88 |
if (al21 > 0) |
| 89 |
al21 += PI; |
| 90 |
if (al12 < 0.) |
| 91 |
al21 -= PI; |
| 92 |
al21 = adjlon(al21); |
| 93 |
phi2 = atan(-(sinth1 * cosds + N * sinds) * sin(al21) / |
| 94 |
(ellipse ? onef * M : M)); |
| 95 |
de = atan2(sinds * sina12 , |
| 96 |
(costh1 * cosds - sinth1 * sinds * cosa12)); |
| 97 |
if (ellipse)
|
| 98 |
if (signS)
|
| 99 |
de += c1 * ((1. - c2) * ds +
|
| 100 |
c2 * sinds * cos(ss)); |
| 101 |
else
|
| 102 |
de -= c1 * ((1. - c2) * ds -
|
| 103 |
c2 * sinds * cos(ss)); |
| 104 |
} |
| 105 |
lam2 = adjlon( lam1 + de ); |
| 106 |
} |