Converted laea

2 files changed, 331 insertions, 225 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 53a3d49c..382116b2 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -362,7 +362,6 @@ int pj_geocent_selftest (void) {return 10000;}
 int pj_gs48_selftest (void) {return 10000;}
 int pj_gs50_selftest (void) {return 10000;}
 
-int pj_laea_selftest (void) {return 10000;}
 int pj_lagrng_selftest (void) {return 10000;}
 int pj_larr_selftest (void) {return 10000;}
 int pj_lask_selftest (void) {return 10000;}
diff --git a/src/PJ_laea.c b/src/PJ_laea.c
index 445e39c3..b58b07e2 100644
--- a/src/PJ_laea.c
+++ b/src/PJ_laea.c
@@ -1,233 +1,340 @@
-#define PROJ_PARMS__ \
-	double	sinb1; \
-	double	cosb1; \
-	double	xmf; \
-	double	ymf; \
-	double	mmf; \
-	double	qp; \
-	double	dd; \
-	double	rq; \
-	double	*apa; \
-	int		mode;
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(laea, "Lambert Azimuthal Equal Area") "\n\tAzi, Sph&Ell";
-#define sinph0	P->sinb1
-#define cosph0	P->cosb1
-#define EPS10	1.e-10
-#define NITER	20
-#define CONV	1.e-10
-#define N_POLE	0
-#define S_POLE	1
-#define EQUIT	2
-#define OBLIQ	3
-FORWARD(e_forward); /* ellipsoid */
-	double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
-
-	coslam = cos(lp.lam);
-	sinlam = sin(lp.lam);
-	sinphi = sin(lp.phi);
-	q = pj_qsfn(sinphi, P->e, P->one_es);
-	if (P->mode == OBLIQ || P->mode == EQUIT) {
-		sinb = q / P->qp;
-		cosb = sqrt(1. - sinb * sinb);
-	}
-	switch (P->mode) {
-	case OBLIQ:
-		b = 1. + P->sinb1 * sinb + P->cosb1 * cosb * coslam;
-		break;
-	case EQUIT:
-		b = 1. + cosb * coslam;
-		break;
-	case N_POLE:
-		b = HALFPI + lp.phi;
-		q = P->qp - q;
-		break;
-	case S_POLE:
-		b = lp.phi - HALFPI;
-		q = P->qp + q;
-		break;
-	}
-	if (fabs(b) < EPS10) F_ERROR;
-	switch (P->mode) {
-	case OBLIQ:
-		xy.y = P->ymf * ( b = sqrt(2. / b) )
-		   * (P->cosb1 * sinb - P->sinb1 * cosb * coslam);
-		goto eqcon;
-		break;
-	case EQUIT:
-		xy.y = (b = sqrt(2. / (1. + cosb * coslam))) * sinb * P->ymf; 
+
+struct pj_opaque {
+    double sinb1;
+    double cosb1;
+    double xmf;
+    double ymf;
+    double mmf;
+    double qp;
+    double dd;
+    double rq;
+    double *apa;
+    int mode;
+};
+
+#define EPS10   1.e-10
+#define NITER   20
+#define CONV    1.e-10
+#define N_POLE  0
+#define S_POLE  1
+#define EQUIT   2
+#define OBLIQ   3
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double coslam, sinlam, sinphi, q, sinb=0.0, cosb=0.0, b=0.0;
+
+    coslam = cos(lp.lam);
+    sinlam = sin(lp.lam);
+    sinphi = sin(lp.phi);
+    q = pj_qsfn(sinphi, P->e, P->one_es);
+
+    if (Q->mode == OBLIQ || Q->mode == EQUIT) {
+        sinb = q / Q->qp;
+        cosb = sqrt(1. - sinb * sinb);
+    }
+
+    switch (Q->mode) {
+    case OBLIQ:
+        b = 1. + Q->sinb1 * sinb + Q->cosb1 * cosb * coslam;
+        break;
+    case EQUIT:
+        b = 1. + cosb * coslam;
+        break;
+    case N_POLE:
+        b = HALFPI + lp.phi;
+        q = Q->qp - q;
+        break;
+    case S_POLE:
+        b = lp.phi - HALFPI;
+        q = Q->qp + q;
+        break;
+    }
+    if (fabs(b) < EPS10) F_ERROR;
+
+    switch (Q->mode) {
+    case OBLIQ:
+        b = sqrt(2. / b);
+        xy.y = Q->ymf * b * (Q->cosb1 * sinb - Q->sinb1 * cosb * coslam);
+        goto eqcon;
+        break;
+    case EQUIT:
+        b = sqrt(2. / (1. + cosb * coslam));
+        xy.y = b * sinb * Q->ymf;
 eqcon:
-		xy.x = P->xmf * b * cosb * sinlam;
-		break;
-	case N_POLE:
-	case S_POLE:
-		if (q >= 0.) {
-			xy.x = (b = sqrt(q)) * sinlam;
-			xy.y = coslam * (P->mode == S_POLE ? b : -b);
-		} else
-			xy.x = xy.y = 0.;
-		break;
-	}
-	return (xy);
+        xy.x = Q->xmf * b * cosb * sinlam;
+        break;
+    case N_POLE:
+    case S_POLE:
+        if (q >= 0.) {
+            b = sqrt(q);
+            xy.x = b * sinlam;
+            xy.y = coslam * (Q->mode == S_POLE ? b : -b);
+        } else
+            xy.x = xy.y = 0.;
+        break;
+    }
+    return xy;
 }
