svn-gvsig-desktop / trunk / libraries / libjni-proj4 / src / PJ_ortho.c @ 7098
History | View | Annotate | Download (2.45 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_ortho.c 4.1 94/02/15 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double sinph0; \
|
6 |
double cosph0; \
|
7 |
int mode;
|
8 |
#define PJ_LIB__
|
9 |
#include <projects.h> |
10 |
PROJ_HEAD(ortho, "Orthographic") "\n\tAzi, Sph."; |
11 |
#define EPS10 1.e-10 |
12 |
#define N_POLE 0 |
13 |
#define S_POLE 1 |
14 |
#define EQUIT 2 |
15 |
#define OBLIQ 3 |
16 |
FORWARD(s_forward); /* spheroid */
|
17 |
double coslam, cosphi, sinphi;
|
18 |
|
19 |
cosphi = cos(lp.phi); |
20 |
coslam = cos(lp.lam); |
21 |
switch (P->mode) {
|
22 |
case EQUIT:
|
23 |
if (cosphi * coslam < - EPS10) F_ERROR;
|
24 |
xy.y = sin(lp.phi); |
25 |
break;
|
26 |
case OBLIQ:
|
27 |
if (P->sinph0 * (sinphi = sin(lp.phi)) +
|
28 |
P->cosph0 * cosphi * coslam < - EPS10) F_ERROR; |
29 |
xy.y = P->cosph0 * sinphi - P->sinph0 * cosphi * coslam; |
30 |
break;
|
31 |
case N_POLE:
|
32 |
coslam = - coslam; |
33 |
case S_POLE:
|
34 |
if (fabs(lp.phi - P->phi0) - EPS10 > HALFPI) F_ERROR;
|
35 |
xy.y = cosphi * coslam; |
36 |
break;
|
37 |
} |
38 |
xy.x = cosphi * sin(lp.lam); |
39 |
return (xy);
|
40 |
} |
41 |
|
42 |
INVERSE(s_inverse); /* spheroid */
|
43 |
double rh, cosc, sinc;
|
44 |
|
45 |
if ((sinc = (rh = hypot(xy.x, xy.y))) > 1.) { |
46 |
if ((sinc - 1.) > EPS10) I_ERROR; |
47 |
sinc = 1.;
|
48 |
} |
49 |
cosc = sqrt(1. - sinc * sinc); /* in this range OK */ |
50 |
if (fabs(rh) <= EPS10) {
|
51 |
lp.phi = P->phi0; |
52 |
lp.lam = 0.0; |
53 |
} else {
|
54 |
switch (P->mode) {
|
55 |
case N_POLE:
|
56 |
xy.y = -xy.y; |
57 |
lp.phi = acos(sinc); |
58 |
break;
|
59 |
case S_POLE:
|
60 |
lp.phi = - acos(sinc); |
61 |
break;
|
62 |
case EQUIT:
|
63 |
lp.phi = xy.y * sinc / rh; |
64 |
xy.x *= sinc; |
65 |
xy.y = cosc * rh; |
66 |
goto sinchk;
|
67 |
case OBLIQ:
|
68 |
lp.phi = cosc * P->sinph0 + xy.y * sinc * P->cosph0 /rh; |
69 |
xy.y = (cosc - P->sinph0 * lp.phi) * rh; |
70 |
xy.x *= sinc * P->cosph0; |
71 |
sinchk:
|
72 |
if (fabs(lp.phi) >= 1.) |
73 |
lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
|
74 |
else
|
75 |
lp.phi = asin(lp.phi); |
76 |
break;
|
77 |
} |
78 |
lp.lam = (xy.y == 0. && (P->mode == OBLIQ || P->mode == EQUIT))
|
79 |
? (xy.x == 0. ? 0. : xy.x < 0. ? -HALFPI : HALFPI) |
80 |
: atan2(xy.x, xy.y); |
81 |
} |
82 |
return (lp);
|
83 |
} |
84 |
|
85 |
FREEUP; if (P) pj_dalloc(P); }
|
86 |
ENTRY0(ortho) |
87 |
if (fabs(fabs(P->phi0) - HALFPI) <= EPS10)
|
88 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
89 |
else if (fabs(P->phi0) > EPS10) { |
90 |
P->mode = OBLIQ; |
91 |
P->sinph0 = sin(P->phi0); |
92 |
P->cosph0 = cos(P->phi0); |
93 |
} else
|
94 |
P->mode = EQUIT; |
95 |
P->inv = s_inverse; |
96 |
P->fwd = s_forward; |
97 |
P->es = 0.;
|
98 |
ENDENTRY(P) |