svn-gvsig-desktop / tags / v2_0_0_Build_2006 / libraries / libjni-proj4 / src / PJ_vandg.c @ 42148
History | View | Annotate | Download (2.36 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_vandg.c 4.1 94/02/15 GIE REL"; |
3 |
#endif
|
4 |
#define PJ_LIB__
|
5 |
# include <projects.h> |
6 |
PROJ_HEAD(vandg, "van der Grinten (I)") "\n\tMisc Sph"; |
7 |
# define TOL 1.e-10 |
8 |
# define THIRD .33333333333333333333 |
9 |
# define TWO_THRD .66666666666666666666 |
10 |
# define C2_27 .07407407407407407407 |
11 |
# define PI4_3 4.18879020478639098458 |
12 |
# define PISQ 9.86960440108935861869 |
13 |
# define TPISQ 19.73920880217871723738 |
14 |
# define HPISQ 4.93480220054467930934 |
15 |
FORWARD(s_forward); /* spheroid */
|
16 |
double al, al2, g, g2, p2;
|
17 |
|
18 |
p2 = fabs(lp.phi / HALFPI); |
19 |
if ((p2 - TOL) > 1.) F_ERROR; |
20 |
if (p2 > 1.) |
21 |
p2 = 1.;
|
22 |
if (fabs(lp.phi) <= TOL) {
|
23 |
xy.x = lp.lam; |
24 |
xy.y = 0.;
|
25 |
} else if (fabs(lp.lam) <= TOL || fabs(p2 - 1.) < TOL) { |
26 |
xy.x = 0.;
|
27 |
xy.y = PI * tan(.5 * asin(p2));
|
28 |
if (lp.phi < 0.) xy.y = -xy.y; |
29 |
} else {
|
30 |
al = .5 * fabs(PI / lp.lam - lp.lam / PI);
|
31 |
al2 = al * al; |
32 |
g = sqrt(1. - p2 * p2);
|
33 |
g = g / (p2 + g - 1.);
|
34 |
g2 = g * g; |
35 |
p2 = g * (2. / p2 - 1.); |
36 |
p2 = p2 * p2; |
37 |
xy.x = g - p2; g = p2 + al2; |
38 |
xy.x = PI * (al * xy.x + sqrt(al2 * xy.x * xy.x - g * (g2 - p2))) / g; |
39 |
if (lp.lam < 0.) xy.x = -xy.x; |
40 |
xy.y = fabs(xy.x / PI); |
41 |
xy.y = 1. - xy.y * (xy.y + 2. * al); |
42 |
if (xy.y < -TOL) F_ERROR;
|
43 |
if (xy.y < 0.) xy.y = 0.; |
44 |
else xy.y = sqrt(xy.y) * (lp.phi < 0. ? -PI : PI); |
45 |
} |
46 |
return (xy);
|
47 |
} |
48 |
INVERSE(s_inverse); /* spheroid */
|
49 |
double t, c0, c1, c2, c3, al, r2, r, m, d, ay, x2, y2;
|
50 |
|
51 |
x2 = xy.x * xy.x; |
52 |
if ((ay = fabs(xy.y)) < TOL) {
|
53 |
lp.phi = 0.;
|
54 |
t = x2 * x2 + TPISQ * (x2 + HPISQ); |
55 |
lp.lam = fabs(xy.x) <= TOL ? 0. :
|
56 |
.5 * (x2 - PISQ + sqrt(t)) / xy.x;
|
57 |
return (lp);
|
58 |
} |
59 |
y2 = xy.y * xy.y; |
60 |
r = x2 + y2; r2 = r * r; |
61 |
c1 = - PI * ay * (r + PISQ); |
62 |
c3 = r2 + TWOPI * (ay * r + PI * (y2 + PI * (ay + HALFPI))); |
63 |
c2 = c1 + PISQ * (r - 3. * y2);
|
64 |
c0 = PI * ay; |
65 |
c2 /= c3; |
66 |
al = c1 / c3 - THIRD * c2 * c2; |
67 |
m = 2. * sqrt(-THIRD * al);
|
68 |
d = C2_27 * c2 * c2 * c2 + (c0 * c0 - THIRD * c2 * c1) / c3; |
69 |
if (((t = fabs(d = 3. * d / (al * m))) - TOL) <= 1.) { |
70 |
d = t > 1. ? (d > 0. ? 0. : PI) : acos(d); |
71 |
lp.phi = PI * (m * cos(d * THIRD + PI4_3) - THIRD * c2); |
72 |
if (xy.y < 0.) lp.phi = -lp.phi; |
73 |
t = r2 + TPISQ * (x2 - y2 + HPISQ); |
74 |
lp.lam = fabs(xy.x) <= TOL ? 0. :
|
75 |
.5 * (r - PISQ + (t <= 0. ? 0. : sqrt(t))) / xy.x; |
76 |
} else
|
77 |
I_ERROR; |
78 |
return (lp);
|
79 |
} |
80 |
FREEUP; if (P) pj_dalloc(P); }
|
81 |
ENTRY0(vandg) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
|