-FORWARD(s_forward); /* spheroid */
-	double  coslam, cosphi, sinphi;
-
-	sinphi = sin(lp.phi);
-	cosphi = cos(lp.phi);
-	coslam = cos(lp.lam);
-	switch (P->mode) {
-	case EQUIT:
-		xy.y = 1. + cosphi * coslam;
-		goto oblcon;
-	case OBLIQ:
-		xy.y = 1. + sinph0 * sinphi + cosph0 * cosphi * coslam;
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double  coslam, cosphi, sinphi;
+
+    sinphi = sin(lp.phi);
+    cosphi = cos(lp.phi);
+    coslam = cos(lp.lam);
+    switch (Q->mode) {
+    case EQUIT:
+        xy.y = 1. + cosphi * coslam;
+        goto oblcon;
+    case OBLIQ:
+        xy.y = 1. + Q->sinb1 * sinphi + Q->cosb1 * cosphi * coslam;
 oblcon:
-		if (xy.y <= EPS10) F_ERROR;
-		xy.x = (xy.y = sqrt(2. / xy.y)) * cosphi * sin(lp.lam);
-		xy.y *= P->mode == EQUIT ? sinphi :
-		   cosph0 * sinphi - sinph0 * cosphi * coslam;
-		break;
-	case N_POLE:
-		coslam = -coslam;
-	case S_POLE:
-		if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
-		xy.y = FORTPI - lp.phi * .5;
-		xy.y = 2. * (P->mode == S_POLE ? cos(xy.y) : sin(xy.y));
-		xy.x = xy.y * sin(lp.lam);
-		xy.y *= coslam;
-		break;
-	}
-	return (xy);
+        if (xy.y <= EPS10) F_ERROR;
+        xy.y = sqrt(2. / xy.y);
+        xy.x = xy.y * cosphi * sin(lp.lam);
+        xy.y *= Q->mode == EQUIT ? sinphi :
+           Q->cosb1 * sinphi - Q->sinb1 * cosphi * coslam;
+        break;
+    case N_POLE:
+        coslam = -coslam;
+    case S_POLE:
+        if (fabs(lp.phi + P->phi0) < EPS10) F_ERROR;
+        xy.y = FORTPI - lp.phi * .5;
+        xy.y = 2. * (Q->mode == S_POLE ? cos(xy.y) : sin(xy.y));
+        xy.x = xy.y * sin(lp.lam);
+        xy.y *= coslam;
+        break;
+    }
+    return xy;
+}
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double cCe, sCe, q, rho, ab=0.0;
+
+    switch (Q->mode) {
+    case EQUIT:
+    case 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;
+        }
+        sCe = 2. * asin(.5 * rho / Q->rq);
+        cCe = cos(sCe);
+        sCe = sin(sCe);
+        xy.x *= sCe;
+        if (Q->mode == 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 N_POLE:
+        xy.y = -xy.y;
+    case S_POLE:
+        q = (xy.x * xy.x + xy.y * xy.y);
+        if (!q) {
+            lp.lam = 0.;
+            lp.phi = P->phi0;
+            return (lp);
+        }
+        ab = 1. - q / Q->qp;
+        if (Q->mode == S_POLE)
+            ab = - ab;
+        break;
+    }
+    lp.lam = atan2(xy.x, xy.y);
+    lp.phi = pj_authlat(asin(ab), Q->apa);
+    return lp;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double  cosz=0.0, rh, sinz=0.0;
+
+    rh = hypot(xy.x, xy.y);
+    if ((lp.phi = rh * .5 ) > 1.) I_ERROR;
+    lp.phi = 2. * asin(lp.phi);
+    if (Q->mode == OBLIQ || Q->mode == EQUIT) {
+        sinz = sin(lp.phi);
+        cosz = cos(lp.phi);
+    }
+    switch (Q->mode) {
+    case EQUIT:
+        lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
+        xy.x *= sinz;
+        xy.y = cosz * rh;
+        break;
+    case 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 N_POLE:
+        xy.y = -xy.y;
+        lp.phi = HALFPI - lp.phi;
+        break;
+    case S_POLE:
+        lp.phi -= HALFPI;
+        break;
+    }
+    lp.lam = (xy.y == 0. && (Q->mode == EQUIT || Q->mode == OBLIQ)) ?
+        0. : atan2(xy.x, xy.y);
+    return (lp);
+}
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    pj_dealloc (P->opaque->apa);
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
 }
