svn-gvsig-desktop / trunk / libraries / libjni-proj4 / src / PJ_aeqd.c @ 7098
History | View | Annotate | Download (7.45 KB)
1 |
/******************************************************************************
|
---|---|
2 |
* $Id: PJ_aeqd.c,v 1.3 2002/12/14 19:27:06 warmerda Exp $
|
3 |
*
|
4 |
* Project: PROJ.4
|
5 |
* Purpose: Implementation of the aeqd (Azimuthal Equidistant) projection.
|
6 |
* Author: Gerald Evenden
|
7 |
*
|
8 |
******************************************************************************
|
9 |
* Copyright (c) 1995, Gerald Evenden
|
10 |
*
|
11 |
* Permission is hereby granted, free of charge, to any person obtaining a
|
12 |
* copy of this software and associated documentation files (the "Software"),
|
13 |
* to deal in the Software without restriction, including without limitation
|
14 |
* the rights to use, copy, modify, merge, publish, distribute, sublicense,
|
15 |
* and/or sell copies of the Software, and to permit persons to whom the
|
16 |
* Software is furnished to do so, subject to the following conditions:
|
17 |
*
|
18 |
* The above copyright notice and this permission notice shall be included
|
19 |
* in all copies or substantial portions of the Software.
|
20 |
*
|
21 |
* THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
|
22 |
* OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
|
23 |
* FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
|
24 |
* THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
|
25 |
* LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
|
26 |
* FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
|
27 |
* DEALINGS IN THE SOFTWARE.
|
28 |
******************************************************************************
|
29 |
*
|
30 |
* $Log: PJ_aeqd.c,v $
|
31 |
* Revision 1.3 2002/12/14 19:27:06 warmerda
|
32 |
* updated header
|
33 |
*
|
34 |
*/
|
35 |
|
36 |
#define PROJ_PARMS__ \
|
37 |
double sinph0; \
|
38 |
double cosph0; \
|
39 |
double *en; \
|
40 |
double M1; \
|
41 |
double N1; \
|
42 |
double Mp; \
|
43 |
double He; \
|
44 |
double G; \
|
45 |
int mode;
|
46 |
#define PJ_LIB__
|
47 |
#include <projects.h> |
48 |
|
49 |
PJ_CVSID("$Id: PJ_aeqd.c,v 1.3 2002/12/14 19:27:06 warmerda Exp $");
|
50 |
|
51 |
PROJ_HEAD(aeqd, "Azimuthal Equidistant") "\n\tAzi, Sph&Ell\n\tlat_0 guam"; |
52 |
|
53 |
#define EPS10 1.e-10 |
54 |
#define TOL 1.e-14 |
55 |
|
56 |
#define N_POLE 0 |
57 |
#define S_POLE 1 |
58 |
#define EQUIT 2 |
59 |
#define OBLIQ 3 |
60 |
FORWARD(e_guam_fwd); /* Guam elliptical */
|
61 |
double cosphi, sinphi, t;
|
62 |
|
63 |
cosphi = cos(lp.phi); |
64 |
sinphi = sin(lp.phi); |
65 |
t = 1. / sqrt(1. - P->es * sinphi * sinphi); |
66 |
xy.x = lp.lam * cosphi * t; |
67 |
xy.y = pj_mlfn(lp.phi, sinphi, cosphi, P->en) - P->M1 + |
68 |
.5 * lp.lam * lp.lam * cosphi * sinphi * t;
|
69 |
return (xy);
|
70 |
} |
71 |
FORWARD(e_forward); /* elliptical */
|
72 |
double coslam, cosphi, sinphi, rho, s, H, H2, c, Az, t, ct, st, cA, sA;
|
73 |
|
74 |
coslam = cos(lp.lam); |
75 |
cosphi = cos(lp.phi); |
76 |
sinphi = sin(lp.phi); |
77 |
switch (P->mode) {
|
78 |
case N_POLE:
|
79 |
coslam = - coslam; |
80 |
case S_POLE:
|
81 |
xy.x = (rho = fabs(P->Mp - pj_mlfn(lp.phi, sinphi, cosphi, P->en))) * |
82 |
sin(lp.lam); |
83 |
xy.y = rho * coslam; |
84 |
break;
|
85 |
case EQUIT:
|
86 |
case OBLIQ:
|
87 |
if (fabs(lp.lam) < EPS10 && fabs(lp.phi - P->phi0) < EPS10) {
|
88 |
xy.x = xy.y = 0.;
|
89 |
break;
|
90 |
} |
91 |
t = atan2(P->one_es * sinphi + P->es * P->N1 * P->sinph0 * |
92 |
sqrt(1. - P->es * sinphi * sinphi), cosphi);
|
93 |
ct = cos(t); st = sin(t); |
94 |
Az = atan2(sin(lp.lam) * ct, P->cosph0 * st - P->sinph0 * coslam * ct); |
95 |
cA = cos(Az); sA = sin(Az); |
96 |
s = aasin( fabs(sA) < TOL ? |
97 |
(P->cosph0 * st - P->sinph0 * coslam * ct) / cA : |
98 |
sin(lp.lam) * ct / sA ); |
99 |
H = P->He * cA; |
100 |
H2 = H * H; |
101 |
c = P->N1 * s * (1. + s * s * (- H2 * (1. - H2)/6. + |
102 |
s * ( P->G * H * (1. - 2. * H2 * H2) / 8. + |
103 |
s * ((H2 * (4. - 7. * H2) - 3. * P->G * P->G * (1. - 7. * H2)) / |
104 |
120. - s * P->G * H / 48.)))); |
105 |
xy.x = c * sA; |
106 |
xy.y = c * cA; |
107 |
break;
|
108 |
} |
109 |
return (xy);
|
110 |
} |
111 |
FORWARD(s_forward); /* spherical */
|
112 |
double coslam, cosphi, sinphi;
|
113 |
|
114 |
sinphi = sin(lp.phi); |
115 |
cosphi = cos(lp.phi); |
116 |
coslam = cos(lp.lam); |
117 |
switch (P->mode) {
|
118 |
case EQUIT:
|
119 |
xy.y = cosphi * coslam; |
120 |
goto oblcon;
|
121 |
case OBLIQ:
|
122 |
xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam; |
123 |
oblcon:
|
124 |
if (fabs(fabs(xy.y) - 1.) < TOL) |
125 |
if (xy.y < 0.) |
126 |
F_ERROR |
127 |
else
|
128 |
xy.x = xy.y = 0.;
|
129 |
else {
|
130 |
xy.y = acos(xy.y); |
131 |
xy.y /= sin(xy.y); |
132 |
xy.x = xy.y * cosphi * sin(lp.lam); |
133 |
xy.y *= (P->mode == EQUIT) ? sinphi : |
134 |
P->cosph0 * sinphi - P->sinph0 * cosphi * coslam; |
135 |
} |
136 |
break;
|
137 |
case N_POLE:
|
138 |
lp.phi = -lp.phi; |
139 |
coslam = -coslam; |
140 |
case S_POLE:
|
141 |
if (fabs(lp.phi - HALFPI) < EPS10) F_ERROR;
|
142 |
xy.x = (xy.y = (HALFPI + lp.phi)) * sin(lp.lam); |
143 |
xy.y *= coslam; |
144 |
break;
|
145 |
} |
146 |
return (xy);
|
147 |
} |
148 |
INVERSE(e_guam_inv); /* Guam elliptical */
|
149 |
double x2, t;
|
150 |
int i;
|
151 |
|
152 |
x2 = 0.5 * xy.x * xy.x; |
153 |
lp.phi = P->phi0; |
154 |
for (i = 0; i < 3; ++i) { |
155 |
t = P->e * sin(lp.phi); |
156 |
lp.phi = pj_inv_mlfn(P->M1 + xy.y - |
157 |
x2 * tan(lp.phi) * (t = sqrt(1. - t * t)), P->es, P->en);
|
158 |
} |
159 |
lp.lam = xy.x * t / cos(lp.phi); |
160 |
return (lp);
|
161 |
} |
162 |
INVERSE(e_inverse); /* elliptical */
|
163 |
double c, Az, cosAz, A, B, D, E, F, psi, t;
|
164 |
|
165 |
if ((c = hypot(xy.x, xy.y)) < EPS10) {
|
166 |
lp.phi = P->phi0; |
167 |
lp.lam = 0.;
|
168 |
return (lp);
|
169 |
} |
170 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
171 |
cosAz = cos(Az = atan2(xy.x, xy.y)); |
172 |
t = P->cosph0 * cosAz; |
173 |
B = P->es * t / P->one_es; |
174 |
A = - B * t; |
175 |
B *= 3. * (1. - A) * P->sinph0; |
176 |
D = c / P->N1; |
177 |
E = D * (1. - D * D * (A * (1. + A) / 6. + B * (1. + 3.*A) * D / 24.)); |
178 |
F = 1. - E * E * (A / 2. + B * E / 6.); |
179 |
psi = aasin(P->sinph0 * cos(E) + t * sin(E)); |
180 |
lp.lam = aasin(sin(Az) * sin(E) / cos(psi)); |
181 |
if ((t = fabs(psi)) < EPS10)
|
182 |
lp.phi = 0.;
|
183 |
else if (fabs(t - HALFPI) < 0.) |
184 |
lp.phi = HALFPI; |
185 |
else
|
186 |
lp.phi = atan((1. - P->es * F * P->sinph0 / sin(psi)) * tan(psi) /
|
187 |
P->one_es); |
188 |
} else { /* Polar */ |
189 |
lp.phi = pj_inv_mlfn(P->mode == N_POLE ? P->Mp - c : P->Mp + c, |
190 |
P->es, P->en); |
191 |
lp.lam = atan2(xy.x, P->mode == N_POLE ? -xy.y : xy.y); |
192 |
} |
193 |
return (lp);
|
194 |
} |
195 |
INVERSE(s_inverse); /* spherical */
|
196 |
double cosc, c_rh, sinc;
|
197 |
|
198 |
if ((c_rh = hypot(xy.x, xy.y)) > PI) {
|
199 |
if (c_rh - EPS10 > PI) I_ERROR;
|
200 |
c_rh = PI; |
201 |
} else if (c_rh < EPS10) { |
202 |
lp.phi = P->phi0; |
203 |
lp.lam = 0.;
|
204 |
return (lp);
|
205 |
} |
206 |
if (P->mode == OBLIQ || P->mode == EQUIT) {
|
207 |
sinc = sin(c_rh); |
208 |
cosc = cos(c_rh); |
209 |
if (P->mode == EQUIT) {
|
210 |
lp.phi = aasin(xy.y * sinc / c_rh); |
211 |
xy.x *= sinc; |
212 |
xy.y = cosc * c_rh; |
213 |
} else {
|
214 |
lp.phi = aasin(cosc * P->sinph0 + xy.y * sinc * P->cosph0 / |
215 |
c_rh); |
216 |
xy.y = (cosc - P->sinph0 * sin(lp.phi)) * c_rh; |
217 |
xy.x *= sinc * P->cosph0; |
218 |
} |
219 |
lp.lam = xy.y == 0. ? 0. : atan2(xy.x, xy.y); |
220 |
} else if (P->mode == N_POLE) { |
221 |
lp.phi = HALFPI - c_rh; |
222 |
lp.lam = atan2(xy.x, -xy.y); |
223 |
} else {
|
224 |
lp.phi = c_rh - HALFPI; |
225 |
lp.lam = atan2(xy.x, xy.y); |
226 |
} |
227 |
return (lp);
|
228 |
} |
229 |
FREEUP; |
230 |
if (P) {
|
231 |
if (P->en)
|
232 |
pj_dalloc(P->en); |
233 |
pj_dalloc(P); |
234 |
} |
235 |
} |
236 |
ENTRY1(aeqd, en) |
237 |
P->phi0 = pj_param(P->params, "rlat_0").f;
|
238 |
if (fabs(fabs(P->phi0) - HALFPI) < EPS10) {
|
239 |
P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
|
240 |
P->sinph0 = P->phi0 < 0. ? -1. : 1.; |
241 |
P->cosph0 = 0.;
|
242 |
} else if (fabs(P->phi0) < EPS10) { |
243 |
P->mode = EQUIT; |
244 |
P->sinph0 = 0.;
|
245 |
P->cosph0 = 1.;
|
246 |
} else {
|
247 |
P->mode = OBLIQ; |
248 |
P->sinph0 = sin(P->phi0); |
249 |
P->cosph0 = cos(P->phi0); |
250 |
} |
251 |
if (! P->es) {
|
252 |
P->inv = s_inverse; P->fwd = s_forward; |
253 |
} else {
|
254 |
if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
|
255 |
if (pj_param(P->params, "bguam").i) { |
256 |
P->M1 = pj_mlfn(P->phi0, P->sinph0, P->cosph0, P->en); |
257 |
P->inv = e_guam_inv; P->fwd = e_guam_fwd; |
258 |
} else {
|
259 |
switch (P->mode) {
|
260 |
case N_POLE:
|
261 |
P->Mp = pj_mlfn(HALFPI, 1., 0., P->en); |
262 |
break;
|
263 |
case S_POLE:
|
264 |
P->Mp = pj_mlfn(-HALFPI, -1., 0., P->en); |
265 |
break;
|
266 |
case EQUIT:
|
267 |
case OBLIQ:
|
268 |
P->inv = e_inverse; P->fwd = e_forward; |
269 |
P->N1 = 1. / sqrt(1. - P->es * P->sinph0 * P->sinph0); |
270 |
P->G = P->sinph0 * (P->He = P->e / sqrt(P->one_es)); |
271 |
P->He *= P->cosph0; |
272 |
break;
|
273 |
} |
274 |
P->inv = e_inverse; P->fwd = e_forward; |
275 |
} |
276 |
} |
277 |
ENDENTRY(P) |