svn-gvsig-desktop / tags / JCRS_V02_BN11 / libjni-proj4 / src / PJ_laea.c @ 32257
History | View | Annotate | Download (5 KB)
1 |
#ifndef lint
|
---|---|
2 |
static const char SCCSID[]="@(#)PJ_laea.c 4.1 94/02/15 GIE REL"; |
3 |
#endif
|
4 |
#define PROJ_PARMS__ \
|
5 |
double sinb1; \
|
6 |
double cosb1; \
|
7 |
double xmf; \
|
8 |
double ymf; \
|
9 |
double mmf; \
|
10 |
double qp; \
|
11 |
double dd; \
|
12 |
double rq; \
|
13 |
double *apa; \
|
14 |
int mode;
|
15 |
#define PJ_LIB__
|
16 |
#include <projects.h> |
17 |
PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell"; |
18 |
#define sinph0 P->sinb1
|
19 |
#define cosph0 P->cosb1
|
20 |
#define EPS10 1.e-10 |
21 |
#define NITER 20 |
22 |
#define CONV 1.e-10 |
23 |
#define N_POLE 0 |
24 |
#define S_POLE 1 |
25 |
#define EQUIT 2 |
26 |
#define OBLIQ 3 |
27 |
FORWARD(e_forward); /* ellipsoid */
|
28 |
double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0; |
29 |
|
30 |
coslam = cos(lp.lam); |
31 |
sinlam = sin(lp.lam); |
32 |
sinphi = sin(lp.phi); |
33 |
q = pj_qsfn(sinphi, P->e, P->one_es); |
34 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
35 |
sinb = q / P->qp; |
36 |
cosb = sqrt(1. - sinb * sinb);
|
37 |
} |
38 |
switch (P->mode) {
|
39 |
case OBLIQ:
|
40 |
b = 1. + P->sinb1 * sinb + P->cosb1 * cosb * coslam;
|
41 |
break;
|
42 |
case EQUIT:
|
43 |
b = 1. + cosb * coslam;
|
44 |
break;
|
45 |
case N_POLE:
|
46 |
b = HALFPI + lp.phi; |
47 |
q = P->qp - q; |
48 |
break;
|
49 |
case S_POLE:
|
50 |
b = lp.phi - HALFPI; |
51 |
q = P->qp + q; |
52 |
break;
|
53 |
} |
54 |
if (fabs(b) < EPS10) F_ERROR;
|
55 |
switch (P->mode) {
|
56 |
case OBLIQ:
|
57 |
xy.y = P->ymf * ( b = sqrt(2. / b) )
|
58 |
* (P->cosb1 * sinb - P->sinb1 * cosb * coslam); |
59 |
goto eqcon;
|
60 |
break;
|
61 |
case EQUIT:
|
62 |
xy.y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * P->ymf; |
63 |
eqcon:
|
64 |
xy.x = P->xmf * b * cosb * sinlam; |
65 |
break;
|
66 |
case N_POLE:
|
67 |
case S_POLE:
|
68 |
if (q >= 0.) { |
69 |
xy.x = (b = sqrt(q)) * sinlam; |
70 |
xy.y = coslam * (P->mode == S_POLE ? b : -b); |
71 |
} else
|
72 |
xy.x = xy.y = 0.;
|
73 |
break;
|
74 |
} |
75 |
return (xy);
|
76 |
} |
77 |
FORWARD(s_forward); /* spheroid */
|
78 |
double coslam, cosphi, sinphi;
|
79 |
|
80 |
sinphi = sin(lp.phi); |
81 |
cosphi = cos(lp.phi); |
82 |
coslam = cos(lp.lam); |
83 |
switch (P->mode) {
|
84 |
case EQUIT:
|
85 |
xy.y = 1. + cosphi * coslam;
|
86 |
goto oblcon;
|
87 |
case OBLIQ:
|
88 |
xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
|
89 |
oblcon:
|
90 |
if (xy.y <= EPS10) F_ERROR;
|
91 |
xy.x = (xy.y = sqrt(2. / xy.y)) * cosphi * sin(lp.lam);
|
92 |
xy.y *= P->mode == EQUIT ? sinphi : |
93 |
cosph0 * sinphi - sinph0 * cosphi * coslam; |
94 |
break;
|
95 |
case N_POLE:
|
96 |
coslam = -coslam; |
97 |
case S_POLE:
|
98 |
if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
|
99 |
xy.y = FORTPI - lp.phi * .5;
|
100 |
xy.y = 2. * (P->mode == S_POLE ? cos(xy.y) : sin(xy.y));
|
101 |
xy.x = xy.y * sin(lp.lam); |
102 |
xy.y *= coslam; |
103 |
break;
|
104 |
} |
105 |
return (xy);
|
106 |
} |
107 |
INVERSE(e_inverse); /* ellipsoid */
|
108 |
double cCe, sCe, q, rho, ab=0.0; |
109 |
|
110 |
switch (P->mode) {
|
111 |
case EQUIT:
|
112 |
case OBLIQ:
|
113 |
if ((rho = hypot(xy.x /= P->dd, xy.y *= P->dd)) < EPS10) {
|
114 |
lp.lam = 0.;
|
115 |
lp.phi = P->phi0; |
116 |
return (lp);
|
117 |
} |
118 |
cCe = cos(sCe = 2. * asin(.5 * rho / P->rq)); |
119 |
xy.x *= (sCe = sin(sCe)); |
120 |
if (P->mode == OBLIQ) {
|
121 |
q = P->qp * (ab = cCe * P->sinb1 + xy.y * sCe * P->cosb1 / rho); |
122 |
xy.y = rho * P->cosb1 * cCe - xy.y * P->sinb1 * sCe; |
123 |
} else {
|
124 |
q = P->qp * (ab = xy.y * sCe / rho); |
125 |
xy.y = rho * cCe; |
126 |
} |
127 |
break;
|
128 |
case N_POLE:
|
129 |
xy.y = -xy.y; |
130 |
case S_POLE:
|
131 |
if (!(q = (xy.x * xy.x + xy.y * xy.y)) ) {
|
132 |
lp.lam = 0.;
|
133 |
lp.phi = P->phi0; |
134 |
return (lp);
|
135 |
} |
136 |
/*
|
137 |
q = P->qp - q;
|
138 |
*/
|
139 |
ab = 1. - q / P->qp;
|
140 |
if (P->mode == S_POLE)
|
141 |
ab = - ab; |
142 |
break;
|
143 |
} |
144 |
lp.lam = atan2(xy.x, xy.y); |
145 |
lp.phi = pj_authlat(asin(ab), P->apa); |
146 |
return (lp);
|
147 |
} |
148 |
INVERSE(s_inverse); /* spheroid */
|
149 |
double cosz=0.0, rh, sinz=0.0; |
150 |
|
151 |
rh = hypot(xy.x, xy.y); |
152 |
if ((lp.phi = rh * .5 ) > 1.) I_ERROR; |
153 |
lp.phi = 2. * asin(lp.phi);
|
154 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
155 |
sinz = sin(lp.phi); |
156 |
cosz = cos(lp.phi); |
157 |
} |
158 |
switch (P->mode) {
|
159 |
case EQUIT:
|
160 |
lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
|
161 |
xy.x *= sinz; |
162 |
xy.y = cosz * rh; |
163 |
break;
|
164 |
case OBLIQ:
|
165 |
lp.phi = fabs(rh) <= EPS10 ? P->phi0 : |
166 |
asin(cosz * sinph0 + xy.y * sinz * cosph0 / rh); |
167 |
xy.x *= sinz * cosph0; |
168 |
xy.y = (cosz - sin(lp.phi) * sinph0) * rh; |
169 |
break;
|
170 |
case N_POLE:
|
171 |
xy.y = -xy.y; |
172 |
lp.phi = HALFPI - lp.phi; |
173 |
break;
|
174 |
case S_POLE:
|
175 |
lp.phi -= HALFPI; |
176 |
break;
|
177 |
} |
178 |
lp.lam = (xy.y == 0. && (P->mode == EQUIT || P->mode == OBLIQ)) ?
|
179 |
0. : atan2(xy.x, xy.y);
|
180 |
return (lp);
|
181 |
} |
182 |
FREEUP; if (P) pj_dalloc(P); }
|
183 |
ENTRY0(laea) |
184 |
double t;
|
185 |
|
186 |
if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
|
187 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
188 |
else if (fabs(t) < EPS10) |
189 |
P->mode = EQUIT; |
190 |
else
|
191 |
P->mode = OBLIQ; |
192 |
if (P->es) {
|
193 |
double sinphi;
|
194 |
|
195 |
P->e = sqrt(P->es); |
196 |
P->qp = pj_qsfn(1., P->e, P->one_es);
|
197 |
P->mmf = .5 / (1. - P->es); |
198 |
P->apa = pj_authset(P->es); |
199 |
switch (P->mode) {
|
200 |
case N_POLE:
|
201 |
case S_POLE:
|
202 |
P->dd = 1.;
|
203 |
break;
|
204 |
case EQUIT:
|
205 |
P->dd = 1. / (P->rq = sqrt(.5 * P->qp)); |
206 |
P->xmf = 1.;
|
207 |
P->ymf = .5 * P->qp;
|
208 |
break;
|
209 |
case OBLIQ:
|
210 |
P->rq = sqrt(.5 * P->qp);
|
211 |
sinphi = sin(P->phi0); |
212 |
P->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / P->qp; |
213 |
P->cosb1 = sqrt(1. - P->sinb1 * P->sinb1);
|
214 |
P->dd = cos(P->phi0) / (sqrt(1. - P->es * sinphi * sinphi) *
|
215 |
P->rq * P->cosb1); |
216 |
P->ymf = (P->xmf = P->rq) / P->dd; |
217 |
P->xmf *= P->dd; |
218 |
break;
|
219 |
} |
220 |
P->inv = e_inverse; |
221 |
P->fwd = e_forward; |
222 |
} else {
|
223 |
if (P->mode == OBLIQ) {
|
224 |
sinph0 = sin(P->phi0); |
225 |
cosph0 = cos(P->phi0); |
226 |
} |
227 |
P->inv = s_inverse; |
228 |
P->fwd = s_forward; |
229 |
} |
230 |
ENDENTRY(P) |