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