/src/gdal/proj/src/projections/laea.cpp

Source (jump to first uncovered line)
#include "proj.h"
#include "proj_internal.h"
#include <errno.h>
#include <math.h>

PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";

namespace pj_laea_ns {
enum Mode { N_POLE = 0, S_POLE = 1, EQUIT = 2, OBLIQ = 3 };
}

namespace { // anonymous namespace
struct pj_laea_data {
    double sinb1;
    double cosb1;
    double xmf;
    double ymf;
    double mmf;
    double qp;
    double dd;
    double rq;
    double *apa;
    enum pj_laea_ns::Mode mode;
};
} // anonymous namespace

#define EPS10 1.e-10

static PJ_XY laea_e_forward(PJ_LP lp, PJ *P) { /* Ellipsoidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_laea_data *Q = static_cast<struct pj_laea_data *>(P->opaque);
    double coslam, sinlam, sinphi, cosphi, q, sinb = 0.0, cosb = 0.0, b = 0.0;

    coslam = cos(lp.lam);
    sinlam = sin(lp.lam);
    sinphi = sin(lp.phi);
    cosphi = cos(lp.phi);
    const double xi = pj_authalic_lat(lp.phi, sinphi, cosphi,
                                      Q->apa, P, Q->qp);
    q = sin(xi) * Q->qp;

    if (Q->mode == pj_laea_ns::OBLIQ || Q->mode == pj_laea_ns::EQUIT) {
        sinb = sin(xi); cosb = cos(xi);
    }

    switch (Q->mode) {
    case pj_laea_ns::OBLIQ:
        b = 1. + Q->sinb1 * sinb + Q->cosb1 * cosb * coslam;
        break;
    case pj_laea_ns::EQUIT:
        b = 1. + cosb * coslam;
        break;
    case pj_laea_ns::N_POLE:
        b = M_HALFPI + lp.phi;
        q = Q->qp - q;
        break;
    case pj_laea_ns::S_POLE:
        b = lp.phi - M_HALFPI;
        q = Q->qp + q;
        break;
    }
    if (fabs(b) < EPS10) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return xy;
    }

    switch (Q->mode) {
    case pj_laea_ns::OBLIQ:
        b = sqrt(2. / b);
        xy.y = Q->ymf * b * (Q->cosb1 * sinb - Q->sinb1 * cosb * coslam);
        goto eqcon;
        break;
    case pj_laea_ns::EQUIT:
        b = sqrt(2. / (1. + cosb * coslam));
        xy.y = b * sinb * Q->ymf;
    eqcon:
        xy.x = Q->xmf * b * cosb * sinlam;
        break;
    case pj_laea_ns::N_POLE:
    case pj_laea_ns::S_POLE:
        if (q >= 1e-15) {
            b = sqrt(q);
            xy.x = b * sinlam;
            xy.y = coslam * (Q->mode == pj_laea_ns::S_POLE ? b : -b);
        } else
            xy.x = xy.y = 0.;
        break;
    }
    return xy;
}

static PJ_XY laea_s_forward(PJ_LP lp, PJ *P) { /* Spheroidal, forward */
    PJ_XY xy = {0.0, 0.0};
    struct pj_laea_data *Q = static_cast<struct pj_laea_data *>(P->opaque);
    double coslam, cosphi, sinphi;

    sinphi = sin(lp.phi);
    cosphi = cos(lp.phi);
    coslam = cos(lp.lam);
    switch (Q->mode) {
    case pj_laea_ns::EQUIT:
        xy.y = 1. + cosphi * coslam;
        goto oblcon;
    case pj_laea_ns::OBLIQ:
        xy.y = 1. + Q->sinb1 * sinphi + Q->cosb1 * cosphi * coslam;
    oblcon:
        if (xy.y <= EPS10) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
        xy.y = sqrt(2. / xy.y);
        xy.x = xy.y * cosphi * sin(lp.lam);
        xy.y *= Q->mode == pj_laea_ns::EQUIT
                    ? sinphi
                    : Q->cosb1 * sinphi - Q->sinb1 * cosphi * coslam;
        break;
    case pj_laea_ns::N_POLE:
        coslam = -coslam;
        PROJ_FALLTHROUGH;
    case pj_laea_ns::S_POLE:
        if (fabs(lp.phi + P->phi0) < EPS10) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return xy;
        }
        xy.y = M_FORTPI - lp.phi * .5;
        xy.y = 2. * (Q->mode == pj_laea_ns::S_POLE ? cos(xy.y) : sin(xy.y));
        xy.x = xy.y * sin(lp.lam);
        xy.y *= coslam;
        break;
    }
    return xy;
}

