gvsig-3d / 1.10 / trunk / libraries / libjni-proj4 / src / PJ_aeqd.c @ 72
History | View | Annotate | Download (7.45 KB)
1 | 5 | jzarzoso | /******************************************************************************
|
---|---|---|---|
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) |