svn-gvsig-desktop / tags / v1_9_Build_1222 / libraries / libjni-proj4 / src / PJ_gnom.c @ 41849
History | View | Annotate | Download (2.19 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_gnom.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(gnom, "Gnomonic") "\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 |
sinphi = sin(lp.phi); |
20 |
cosphi = cos(lp.phi); |
21 |
coslam = cos(lp.lam); |
22 |
switch (P->mode) {
|
23 |
case EQUIT:
|
24 |
xy.y = cosphi * coslam; |
25 |
break;
|
26 |
case OBLIQ:
|
27 |
xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam; |
28 |
break;
|
29 |
case S_POLE:
|
30 |
xy.y = - sinphi; |
31 |
break;
|
32 |
case N_POLE:
|
33 |
xy.y = sinphi; |
34 |
break;
|
35 |
} |
36 |
if (xy.y <= EPS10) F_ERROR;
|
37 |
xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam);
|
38 |
switch (P->mode) {
|
39 |
case EQUIT:
|
40 |
xy.y *= sinphi; |
41 |
break;
|
42 |
case OBLIQ:
|
43 |
xy.y *= P->cosph0 * sinphi - P->sinph0 * cosphi * coslam; |
44 |
break;
|
45 |
case N_POLE:
|
46 |
coslam = - coslam; |
47 |
case S_POLE:
|
48 |
xy.y *= cosphi * coslam; |
49 |
break;
|
50 |
} |
51 |
return (xy);
|
52 |
} |
53 |
INVERSE(s_inverse); /* spheroid */
|
54 |
double rh, cosz, sinz;
|
55 |
|
56 |
rh = hypot(xy.x, xy.y); |
57 |
sinz = sin(lp.phi = atan(rh)); |
58 |
cosz = sqrt(1. - sinz * sinz);
|
59 |
if (fabs(rh) <= EPS10) {
|
60 |
lp.phi = P->phi0; |
61 |
lp.lam = 0.;
|
62 |
} else {
|
63 |
switch (P->mode) {
|
64 |
case OBLIQ:
|
65 |
lp.phi = cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh; |
66 |
if (fabs(lp.phi) >= 1.) |
67 |
lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
|
68 |
else
|
69 |
lp.phi = asin(lp.phi); |
70 |
xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh; |
71 |
xy.x *= sinz * P->cosph0; |
72 |
break;
|
73 |
case EQUIT:
|
74 |
lp.phi = xy.y * sinz / rh; |
75 |
if (fabs(lp.phi) >= 1.) |
76 |
lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
|
77 |
else
|
78 |
lp.phi = asin(lp.phi); |
79 |
xy.y = cosz * rh; |
80 |
xy.x *= sinz; |
81 |
break;
|
82 |
case S_POLE:
|
83 |
lp.phi -= HALFPI; |
84 |
break;
|
85 |
case N_POLE:
|
86 |
lp.phi = HALFPI - lp.phi; |
87 |
xy.y = -xy.y; |
88 |
break;
|
89 |
} |
90 |
lp.lam = atan2(xy.x, xy.y); |
91 |
} |
92 |
return (lp);
|
93 |
} |
94 |
FREEUP; if (P) pj_dalloc(P); }
|
95 |
ENTRY0(gnom) |
96 |
if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
|
97 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
98 |
else if (fabs(P->phi0) < EPS10) |
99 |
P->mode = EQUIT; |
100 |
else {
|
101 |
P->mode = OBLIQ; |
102 |
P->sinph0 = sin(P->phi0); |
103 |
P->cosph0 = cos(P->phi0); |
104 |
} |
105 |
P->inv = s_inverse; |
106 |
P->fwd = s_forward; |
107 |
P->es = 0.;
|
108 |
ENDENTRY(P) |