static PJ_LP laea_e_inverse(PJ_XY xy, PJ *P) { /* Ellipsoidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_laea_data *Q = static_cast<struct pj_laea_data *>(P->opaque);
    double cCe, sCe, q, rho, ab = 0.0;

    switch (Q->mode) {
    case pj_laea_ns::EQUIT:
    case pj_laea_ns::OBLIQ: {
        xy.x /= Q->dd;
        xy.y *= Q->dd;
        rho = hypot(xy.x, xy.y);
        if (rho < EPS10) {
            lp.lam = 0.;
            lp.phi = P->phi0;
            return lp;
        }
        const double asin_argument = .5 * rho / Q->rq;
        if (asin_argument > 1) {
            proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
            return lp;
        }
        sCe = 2. * asin(asin_argument);
        cCe = cos(sCe);
        sCe = sin(sCe);
        xy.x *= sCe;
        if (Q->mode == pj_laea_ns::OBLIQ) {
            ab = cCe * Q->sinb1 + xy.y * sCe * Q->cosb1 / rho;
            xy.y = rho * Q->cosb1 * cCe - xy.y * Q->sinb1 * sCe;
        } else {
            ab = xy.y * sCe / rho;
            xy.y = rho * cCe;
        }
        break;
    }
    case pj_laea_ns::N_POLE:
        xy.y = -xy.y;
        PROJ_FALLTHROUGH;
    case pj_laea_ns::S_POLE:
        q = (xy.x * xy.x + xy.y * xy.y);
        if (q == 0.0) {
            lp.lam = 0.;
            lp.phi = P->phi0;
            return (lp);
        }
        ab = 1. - q / Q->qp;
        if (Q->mode == pj_laea_ns::S_POLE)
            ab = -ab;
        break;
    }
    lp.lam = atan2(xy.x, xy.y);
    lp.phi = pj_authalic_lat_inverse(asin(ab), Q->apa, P, Q->qp);
    return lp;
}

static PJ_LP laea_s_inverse(PJ_XY xy, PJ *P) { /* Spheroidal, inverse */
    PJ_LP lp = {0.0, 0.0};
    struct pj_laea_data *Q = static_cast<struct pj_laea_data *>(P->opaque);
    double cosz = 0.0, rh, sinz = 0.0;

    rh = hypot(xy.x, xy.y);
    if ((lp.phi = rh * .5) > 1.) {
        proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
        return lp;
    }
    lp.phi = 2. * asin(lp.phi);
    if (Q->mode == pj_laea_ns::OBLIQ || Q->mode == pj_laea_ns::EQUIT) {
        sinz = sin(lp.phi);
        cosz = cos(lp.phi);
    }
    switch (Q->mode) {
    case pj_laea_ns::EQUIT:
        lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
        xy.x *= sinz;
        xy.y = cosz * rh;
        break;
    case pj_laea_ns::OBLIQ:
        lp.phi = fabs(rh) <= EPS10
                     ? P->phi0
                     : asin(cosz * Q->sinb1 + xy.y * sinz * Q->cosb1 / rh);
        xy.x *= sinz * Q->cosb1;
        xy.y = (cosz - sin(lp.phi) * Q->sinb1) * rh;
        break;
    case pj_laea_ns::N_POLE:
        xy.y = -xy.y;
        lp.phi = M_HALFPI - lp.phi;
        break;
    case pj_laea_ns::S_POLE:
        lp.phi -= M_HALFPI;
        break;
    }
    lp.lam = (xy.y == 0. &&
              (Q->mode == pj_laea_ns::EQUIT || Q->mode == pj_laea_ns::OBLIQ))
                 ? 0.
                 : atan2(xy.x, xy.y);
    return (lp);
}

static PJ *pj_laea_destructor(PJ *P, int errlev) {
    if (nullptr == P)
        return nullptr;

    if (nullptr == P->opaque)
        return pj_default_destructor(P, errlev);

    free(static_cast<struct pj_laea_data *>(P->opaque)->apa);

    return pj_default_destructor(P, errlev);
}

