1 #define PJ_LIB__
2 #include <projects.h>
3
4 PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
5
6 struct pj_opaque {
7 double sinb1;
8 double cosb1;
9 double xmf;
10 double ymf;
11 double mmf;
12 double qp;
13 double dd;
14 double rq;
15 double *apa;
16 int mode;
17 };
18
19 #define EPS10 1.e-10
20 #define NITER 20
21 #define CONV 1.e-10
22 #define N_POLE 0
23 #define S_POLE 1
24 #define EQUIT 2
25 #define OBLIQ 3
26
e_forward(LP lp,PJ * P)27 static XY e_forward (LP lp, PJ *P) { /* Ellipsoidal, forward */
28 XY xy = {0.0,0.0};
29 struct pj_opaque *Q = P->opaque;
30 double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
31
32 coslam = cos(lp.lam);
33 sinlam = sin(lp.lam);
34 sinphi = sin(lp.phi);
35 q = pj_qsfn(sinphi, P->e, P->one_es);
36
37 if (Q->mode == OBLIQ || Q->mode == EQUIT) {
38 sinb = q / Q->qp;
39 cosb = sqrt(1. - sinb * sinb);
40 }
41
42 switch (Q->mode) {
43 case OBLIQ:
44 b = 1. + Q->sinb1 * sinb + Q->cosb1 * cosb * coslam;
45 break;
46 case EQUIT:
47 b = 1. + cosb * coslam;
48 break;
49 case N_POLE:
50 b = M_HALFPI + lp.phi;
51 q = Q->qp - q;
52 break;
53 case S_POLE:
54 b = lp.phi - M_HALFPI;
55 q = Q->qp + q;
56 break;
57 }
58 if (fabs(b) < EPS10) F_ERROR;
59
60 switch (Q->mode) {
61 case OBLIQ:
62 b = sqrt(2. / b);
63 xy.y = Q->ymf * b * (Q->cosb1 * sinb - Q->sinb1 * cosb * coslam);
64 goto eqcon;
65 break;
66 case EQUIT:
67 b = sqrt(2. / (1. + cosb * coslam));
68 xy.y = b * sinb * Q->ymf;
69 eqcon:
70 xy.x = Q->xmf * b * cosb * sinlam;
71 break;
72 case N_POLE:
73 case S_POLE:
74 if (q >= 0.) {
75 b = sqrt(q);
76 xy.x = b * sinlam;
77 xy.y = coslam * (Q->mode == S_POLE ? b : -b);
78 } else
79 xy.x = xy.y = 0.;
80 break;
81 }
82 return xy;
83 }
84
85
s_forward(LP lp,PJ * P)86 static XY s_forward (LP lp, PJ *P) { /* Spheroidal, forward */
87 XY xy = {0.0,0.0};
88 struct pj_opaque *Q = P->opaque;
89 double coslam, cosphi, sinphi;
90
91 sinphi = sin(lp.phi);
92 cosphi = cos(lp.phi);
93 coslam = cos(lp.lam);
94 switch (Q->mode) {
95 case EQUIT:
96 xy.y = 1. + cosphi * coslam;
97 goto oblcon;
98 case OBLIQ:
99 xy.y = 1. + Q->sinb1 * sinphi + Q->cosb1 * cosphi * coslam;
100 oblcon:
101 if (xy.y <= EPS10) F_ERROR;
102 xy.y = sqrt(2. / xy.y);
103 xy.x = xy.y * cosphi * sin(lp.lam);
104 xy.y *= Q->mode == EQUIT ? sinphi :
105 Q->cosb1 * sinphi - Q->sinb1 * cosphi * coslam;
106 break;
107 case N_POLE:
108 coslam = -coslam;
109 case S_POLE:
110 if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
111 xy.y = M_FORTPI - lp.phi * .5;
112 xy.y = 2. * (Q->mode == S_POLE ? cos(xy.y) : sin(xy.y));
113 xy.x = xy.y * sin(lp.lam);
114 xy.y *= coslam;
115 break;
116 }
117 return xy;
118 }
119
120
e_inverse(XY xy,PJ * P)121 static LP e_inverse (XY xy, PJ *P) { /* Ellipsoidal, inverse */
122 LP lp = {0.0,0.0};
123 struct pj_opaque *Q = P->opaque;
124 double cCe, sCe, q, rho, ab=0.0;
125
126 switch (Q->mode) {
127 case EQUIT:
128 case OBLIQ:
129 xy.x /= Q->dd;
130 xy.y *= Q->dd;
131 rho = hypot(xy.x, xy.y);
132 if (rho < EPS10) {
133 lp.lam = 0.;
134 lp.phi = P->phi0;
135 return lp;
136 }
137 sCe = 2. * asin(.5 * rho / Q->rq);
138 cCe = cos(sCe);
139 sCe = sin(sCe);
140 xy.x *= sCe;
141 if (Q->mode == OBLIQ) {
142 ab = cCe * Q->sinb1 + xy.y * sCe * Q->cosb1 / rho;
143 xy.y = rho * Q->cosb1 * cCe - xy.y * Q->sinb1 * sCe;
144 } else {
145 ab = xy.y * sCe / rho;
146 xy.y = rho * cCe;
147 }
148 break;
149 case N_POLE:
150 xy.y = -xy.y;
151 case S_POLE:
152 q = (xy.x * xy.x + xy.y * xy.y);
153 if (!q) {
154 lp.lam = 0.;
155 lp.phi = P->phi0;
156 return (lp);
157 }
158 ab = 1. - q / Q->qp;
159 if (Q->mode == S_POLE)
160 ab = - ab;
161 break;
162 }
163 lp.lam = atan2(xy.x, xy.y);
164 lp.phi = pj_authlat(asin(ab), Q->apa);
165 return lp;
166 }
167
168
s_inverse(XY xy,PJ * P)169 static LP s_inverse (XY xy, PJ *P) { /* Spheroidal, inverse */
170 LP lp = {0.0,0.0};
171 struct pj_opaque *Q = P->opaque;
172 double cosz=0.0, rh, sinz=0.0;
173
174 rh = hypot(xy.x, xy.y);
175 if ((lp.phi = rh * .5 ) > 1.) I_ERROR;
176 lp.phi = 2. * asin(lp.phi);
177 if (Q->mode == OBLIQ || Q->mode == EQUIT) {
178 sinz = sin(lp.phi);
179 cosz = cos(lp.phi);
180 }
181 switch (Q->mode) {
182 case EQUIT:
183 lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
184 xy.x *= sinz;
185 xy.y = cosz * rh;
186 break;
187 case OBLIQ:
188 lp.phi = fabs(rh) <= EPS10 ? P->phi0 :
189 asin(cosz * Q->sinb1 + xy.y * sinz * Q->cosb1 / rh);
190 xy.x *= sinz * Q->cosb1;
191 xy.y = (cosz - sin(lp.phi) * Q->sinb1) * rh;
192 break;
193 case N_POLE:
194 xy.y = -xy.y;
195 lp.phi = M_HALFPI - lp.phi;
196 break;
197 case S_POLE:
198 lp.phi -= M_HALFPI;
199 break;
200 }
201 lp.lam = (xy.y == 0. && (Q->mode == EQUIT || Q->mode == OBLIQ)) ?
202 0. : atan2(xy.x, xy.y);
203 return (lp);
204 }
205
206
freeup_new(PJ * P)207 static void *freeup_new (PJ *P) { /* Destructor */
208 if (0==P)
209 return 0;
210 if (0==P->opaque)
211 return pj_dealloc (P);
212
213 pj_dealloc (P->opaque->apa);
214 pj_dealloc (P->opaque);
215 return pj_dealloc(P);
216 }
217
freeup(PJ * P)218 static void freeup (PJ *P) {
219 freeup_new (P);
220 return;
221 }
222
223
PROJECTION(laea)224 PJ *PROJECTION(laea) {
225 double t;
226 struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
227 if (0==Q)
228 return freeup_new (P);
229 P->opaque = Q;
230
231 t = fabs(P->phi0);
232 if (fabs(t - M_HALFPI) < EPS10)
233 Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
234 else if (fabs(t) < EPS10)
235 Q->mode = EQUIT;
236 else
237 Q->mode = OBLIQ;
238 if (P->es) {
239 double sinphi;
240
241 P->e = sqrt(P->es);
242 Q->qp = pj_qsfn(1., P->e, P->one_es);
243 Q->mmf = .5 / (1. - P->es);
244 Q->apa = pj_authset(P->es);
245 switch (Q->mode) {
246 case N_POLE:
247 case S_POLE:
248 Q->dd = 1.;
249 break;
250 case EQUIT:
251 Q->dd = 1. / (Q->rq = sqrt(.5 * Q->qp));
252 Q->xmf = 1.;
253 Q->ymf = .5 * Q->qp;
254 break;
255 case OBLIQ:
256 Q->rq = sqrt(.5 * Q->qp);
257 sinphi = sin(P->phi0);
258 Q->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / Q->qp;
259 Q->cosb1 = sqrt(1. - Q->sinb1 * Q->sinb1);
260 Q->dd = cos(P->phi0) / (sqrt(1. - P->es * sinphi * sinphi) *
261 Q->rq * Q->cosb1);
262 Q->ymf = (Q->xmf = Q->rq) / Q->dd;
263 Q->xmf *= Q->dd;
264 break;
265 }
266 P->inv = e_inverse;
267 P->fwd = e_forward;
268 } else {
269 if (Q->mode == OBLIQ) {
270 Q->sinb1 = sin(P->phi0);
271 Q->cosb1 = cos(P->phi0);
272 }
273 P->inv = s_inverse;
274 P->fwd = s_forward;
275 }
276
277 return P;
278 }
279
280
281 #ifndef PJ_SELFTEST
pj_laea_selftest(void)282 int pj_laea_selftest (void) {return 0;}
283 #else
284
pj_laea_selftest(void)285 int pj_laea_selftest (void) {
286 double tolerance_lp = 1e-10;
287 double tolerance_xy = 1e-7;
288
289 char e_args[] = {"+proj=laea +ellps=GRS80 +lat_1=0.5 +lat_2=2"};
290 char s_args[] = {"+proj=laea +a=6400000 +lat_1=0.5 +lat_2=2"};
291
292 LP fwd_in[] = {
293 { 2, 1},
294 { 2,-1},
295 {-2, 1},
296 {-2,-1}
297 };
298
299 XY e_fwd_expect[] = {
300 { 222602.471450095181, 110589.82722441027},
301 { 222602.471450095181, -110589.827224408786},
302 {-222602.471450095181, 110589.82722441027},
303 {-222602.471450095181, -110589.827224408786},
304 };
305
306 XY s_fwd_expect[] = {
307 { 223365.281370124663, 111716.668072915665},
308 { 223365.281370124663, -111716.668072915665},
309 {-223365.281370124663, 111716.668072915665},
310 {-223365.281370124663, -111716.668072915665},
311 };
312
313 XY inv_in[] = {
314 { 200, 100},
315 { 200,-100},
316 {-200, 100},
317 {-200,-100}
318 };
319
320 LP e_inv_expect[] = {
321 { 0.00179663056847900867, 0.000904369475966495845},
322 { 0.00179663056847900867, -0.000904369475966495845},
323 {-0.00179663056847900867, 0.000904369475966495845},
324 {-0.00179663056847900867, -0.000904369475966495845},
325 };
326
327 LP s_inv_expect[] = {
328 { 0.00179049311002060264, 0.000895246554791735271},
329 { 0.00179049311002060264, -0.000895246554791735271},
330 {-0.00179049311002060264, 0.000895246554791735271},
331 {-0.00179049311002060264, -0.000895246554791735271},
332 };
333
334 return pj_generic_selftest (e_args, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, s_fwd_expect, inv_in, e_inv_expect, s_inv_expect);
335 }
336
337
338 #endif
339