-INVERSE(e_inverse); /* ellipsoid */
-	double cCe, sCe, q, rho, ab=0.0;
-
-	switch (P->mode) {
-	case EQUIT:
-	case OBLIQ:
-		if ((rho = hypot(xy.x /= P->dd, xy.y *=  P->dd)) < EPS10) {
-			lp.lam = 0.;
-			lp.phi = P->phi0;
-			return (lp);
-		}
-		cCe = cos(sCe = 2. * asin(.5 * rho / P->rq));
-		xy.x *= (sCe = sin(sCe));
-		if (P->mode == OBLIQ) {
-			ab = cCe * P->sinb1 + xy.y * sCe * P->cosb1 / rho;
-			xy.y = rho * P->cosb1 * cCe - xy.y * P->sinb1 * sCe;
-		} else {
-			ab = xy.y * sCe / rho;
-			xy.y = rho * cCe;
-		}
-		break;
-	case N_POLE:
-		xy.y = -xy.y;
-	case S_POLE:
-		if (!(q = (xy.x * xy.x + xy.y * xy.y)) ) {
-			lp.lam = 0.;
-			lp.phi = P->phi0;
-			return (lp);
-		}
-		/*
-		q = P->qp - q;
-		*/
-		ab = 1. - q / P->qp;
-		if (P->mode == S_POLE)
-			ab = - ab;
-		break;
-	}
-	lp.lam = atan2(xy.x, xy.y);
-	lp.phi = pj_authlat(asin(ab), P->apa);
-	return (lp);
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-INVERSE(s_inverse); /* spheroid */
-	double  cosz=0.0, rh, sinz=0.0;
-
-	rh = hypot(xy.x, xy.y);
-	if ((lp.phi = rh * .5 ) > 1.) I_ERROR;
-	lp.phi = 2. * asin(lp.phi);
-	if (P->mode == OBLIQ || P->mode == EQUIT) {
-		sinz = sin(lp.phi);
-		cosz = cos(lp.phi);
-	}
-	switch (P->mode) {
-	case EQUIT:
-		lp.phi = fabs(rh) <= EPS10 ? 0. : asin(xy.y * sinz / rh);
-		xy.x *= sinz;
-		xy.y = cosz * rh;
-		break;
-	case OBLIQ:
-		lp.phi = fabs(rh) <= EPS10 ? P->phi0 :
-		   asin(cosz * sinph0 + xy.y * sinz * cosph0 / rh);
-		xy.x *= sinz * cosph0;
-		xy.y = (cosz - sin(lp.phi) * sinph0) * rh;
-		break;
-	case N_POLE:
-		xy.y = -xy.y;
-		lp.phi = HALFPI - lp.phi;
-		break;
-	case S_POLE:
-		lp.phi -= HALFPI;
-		break;
-	}
-	lp.lam = (xy.y == 0. && (P->mode == EQUIT || P->mode == OBLIQ)) ?
-		0. : atan2(xy.x, xy.y);
-	return (lp);
+
+
+PJ *PROJECTION(laea) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    P->pfree = freeup;
+    P->descr = des_laea;
+    double t;
+    t = fabs(P->phi0);
+    if (fabs(t - HALFPI) < EPS10)
+        Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
+    else if (fabs(t) < EPS10)
+        Q->mode = EQUIT;
+    else
+        Q->mode = OBLIQ;
+    if (P->es) {
+        double sinphi;
+
+        P->e = sqrt(P->es);
+        Q->qp = pj_qsfn(1., P->e, P->one_es);
+        Q->mmf = .5 / (1. - P->es);
+        Q->apa = pj_authset(P->es);
+        switch (Q->mode) {
+        case N_POLE:
+        case S_POLE:
+            Q->dd = 1.;
+            break;
+        case EQUIT:
+            Q->dd = 1. / (Q->rq = sqrt(.5 * Q->qp));
+            Q->xmf = 1.;
+            Q->ymf = .5 * Q->qp;
+            break;
+        case OBLIQ:
+            Q->rq = sqrt(.5 * Q->qp);
+            sinphi = sin(P->phi0);
+            Q->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / Q->qp;
+            Q->cosb1 = sqrt(1. - Q->sinb1 * Q->sinb1);
+            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 = e_inverse;
+        P->fwd = e_forward;
+    } else {
+        if (Q->mode == OBLIQ) {
+            Q->sinb1 = sin(P->phi0);
+            Q->cosb1 = cos(P->phi0);
+        }
+        P->inv = s_inverse;
+        P->fwd = s_forward;
+    }
+
+    return P;
 }