PJ *PJ_PROJECTION(laea) {
    double t;
    struct pj_laea_data *Q = static_cast<struct pj_laea_data *>(
        calloc(1, sizeof(struct pj_laea_data)));
    if (nullptr == Q)
        return pj_default_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
    P->opaque = Q;
    P->destructor = pj_laea_destructor;

    t = fabs(P->phi0);
    if (t > M_HALFPI + EPS10) {
        proj_log_error(P,
                       _("Invalid value for lat_0: |lat_0| should be <= 90°"));
        return pj_laea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
    }
    if (fabs(t - M_HALFPI) < EPS10)
        Q->mode = P->phi0 < 0. ? pj_laea_ns::S_POLE : pj_laea_ns::N_POLE;
    else if (fabs(t) < EPS10)
        Q->mode = pj_laea_ns::EQUIT;
    else
        Q->mode = pj_laea_ns::OBLIQ;
    if (P->es != 0.0) {
        double sinphi, cosphi;

        P->e = sqrt(P->es);
        Q->qp = pj_authalic_lat_q(1.0, P);
        Q->mmf = .5 / (1. - P->es);
        Q->apa = pj_authalic_lat_compute_coeffs(P->n);
        if (nullptr == Q->apa)
            return pj_laea_destructor(P, PROJ_ERR_OTHER /*ENOMEM*/);
        switch (Q->mode) {
        case pj_laea_ns::N_POLE:
        case pj_laea_ns::S_POLE:
            Q->dd = 1.;
            break;
        case pj_laea_ns::EQUIT:
            Q->dd = 1. / (Q->rq = sqrt(.5 * Q->qp));
            Q->xmf = 1.;
            Q->ymf = .5 * Q->qp;
            break;
        case pj_laea_ns::OBLIQ:
            Q->rq = sqrt(.5 * Q->qp);
            sinphi = sin(P->phi0); cosphi = cos(P->phi0);
            const double b1 = pj_authalic_lat(P->phi0, sinphi, cosphi,
                                              Q->apa, P, Q->qp);
            Q->sinb1 = sin(b1);
            Q->cosb1 = cos(b1);
            Q->dd = cos(P->phi0) /
                    (sqrt(1. - P->es * sinphi * sinphi) * Q->rq * Q->cosb1);
            Q->ymf = (Q->xmf = Q->rq) / Q->dd;
            Q->xmf *= Q->dd;
            break;
        }
        P->inv = laea_e_inverse;
        P->fwd = laea_e_forward;
    } else {
        if (Q->mode == pj_laea_ns::OBLIQ) {
            Q->sinb1 = sin(P->phi0);
            Q->cosb1 = cos(P->phi0);
        }
        P->inv = laea_s_inverse;
        P->fwd = laea_s_forward;
    }

    return P;
}

Coverage Report

Created: 2025-06-09 08:44

