svn-gvsig-desktop / tags / v1_9_Build_1226 / libraries / libjni-proj4 / src / PJ_nsper.c @ 47796
History | View | Annotate | Download (3.25 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_nsper.c 4.1 94/02/15 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double height; \
|
6 |
double sinph0; \
|
7 |
double cosph0; \
|
8 |
double p; \
|
9 |
double rp; \
|
10 |
double pn1; \
|
11 |
double pfact; \
|
12 |
double h; \
|
13 |
double cg; \
|
14 |
double sg; \
|
15 |
double sw; \
|
16 |
double cw; \
|
17 |
int mode; \
|
18 |
int tilt;
|
19 |
#define PJ_LIB__
|
20 |
#include <projects.h> |
21 |
PROJ_HEAD(nsper, "Near-sided perspective") "\n\tAzi, Sph\n\th="; |
22 |
PROJ_HEAD(tpers, "Tilted perspective") "\n\tAzi, Sph\n\ttilt= azi= h="; |
23 |
# define EPS10 1.e-10 |
24 |
# define N_POLE 0 |
25 |
# define S_POLE 1 |
26 |
# define EQUIT 2 |
27 |
# define OBLIQ 3 |
28 |
FORWARD(s_forward); /* spheroid */
|
29 |
double coslam, cosphi, sinphi;
|
30 |
|
31 |
sinphi = sin(lp.phi); |
32 |
cosphi = cos(lp.phi); |
33 |
coslam = cos(lp.lam); |
34 |
switch (P->mode) {
|
35 |
case OBLIQ:
|
36 |
xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam; |
37 |
break;
|
38 |
case EQUIT:
|
39 |
xy.y = cosphi * coslam; |
40 |
break;
|
41 |
case S_POLE:
|
42 |
xy.y = - sinphi; |
43 |
break;
|
44 |
case N_POLE:
|
45 |
xy.y = sinphi; |
46 |
break;
|
47 |
} |
48 |
if (xy.y < P->rp) F_ERROR;
|
49 |
xy.y = P->pn1 / (P->p - xy.y); |
50 |
xy.x = xy.y * cosphi * sin(lp.lam); |
51 |
switch (P->mode) {
|
52 |
case OBLIQ:
|
53 |
xy.y *= (P->cosph0 * sinphi - |
54 |
P->sinph0 * cosphi * coslam); |
55 |
break;
|
56 |
case EQUIT:
|
57 |
xy.y *= sinphi; |
58 |
break;
|
59 |
case N_POLE:
|
60 |
coslam = - coslam; |
61 |
case S_POLE:
|
62 |
xy.y *= cosphi * coslam; |
63 |
break;
|
64 |
} |
65 |
if (P->tilt) {
|
66 |
double yt, ba;
|
67 |
|
68 |
yt = xy.y * P->cg + xy.x * P->sg; |
69 |
ba = 1. / (yt * P->sw * P->h + P->cw);
|
70 |
xy.x = (xy.x * P->cg - xy.y * P->sg) * P->cw * ba; |
71 |
xy.y = yt * ba; |
72 |
} |
73 |
return (xy);
|
74 |
} |
75 |
INVERSE(s_inverse); /* spheroid */
|
76 |
double rh, cosz, sinz;
|
77 |
|
78 |
if (P->tilt) {
|
79 |
double bm, bq, yt;
|
80 |
|
81 |
yt = 1./(P->pn1 - xy.y * P->sw);
|
82 |
bm = P->pn1 * xy.x * yt; |
83 |
bq = P->pn1 * xy.y * P->cw * yt; |
84 |
xy.x = bm * P->cg + bq * P->sg; |
85 |
xy.y = bq * P->cg - bm * P->sg; |
86 |
} |
87 |
rh = hypot(xy.x, xy.y); |
88 |
if ((sinz = 1. - rh * rh * P->pfact) < 0.) I_ERROR; |
89 |
sinz = (P->p - sqrt(sinz)) / (P->pn1 / rh + rh / P->pn1); |
90 |
cosz = sqrt(1. - sinz * sinz);
|
91 |
if (fabs(rh) <= EPS10) {
|
92 |
lp.lam = 0.;
|
93 |
lp.phi = P->phi0; |
94 |
} else {
|
95 |
switch (P->mode) {
|
96 |
case OBLIQ:
|
97 |
lp.phi = asin(cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh); |
98 |
xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh; |
99 |
xy.x *= sinz * P->cosph0; |
100 |
break;
|
101 |
case EQUIT:
|
102 |
lp.phi = asin(xy.y * sinz / rh); |
103 |
xy.y = cosz * rh; |
104 |
xy.x *= sinz; |
105 |
break;
|
106 |
case N_POLE:
|
107 |
lp.phi = asin(cosz); |
108 |
xy.y = -xy.y; |
109 |
break;
|
110 |
case S_POLE:
|
111 |
lp.phi = - asin(cosz); |
112 |
break;
|
113 |
} |
114 |
lp.lam = atan2(xy.x, xy.y); |
115 |
} |
116 |
return (lp);
|
117 |
} |
118 |
FREEUP; if (P) pj_dalloc(P); }
|
119 |
static PJ *
|
120 |
setup(PJ *P) { |
121 |
if ((P->height = pj_param(P->params, "dh").f) <= 0.) E_ERROR(-30); |
122 |
if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
|
123 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
124 |
else if (fabs(P->phi0) < EPS10) |
125 |
P->mode = EQUIT; |
126 |
else {
|
127 |
P->mode = OBLIQ; |
128 |
P->sinph0 = sin(P->phi0); |
129 |
P->cosph0 = cos(P->phi0); |
130 |
} |
131 |
P->pn1 = P->height / P->a; /* normalize by radius */
|
132 |
P->p = 1. + P->pn1;
|
133 |
P->rp = 1. / P->p;
|
134 |
P->h = 1. / P->pn1;
|
135 |
P->pfact = (P->p + 1.) * P->h;
|
136 |
P->inv = s_inverse; |
137 |
P->fwd = s_forward; |
138 |
P->es = 0.;
|
139 |
return P;
|
140 |
} |
141 |
ENTRY0(nsper) |
142 |
P->tilt = 0;
|
143 |
ENDENTRY(setup(P)) |
144 |
ENTRY0(tpers) |
145 |
double omega, gamma;
|
146 |
|
147 |
omega = pj_param(P->params, "dtilt").f * DEG_TO_RAD;
|
148 |
gamma = pj_param(P->params, "dazi").f * DEG_TO_RAD;
|
149 |
P->tilt = 1;
|
150 |
P->cg = cos(gamma); P->sg = sin(gamma); |
151 |
P->cw = cos(omega); P->sw = sin(omega); |
152 |
ENDENTRY(setup(P)) |