-FREEUP;
-    if (P) {
-		if (P->apa)
-			pj_dalloc(P->apa);
-		pj_dalloc(P);
-	}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_laea_selftest (void) {return 0;}
+#else
+
+int pj_laea_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=laea   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=laea   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222602.471450095181,  110589.82722441027},
+        { 222602.471450095181, -110589.827224408786},
+        {-222602.471450095181,  110589.82722441027},
+        {-222602.471450095181, -110589.827224408786},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223365.281370124663,  111716.668072915665},
+        { 223365.281370124663, -111716.668072915665},
+        {-223365.281370124663,  111716.668072915665},
+        {-223365.281370124663, -111716.668072915665},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.00179663056847900867,  0.000904369475966495845},
+        { 0.00179663056847900867, -0.000904369475966495845},
+        {-0.00179663056847900867,  0.000904369475966495845},
+        {-0.00179663056847900867, -0.000904369475966495845},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.00179049311002060264,  0.000895246554791735271},
+        { 0.00179049311002060264, -0.000895246554791735271},
+        {-0.00179049311002060264,  0.000895246554791735271},
+        {-0.00179049311002060264, -0.000895246554791735271},
+    };
+
+    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);
 }
-ENTRY1(laea,apa)
-	double t;
-
-	if (fabs((t = fabs(P->phi0)) - HALFPI) < EPS10)
-		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
-	else if (fabs(t) < EPS10)
-		P->mode = EQUIT;
-	else
-		P->mode = OBLIQ;
-	if (P->es) {
-		double sinphi;
-
-		P->e = sqrt(P->es);
-		P->qp = pj_qsfn(1., P->e, P->one_es);
-		P->mmf = .5 / (1. - P->es);
-		P->apa = pj_authset(P->es);
-		switch (P->mode) {
-		case N_POLE:
-		case S_POLE:
-			P->dd = 1.;
-			break;
-		case EQUIT:
-			P->dd = 1. / (P->rq = sqrt(.5 * P->qp));
-			P->xmf = 1.;
-			P->ymf = .5 * P->qp;
-			break;
-		case OBLIQ:
-			P->rq = sqrt(.5 * P->qp);
-			sinphi = sin(P->phi0);
-			P->sinb1 = pj_qsfn(sinphi, P->e, P->one_es) / P->qp;
-			P->cosb1 = sqrt(1. - P->sinb1 * P->sinb1);
-			P->dd = cos(P->phi0) / (sqrt(1. - P->es * sinphi * sinphi) *
-			   P->rq * P->cosb1);
-			P->ymf = (P->xmf = P->rq) / P->dd;
-			P->xmf *= P->dd;
-			break;
-		}
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
-		if (P->mode == OBLIQ) {
-			sinph0 = sin(P->phi0);
-			cosph0 = cos(P->phi0);
-		}
-		P->inv = s_inverse;
-		P->fwd = s_forward;
-	}
-ENDENTRY(P)
+
+
+#endif