Line	Count	Source (jump to first uncovered line)
1		#include "proj.h"
2		#include "proj_internal.h"
3		#include <errno.h>
4		#include <math.h>
5
6		PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
7
8		namespace pj_laea_ns {
9		enum Mode { N_POLE = 0, S_POLE = 1, EQUIT = 2, OBLIQ = 3 };
10		}
11
12		namespace { // anonymous namespace
13		struct pj_laea_data {
14		double sinb1;
15		double cosb1;
16		double xmf;
17		double ymf;
18		double mmf;
19		double qp;
20		double dd;
21		double rq;
22		double *apa;
23		enum pj_laea_ns::Mode mode;
24		};
25		} // anonymous namespace
26
27	109k	#define EPS10 1.e-10
28
29	0	static PJ_XY laea_e_forward(PJ_LP lp, PJ P) { / Ellipsoidal, forward */
30	0	PJ_XY xy = {0.0, 0.0};
31	0	struct pj_laea_data Q = static_cast<struct pj_laea_data >(P->opaque);
32	0	double coslam, sinlam, sinphi, cosphi, q, sinb = 0.0, cosb = 0.0, b = 0.0;
33
34	0	coslam = cos(lp.lam);
35	0	sinlam = sin(lp.lam);
36	0	sinphi = sin(lp.phi);
37	0	cosphi = cos(lp.phi);
38	0	const double xi = pj_authalic_lat(lp.phi, sinphi, cosphi,
39	0	Q->apa, P, Q->qp);
40	0	q = sin(xi) * Q->qp;
41
42	0	if (Q->mode == pj_laea_ns::OBLIQ \|\| Q->mode == pj_laea_ns::EQUIT) {
43	0	sinb = sin(xi); cosb = cos(xi);
44	0	}
45
46	0	switch (Q->mode) {
47	0	case pj_laea_ns::OBLIQ:
48	0	b = 1. + Q->sinb1 * sinb + Q->cosb1 * cosb * coslam;
49	0	break;
50	0	case pj_laea_ns::EQUIT:
51	0	b = 1. + cosb * coslam;
52	0	break;
53	0	case pj_laea_ns::N_POLE:
54	0	b = M_HALFPI + lp.phi;
55	0	q = Q->qp - q;
56	0	break;
57	0	case pj_laea_ns::S_POLE:
58	0	b = lp.phi - M_HALFPI;
59	0	q = Q->qp + q;
60	0	break;
61	0	}
62	0	if (fabs(b) < EPS10) {
63	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
64	0	return xy;
65	0	}
66
67	0	switch (Q->mode) {
68	0	case pj_laea_ns::OBLIQ:
69	0	b = sqrt(2. / b);
70	0	xy.y = Q->ymf * b * (Q->cosb1 * sinb - Q->sinb1 * cosb * coslam);
71	0	goto eqcon;
72	0	break;
73	0	case pj_laea_ns::EQUIT:
74	0	b = sqrt(2. / (1. + cosb * coslam));
75	0	xy.y = b * sinb * Q->ymf;
76	0	eqcon:
77	0	xy.x = Q->xmf * b * cosb * sinlam;
78	0	break;
79	0	case pj_laea_ns::N_POLE:
80	0	case pj_laea_ns::S_POLE:
81	0	if (q >= 1e-15) {
82	0	b = sqrt(q);
83	0	xy.x = b * sinlam;
84	0	xy.y = coslam * (Q->mode == pj_laea_ns::S_POLE ? b : -b);
85	0	} else
86	0	xy.x = xy.y = 0.;
87	0	break;
88	0	}
89	0	return xy;
90	0	}
91
92	0	static PJ_XY laea_s_forward(PJ_LP lp, PJ P) { / Spheroidal, forward */
93	0	PJ_XY xy = {0.0, 0.0};
94	0	struct pj_laea_data Q = static_cast<struct pj_laea_data >(P->opaque);
95	0	double coslam, cosphi, sinphi;
96
97	0	sinphi = sin(lp.phi);
98	0	cosphi = cos(lp.phi);
99	0	coslam = cos(lp.lam);
100	0	switch (Q->mode) {
101	0	case pj_laea_ns::EQUIT:
102	0	xy.y = 1. + cosphi * coslam;
103	0	goto oblcon;
104	0	case pj_laea_ns::OBLIQ:
105	0	xy.y = 1. + Q->sinb1 * sinphi + Q->cosb1 * cosphi * coslam;
106	0	oblcon:
107	0	if (xy.y <= EPS10) {
108	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
109	0	return xy;
110	0	}
111	0	xy.y = sqrt(2. / xy.y);
112	0	xy.x = xy.y * cosphi * sin(lp.lam);
113	0	xy.y *= Q->mode == pj_laea_ns::EQUIT
114	0	? sinphi
115	0	: Q->cosb1 * sinphi - Q->sinb1 * cosphi * coslam;
116	0	break;
117	0	case pj_laea_ns::N_POLE:
118	0	coslam = -coslam;
119	0	PROJ_FALLTHROUGH;
120	0	case pj_laea_ns::S_POLE:
121	0	if (fabs(lp.phi + P->phi0) < EPS10) {
122	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
123	0	return xy;
124	0	}
125	0	xy.y = M_FORTPI - lp.phi * .5;
126	0	xy.y = 2. * (Q->mode == pj_laea_ns::S_POLE ? cos(xy.y) : sin(xy.y));
127	0	xy.x = xy.y * sin(lp.lam);
128	0	xy.y *= coslam;
129	0	break;
130	0	}
131	0	return xy;
132	0	}
133
134	0	static PJ_LP laea_e_inverse(PJ_XY xy, PJ P) { / Ellipsoidal, inverse */
135	0	PJ_LP lp = {0.0, 0.0};
136	0	struct pj_laea_data Q = static_cast<struct pj_laea_data >(P->opaque);
137	0	double cCe, sCe, q, rho, ab = 0.0;
138
139	0	switch (Q->mode) {
140	0	case pj_laea_ns::EQUIT:
141	0	case pj_laea_ns::OBLIQ: {
142	0	xy.x /= Q->dd;
143	0	xy.y *= Q->dd;
144	0	rho = hypot(xy.x, xy.y);
145	0	if (rho < EPS10) {
146	0	lp.lam = 0.;
147	0	lp.phi = P->phi0;
148	0	return lp;
149	0	}
150	0	const double asin_argument = .5 * rho / Q->rq;
151	0	if (asin_argument > 1) {
152	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
153	0	return lp;
154	0	}
155	0	sCe = 2. * asin(asin_argument);
156	0	cCe = cos(sCe);
157	0	sCe = sin(sCe);
158	0	xy.x *= sCe;
159	0	if (Q->mode == pj_laea_ns::OBLIQ) {
160	0	ab = cCe * Q->sinb1 + xy.y * sCe * Q->cosb1 / rho;
161	0	xy.y = rho * Q->cosb1 * cCe - xy.y * Q->sinb1 * sCe;
162	0	} else {
163	0	ab = xy.y * sCe / rho;
164	0	xy.y = rho * cCe;
165	0	}
166	0	break;
167	0	}
168	0	case pj_laea_ns::N_POLE:
169	0	xy.y = -xy.y;
170	0	PROJ_FALLTHROUGH;
171	0	case pj_laea_ns::S_POLE:
172	0	q = (xy.x * xy.x + xy.y * xy.y);
173	0	if (q == 0.0) {
174	0	lp.lam = 0.;
175	0	lp.phi = P->phi0;
176	0	return (lp);
177	0	}
178	0	ab = 1. - q / Q->qp;
179	0	if (Q->mode == pj_laea_ns::S_POLE)
180	0	ab = -ab;
181	0	break;
182	0	}
183	0	lp.lam = atan2(xy.x, xy.y);
184	0	lp.phi = pj_authalic_lat_inverse(asin(ab), Q->apa, P, Q->qp);
185	0	return lp;
186	0	}
187
188	0	static PJ_LP laea_s_inverse(PJ_XY xy, PJ P) { / Spheroidal, inverse */
189	0	PJ_LP lp = {0.0, 0.0};
190	0	struct pj_laea_data Q = static_cast<struct pj_laea_data >(P->opaque);
191	0	double cosz = 0.0, rh, sinz = 0.0;
192
193	0	rh = hypot(xy.x, xy.y);
194	0	if ((lp.phi = rh * .5) > 1.) {
195	0	proj_errno_set(P, PROJ_ERR_COORD_TRANSFM_OUTSIDE_PROJECTION_DOMAIN);
196	0	return lp;
197	0	}
198	0	lp.phi = 2. * asin(lp.phi);
199	0	if (Q->mode == pj_laea_ns::OBLIQ \|\| Q->mode == pj_laea_ns::EQUIT) {
200	0	sinz = sin(lp.phi);
201	0	cosz = cos(lp.phi);
202	0	}
203	0	switch (Q->mode) {
204	0	case pj_laea_ns::EQUIT:
205	0	lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
206	0	xy.x *= sinz;
207	0	xy.y = cosz * rh;
208	0	break;
209	0	case pj_laea_ns::OBLIQ:
210	0	lp.phi = fabs(rh) <= EPS10
211	0	? P->phi0
212	0	: asin(cosz * Q->sinb1 + xy.y * sinz * Q->cosb1 / rh);
213	0	xy.x = sinz Q->cosb1;
214	0	xy.y = (cosz - sin(lp.phi) * Q->sinb1) * rh;
215	0	break;
216	0	case pj_laea_ns::N_POLE:
217	0	xy.y = -xy.y;
218	0	lp.phi = M_HALFPI - lp.phi;
219	0	break;
220	0	case pj_laea_ns::S_POLE:
221	0	lp.phi -= M_HALFPI;
222	0	break;
223	0	}
224	0	lp.lam = (xy.y == 0. &&
225	0	(Q->mode == pj_laea_ns::EQUIT \|\| Q->mode == pj_laea_ns::OBLIQ))
226	0	? 0.
227	0	: atan2(xy.x, xy.y);
228	0	return (lp);
229	0	}
230
231	40.3k	static PJ pj_laea_destructor(PJ P, int errlev) {
232	40.3k	if (nullptr == P)
233	0	return nullptr;
234
235	40.3k	if (nullptr == P->opaque)
236	0	return pj_default_destructor(P, errlev);
237
238	40.3k	free(static_cast<struct pj_laea_data *>(P->opaque)->apa);
239
240	40.3k	return pj_default_destructor(P, errlev);
241	40.3k	}
242
243	40.3k	PJ *PJ_PROJECTION(laea) {
244	40.3k	double t;
245	40.3k	struct pj_laea_data Q = static_cast<struct pj_laea_data >(
246	40.3k	calloc(1, sizeof(struct pj_laea_data)));
247	40.3k	if (nullptr == Q)
248	0	return pj_default_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
249	40.3k	P->opaque = Q;
250	40.3k	P->destructor = pj_laea_destructor;
251
252	40.3k	t = fabs(P->phi0);
253	40.3k	if (t > M_HALFPI + EPS10) {
254	0	proj_log_error(P,
255	0	_("Invalid value for lat_0: \|lat_0\| should be <= 90°"));
256	0	return pj_laea_destructor(P, PROJ_ERR_INVALID_OP_ILLEGAL_ARG_VALUE);
257	0	}
258	40.3k	if (fabs(t - M_HALFPI) < EPS10)
259	11.3k	Q->mode = P->phi0 < 0. ? pj_laea_ns::S_POLE : pj_laea_ns::N_POLE;
260	29.0k	else if (fabs(t) < EPS10)
261	18.0k	Q->mode = pj_laea_ns::EQUIT;
262	10.9k	else
263	10.9k	Q->mode = pj_laea_ns::OBLIQ;
264	40.3k	if (P->es != 0.0) {
265	27.1k	double sinphi, cosphi;
266
267	27.1k	P->e = sqrt(P->es);
268	27.1k	Q->qp = pj_authalic_lat_q(1.0, P);
269	27.1k	Q->mmf = .5 / (1. - P->es);
270	27.1k	Q->apa = pj_authalic_lat_compute_coeffs(P->n);
271	27.1k	if (nullptr == Q->apa)
272	0	return pj_laea_destructor(P, PROJ_ERR_OTHER /ENOMEM/);
273	27.1k	switch (Q->mode) {
274	7.11k	case pj_laea_ns::N_POLE:
275	9.78k	case pj_laea_ns::S_POLE:
276	9.78k	Q->dd = 1.;
277	9.78k	break;
278	11.4k	case pj_laea_ns::EQUIT:
279	11.4k	Q->dd = 1. / (Q->rq = sqrt(.5 * Q->qp));
280	11.4k	Q->xmf = 1.;
281	11.4k	Q->ymf = .5 * Q->qp;
282	11.4k	break;
283	5.96k	case pj_laea_ns::OBLIQ:
284	5.96k	Q->rq = sqrt(.5 * Q->qp);
285	5.96k	sinphi = sin(P->phi0); cosphi = cos(P->phi0);
286	5.96k	const double b1 = pj_authalic_lat(P->phi0, sinphi, cosphi,
287	5.96k	Q->apa, P, Q->qp);
288	5.96k	Q->sinb1 = sin(b1);
289	5.96k	Q->cosb1 = cos(b1);
290	5.96k	Q->dd = cos(P->phi0) /
291	5.96k	(sqrt(1. - P->es * sinphi * sinphi) * Q->rq * Q->cosb1);
292	5.96k	Q->ymf = (Q->xmf = Q->rq) / Q->dd;
293	5.96k	Q->xmf *= Q->dd;
294	5.96k	break;
295	27.1k	}
296	27.1k	P->inv = laea_e_inverse;
297	27.1k	P->fwd = laea_e_forward;
298	27.1k	} else {
299	13.1k	if (Q->mode == pj_laea_ns::OBLIQ) {
300	5.01k	Q->sinb1 = sin(P->phi0);
301	5.01k	Q->cosb1 = cos(P->phi0);
302	5.01k	}
303	13.1k	P->inv = laea_s_inverse;
304	13.1k	P->fwd = laea_s_forward;
305	13.1k	}
306
307	40.3k	return P;
308	40.3k	}