svn-gvsig-desktop / tags / v1_9_Build_1246 / libraries / libjni-proj4 / src / PJ_vandg.c @ 33782
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)
|