Merge pull request #5 from kbevers/fix-projs-with-g-h-i-k

Converted projections starting with g h i and k

11 files changed, 3197 insertions, 2109 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 4db5e386..72091132 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -355,24 +355,13 @@ source files
 
 int pj_aeqd_selftest     (void)     {return 10000;}
 int pj_alsk_selftest     (void)     {return 10000;}
-int pj_eck6_selftest     (void)     {return 10000;}
 int pj_eqc_selftest (void) {return 10000;}
 int pj_eqdc_selftest (void) {return 10000;}
 int pj_etmerc_selftest (void) {return 10000;}
 int pj_geocent_selftest (void) {return 10000;}
-int pj_gn_sinu_selftest (void) {return 10000;}
-int pj_goode_selftest (void) {return 10000;}
 int pj_gs48_selftest (void) {return 10000;}
 int pj_gs50_selftest (void) {return 10000;}
-int pj_hammer_selftest (void) {return 10000;}
-int pj_hatano_selftest (void) {return 10000;}
-int pj_healpix_selftest (void) {return 10000;}
-int pj_rhealpix_selftest (void) {return 10000;}
-int pj_igh_selftest (void) {return 10000;}
-int pj_imw_p_selftest (void) {return 10000;}
-int pj_isea_selftest (void) {return 10000;}
-
-int pj_krovak_selftest (void) {return 10000;}
+
 int pj_labrd_selftest (void) {return 10000;}
 int pj_laea_selftest (void) {return 10000;}
 int pj_lagrng_selftest (void) {return 10000;}
@@ -392,7 +381,6 @@ int pj_lsat_selftest (void) {return 10000;}
 int pj_mbt_fps_selftest (void) {return 10000;}
 int pj_mbtfpp_selftest (void) {return 10000;}
 int pj_mbtfpq_selftest (void) {return 10000;}
-int pj_mbtfps_selftest (void) {return 10000;}
 int pj_merc_selftest (void) {return 10000;}
 int pj_mil_os_selftest (void) {return 10000;}
 int pj_mill_selftest (void) {return 10000;}
@@ -426,9 +414,7 @@ int pj_robin_selftest (void) {return 10000;}
 int pj_rouss_selftest (void) {return 10000;}
 int pj_rpoly_selftest (void) {return 10000;}
 int pj_sch_selftest (void) {return 10000;}
-int pj_sinu_selftest (void) {return 10000;}
 
-int pj_gstmerc_selftest (void) {return 10000;}
 int pj_tpers_selftest (void) {return 10000;}
 
 int pj_utm_selftest (void) {return 10000;}
diff --git a/src/PJ_gn_sinu.c b/src/PJ_gn_sinu.c
index bfd8bc2d..4b33b185 100644
--- a/src/PJ_gn_sinu.c
+++ b/src/PJ_gn_sinu.c
@@ -1,98 +1,373 @@
-#define PROJ_PARMS__ \
-	double	*en; \
-	double	m, n, C_x, C_y;
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(gn_sinu, "General Sinusoidal Series") "\n\tPCyl, Sph.\n\tm= n=";
 PROJ_HEAD(sinu, "Sinusoidal (Sanson-Flamsteed)") "\n\tPCyl, Sph&Ell";
 PROJ_HEAD(eck6, "Eckert VI") "\n\tPCyl, Sph.";
 PROJ_HEAD(mbtfps, "McBryde-Thomas Flat-Polar Sinusoidal") "\n\tPCyl, Sph.";
-#define EPS10	1e-10
+
+#define EPS10    1e-10
 #define MAX_ITER 8
 #define LOOP_TOL 1e-7
-/* Ellipsoidal Sinusoidal only */
-FORWARD(e_forward); /* ellipsoid */
-	double s, c;
 
-	xy.y = pj_mlfn(lp.phi, s = sin(lp.phi), c = cos(lp.phi), P->en);
-	xy.x = lp.lam * c / sqrt(1. - P->es * s * s);
-	return (xy);
+struct pj_opaque {
+    double *en;
+    double m, n, C_x, C_y;
+};
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    double s, c;
+
+    xy.y = pj_mlfn(lp.phi, s = sin(lp.phi), c = cos(lp.phi), P->opaque->en);
+    xy.x = lp.lam * c / sqrt(1. - P->es * s * s);
+    return xy;
 }
-INVERSE(e_inverse); /* ellipsoid */
-	double s;
-
-	if ((s = fabs(lp.phi = pj_inv_mlfn(P->ctx, xy.y, P->es, P->en))) < HALFPI) {
-		s = sin(lp.phi);
-		lp.lam = xy.x * sqrt(1. - P->es * s * s) / cos(lp.phi);
-	} else if ((s - EPS10) < HALFPI)
-		lp.lam = 0.;
-	else I_ERROR;
-	return (lp);
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    double s;
+
+    if ((s = fabs(lp.phi = pj_inv_mlfn(P->ctx, xy.y, P->es, P->opaque->en))) < HALFPI) {
+        s = sin(lp.phi);
+        lp.lam = xy.x * sqrt(1. - P->es * s * s) / cos(lp.phi);
+    } else if ((s - EPS10) < HALFPI) {
+        lp.lam = 0.;
+    } else {
+        I_ERROR;
+    }
+
+    return lp;
 }
-/* General spherical sinusoidals */
-FORWARD(s_forward); /* sphere */
-	if (!P->m)
-		lp.phi = P->n != 1. ? aasin(P->ctx,P->n * sin(lp.phi)): lp.phi;
-	else {
-		double k, V;
-		int i;
-
-		k = P->n * sin(lp.phi);
-		for (i = MAX_ITER; i ; --i) {
-			lp.phi -= V = (P->m * lp.phi + sin(lp.phi) - k) /
-				(P->m + cos(lp.phi));
-			if (fabs(V) < LOOP_TOL)
-				break;
-		}
-		if (!i)
-			F_ERROR
-	}
-	xy.x = P->C_x * lp.lam * (P->m + cos(lp.phi));
-	xy.y = P->C_y * lp.phi;
-	return (xy);
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    if (!Q->m)
+        lp.phi = Q->n != 1. ? aasin(P->ctx,Q->n * sin(lp.phi)): lp.phi;
+    else {
+        double k, V;
+        int i;
+
+        k = Q->n * sin(lp.phi);
+        for (i = MAX_ITER; i ; --i) {
+            lp.phi -= V = (Q->m * lp.phi + sin(lp.phi) - k) /
+                (Q->m + cos(lp.phi));
+            if (fabs(V) < LOOP_TOL)
+                break;
+        }
+        if (!i)
+            F_ERROR
+    }
+    xy.x = Q->C_x * lp.lam * (Q->m + cos(lp.phi));
+    xy.y = Q->C_y * lp.phi;
+
+    return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    xy.y /= Q->C_y;
+    lp.phi = Q->m ? aasin(P->ctx,(Q->m * xy.y + sin(xy.y)) / Q->n) :
+        ( Q->n != 1. ? aasin(P->ctx,sin(xy.y) / Q->n) : xy.y );
+    lp.lam = xy.x / (Q->C_x * (Q->m + cos(xy.y)));
+    return lp;
+}
+
+
+static void *freeup_new (PJ *P) {              /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    if (P->opaque->en)
+        pj_dalloc(P->opaque->en);
+
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
 }
-INVERSE(s_inverse); /* sphere */
-	xy.y /= P->C_y;
-	lp.phi = P->m ? aasin(P->ctx,(P->m * xy.y + sin(xy.y)) / P->n) :
-		( P->n != 1. ? aasin(P->ctx,sin(xy.y) / P->n) : xy.y );
-	lp.lam = xy.x / (P->C_x * (P->m + cos(xy.y)));
-	return (lp);
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
-	static void /* for spheres, only */
-setup(PJ *P) {
-	P->es = 0;
-	P->C_x = (P->C_y = sqrt((P->m + 1.) / P->n))/(P->m + 1.);
-	P->inv = s_inverse;
-	P->fwd = s_forward;
+
+
+/* for spheres, only */
+static void setup(PJ *P) {
+    struct pj_opaque *Q = P->opaque;
+    P->es = 0;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    Q->C_x = (Q->C_y = sqrt((Q->m + 1.) / Q->n))/(Q->m + 1.);
+}
+
+
+PJ *PROJECTION(sinu) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    if (!(Q->en = pj_enfn(P->es)))
+        E_ERROR_0;
+
+    if (P->es) {
+        P->inv = e_inverse;
+        P->fwd = e_forward;
+    } else {
+        Q->n = 1.;
+        Q->m = 0.;
+        setup(P);
+    }
+    return P;
+}
+
+
+PJ *PROJECTION(eck6) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->m = 1.;
+    Q->n = 2.570796326794896619231321691;
+    setup(P);
+
+    return P;
+}
+
+
+PJ *PROJECTION(mbtfps) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->m = 0.5;
+    Q->n = 1.785398163397448309615660845;
+    setup(P);
+
+    return P;
+}
+
+
+PJ *PROJECTION(gn_sinu) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    if (pj_param(P->ctx, P->params, "tn").i && pj_param(P->ctx, P->params, "tm").i) {
+        Q->n = pj_param(P->ctx, P->params, "dn").f;
+        Q->m = pj_param(P->ctx, P->params, "dm").f;
+    } else
+        E_ERROR(-99)
+
+    setup(P);
+
+    return P;
 }
-ENTRY1(sinu, en)
-	if (!(P->en = pj_enfn(P->es)))
-		E_ERROR_0;
-	if (P->es) {
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
-		P->n = 1.;
-		P->m = 0.;
-		setup(P);
-	}
-ENDENTRY(P)
-ENTRY1(eck6, en)
-	P->m = 1.;
-	P->n = 2.570796326794896619231321691;
-	setup(P);
-ENDENTRY(P)
-ENTRY1(mbtfps, en)
-	P->m = 0.5;
-	P->n = 1.785398163397448309615660845;
-	setup(P);
-ENDENTRY(P)
-ENTRY1(gn_sinu, en)
-	if (pj_param(P->ctx, P->params, "tn").i && pj_param(P->ctx, P->params, "tm").i) {
-		P->n = pj_param(P->ctx, P->params, "dn").f;
-		P->m = pj_param(P->ctx, P->params, "dm").f;
-	} else
-		E_ERROR(-99)
-	setup(P);
-ENDENTRY(P)
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_sinu_selftest (void) {return 0;}
+#else
+
+int pj_sinu_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=sinu   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=sinu   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222605.29953946592,  110574.38855415257},
+        { 222605.29953946592, -110574.38855415257},
+        {-222605.29953946592,  110574.38855415257},
+        {-222605.29953946592, -110574.38855415257},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223368.11902663155,  111701.07212763709},
+        { 223368.11902663155, -111701.07212763709},
+        {-223368.11902663155,  111701.07212763709},
+        {-223368.11902663155, -111701.07212763709},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017966305684613522,  0.00090436947707945409},
+        { 0.0017966305684613522, -0.00090436947707945409},
+        {-0.0017966305684613522,  0.00090436947707945409},
+        {-0.0017966305684613522, -0.00090436947707945409},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931100023887,  0.00089524655489191132},
+        { 0.0017904931100023887, -0.00089524655489191132},
+        {-0.0017904931100023887,  0.00089524655489191132},
+        {-0.0017904931100023887, -0.00089524655489191132},
+    };
+
+    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);
+}
+
+
+#endif
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_eck6_selftest (void) {return 0;}
+#else
+
+int pj_eck6_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=eck6   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 197021.60562899226,  126640.42073317352},
+        { 197021.60562899226, -126640.42073317352},
+        {-197021.60562899226,  126640.42073317352},
+        {-197021.60562899226, -126640.42073317352},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.002029978749734037,  0.00078963032910382171},
+        { 0.002029978749734037, -0.00078963032910382171},
+        {-0.002029978749734037,  0.00078963032910382171},
+        {-0.002029978749734037, -0.00078963032910382171},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_mbtfps_selftest (void) {return 0;}
+#else
+
+int pj_mbtfps_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=mbtfps   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 204740.11747857218,  121864.72971934026},
+        { 204740.11747857218, -121864.72971934026},
+        {-204740.11747857218,  121864.72971934026},
+        {-204740.11747857218, -121864.72971934026},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0019534152166442065,  0.00082057965689633387},
+        { 0.0019534152166442065, -0.00082057965689633387},
+        {-0.0019534152166442065,  0.00082057965689633387},
+        {-0.0019534152166442065, -0.00082057965689633387},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_gn_sinu_selftest (void) {return 0;}
+#else
+
+int pj_gn_sinu_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=gn_sinu   +a=6400000    +lat_1=0.5 +lat_2=2 +m=1 +n=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223385.13250469571,  111698.23644718733},
+        { 223385.13250469571, -111698.23644718733},
+        {-223385.13250469571,  111698.23644718733},
+        {-223385.13250469571, -111698.23644718733},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931098931057,  0.00089524655491012516},
+        { 0.0017904931098931057, -0.00089524655491012516},
+        {-0.0017904931098931057,  0.00089524655491012516},
+        {-0.0017904931098931057, -0.00089524655491012516},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_goode.c b/src/PJ_goode.c
index 387557e6..d577f901 100644
--- a/src/PJ_goode.c
+++ b/src/PJ_goode.c
@@ -1,49 +1,128 @@
-#define PROJ_PARMS__ \
-	struct PJconsts	*sinu; \
-	struct PJconsts	*moll;
-#define PJ_LIB__
-#include	<projects.h>
-PROJ_HEAD(goode, "Goode Homolosine") "\n\tPCyl, Sph.";
-	C_NAMESPACE PJ
-*pj_sinu(PJ *), *pj_moll(PJ *);
-#define Y_COR		0.05280
-#define PHI_LIM	.71093078197902358062
-FORWARD(s_forward); /* spheroid */
-	if (fabs(lp.phi) <= PHI_LIM)
-		xy = P->sinu->fwd(lp, P->sinu);
-	else {
-		xy = P->moll->fwd(lp, P->moll);
-		xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;
-	}
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	if (fabs(xy.y) <= PHI_LIM)
-		lp = P->sinu->inv(xy, P->sinu);
-	else {
-		xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;
-		lp = P->moll->inv(xy, P->moll);
-	}
-	return (lp);
-}
-FREEUP;
-	if (P) {
-		if (P->sinu)
-			(*(P->sinu->pfree))(P->sinu);
-		if (P->moll)
-			(*(P->moll->pfree))(P->moll);
-		pj_dalloc(P);
-	}
-}
-ENTRY2(goode, sinu, moll)
-	P->es = 0.;
-	if (!(P->sinu = pj_sinu(0)) || !(P->moll = pj_moll(0)))
-		E_ERROR_0;
-	P->sinu->es = 0.;
-        P->sinu->ctx = P->ctx;
-        P->moll->ctx = P->ctx;
-	if (!(P->sinu = pj_sinu(P->sinu)) || !(P->moll = pj_moll(P->moll)))
-		E_ERROR_0;
-	P->fwd = s_forward;
-	P->inv = s_inverse;
-ENDENTRY(P)
+#define PJ_LIB__

+#include <projects.h>

+

+PROJ_HEAD(goode, "Goode Homolosine") "\n\tPCyl, Sph.";

+

+#define Y_COR   0.05280

+#define PHI_LIM 0.71093078197902358062

+

+struct pj_opaque {

+    struct PJconsts *sinu; \

+    struct PJconsts *moll;

+};

+

+

+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */

+    XY xy = {0.0,0.0};

+    struct pj_opaque *Q = P->opaque;

+

+    if (fabs(lp.phi) <= PHI_LIM)

+        xy = Q->sinu->fwd(lp, Q->sinu);

+    else {

+        xy = Q->moll->fwd(lp, Q->moll);

+        xy.y -= lp.phi >= 0.0 ? Y_COR : -Y_COR;

+    }

+    return xy;

+}

+

+

+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */

+    LP lp = {0.0,0.0};

+    struct pj_opaque *Q = P->opaque;

+

+    if (fabs(xy.y) <= PHI_LIM)

+        lp = Q->sinu->inv(xy, Q->sinu);

+    else {

+        xy.y += xy.y >= 0.0 ? Y_COR : -Y_COR;

+        lp = Q->moll->inv(xy, Q->moll);

+    }

+    return lp;

+}

+

+

+static void *freeup_new (PJ *P) {              /* Destructor */

+    if (0==P)

+        return 0;

+    if (0==P->opaque)

+        return pj_dealloc(P);

+    if (P->opaque->sinu)

+        pj_dealloc(P->opaque->sinu);

+    if (P->opaque->moll)

+        pj_dealloc(P->opaque->moll);

+    pj_dealloc (P->opaque);

+    return pj_dealloc(P);

+

+}

+

+

+static void freeup (PJ *P) {

+    freeup_new (P);

+    return;

+}

+

+

+PJ *PROJECTION(goode) {

+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));

+    if (0==Q)

+        return freeup_new (P);

+    P->opaque = Q;

+

+    P->es = 0.;

+    if (!(Q->sinu = pj_sinu(0)) || !(Q->moll = pj_moll(0)))

+        E_ERROR_0;

+    Q->sinu->es = 0.;

+        Q->sinu->ctx = P->ctx;

+        Q->moll->ctx = P->ctx;

+    if (!(Q->sinu = pj_sinu(Q->sinu)) || !(Q->moll = pj_moll(Q->moll)))

+        E_ERROR_0;

+

+    P->fwd = s_forward;

+    P->inv = s_inverse;

+

+    return P;

+}

+

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_goode_selftest (void) {return 0;}

+#else

+

+int pj_goode_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

+    char s_args[] = {"+proj=goode   +a=6400000    +lat_1=0.5 +lat_2=2"};

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+

+    XY s_fwd_expect[] = {

+        { 223368.11902663155,  111701.07212763709
},

+        { 223368.11902663155, -111701.07212763709
},

+        {-223368.11902663155,  111701.07212763709
},

+        {-223368.11902663155, -111701.07212763709
},

+    };

+

+    XY inv_in[] = {

+        { 200, 100},

+        { 200,-100},

+        {-200, 100},

+        {-200,-100}

+    };

+

+    LP s_inv_expect[] = {

+        { 0.0017904931100023887,  0.00089524655489191132
},

+        { 0.0017904931100023887, -0.00089524655489191132
},

+        {-0.0017904931100023887,  0.00089524655489191132
},

+        {-0.0017904931100023887, -0.00089524655489191132
},

+    };

+

+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);

+}

+

+

+#endif

diff --git a/src/PJ_gstmerc.c b/src/PJ_gstmerc.c
index bffe0b26..2720c945 100644
--- a/src/PJ_gstmerc.c
+++ b/src/PJ_gstmerc.c
@@ -1,48 +1,131 @@
-#define PROJ_PARMS__ \
-	double lamc;\
-	double phic;\
-	double c;\
-	double n1;\
-	double n2;\
-    double XS;\
-    double YS;
-
 #define PJ_LIB__
-# include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(gstmerc, "Gauss-Schreiber Transverse Mercator (aka Gauss-Laborde Reunion)")
-	"\n\tCyl, Sph&Ell\n\tlat_0= lon_0= k_0=";
-FORWARD(s_forward); /* spheroid */
-	double L, Ls, sinLs1, Ls1;
-	L= P->n1*lp.lam;
-    Ls= P->c+P->n1*log(pj_tsfn(-1.0*lp.phi,-1.0*sin(lp.phi),P->e));
-    sinLs1= sin(L)/cosh(Ls);
-    Ls1= log(pj_tsfn(-1.0*asin(sinLs1),0.0,0.0));
-    xy.x= (P->XS + P->n2*Ls1)*P->ra;
-    xy.y= (P->YS + P->n2*atan(sinh(Ls)/cos(L)))*P->ra;
-    /*fprintf(stderr,"fwd:\nL      =%16.13f\nLs     =%16.13f\nLs1    =%16.13f\nLP(%16.13f,%16.13f)=XY(%16.4f,%16.4f)\n",L,Ls,Ls1,lp.lam+P->lam0,lp.phi,(xy.x*P->a + P->x0)*P->to_meter,(xy.y*P->a + P->y0)*P->to_meter);*/
-	return (xy);
+    "\n\tCyl, Sph&Ell\n\tlat_0= lon_0= k_0=";
+
+struct pj_opaque {
+    double lamc;
+    double phic;
+    double c;
+    double n1;
+    double n2;
+    double XS;
+    double YS;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double L, Ls, sinLs1, Ls1;
+
+    L = Q->n1*lp.lam;
+    Ls = Q->c + Q->n1 * log(pj_tsfn(-1.0 * lp.phi, -1.0 * sin(lp.phi), P->e));
+    sinLs1 = sin(L) / cosh(Ls);
+    Ls1 = log(pj_tsfn(-1.0 * asin(sinLs1), 0.0, 0.0));
+    xy.x = (Q->XS + Q->n2*Ls1) * P->ra;
+    xy.y = (Q->YS + Q->n2*atan(sinh(Ls) / cos(L))) * P->ra;
+
+    return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double L, LC, sinC;
+
+    L = atan(sinh((xy.x * P->a - Q->XS) / Q->n2) / cos((xy.y * P->a - Q->YS) / Q->n2));
+    sinC = sin((xy.y * P->a - Q->YS) / Q->n2) / cosh((xy.x * P->a - Q->XS) / Q->n2);
+    LC = log(pj_tsfn(-1.0 * asin(sinC), 0.0, 0.0));
+    lp.lam = L / Q->n1;
+    lp.phi = -1.0 * pj_phi2(P->ctx, exp((LC - Q->c) / Q->n1), P->e);
+
+    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);
+    return pj_dealloc(P);
 }
-INVERSE(s_inverse); /* spheroid */
-	double L, LC, sinC;
-	L= atan(sinh((xy.x*P->a - P->XS)/P->n2)/cos((xy.y*P->a - P->YS)/P->n2));
-    sinC= sin((xy.y*P->a - P->YS)/P->n2)/cosh((xy.x*P->a - P->XS)/P->n2);
-    LC= log(pj_tsfn(-1.0*asin(sinC),0.0,0.0));
-    lp.lam= L/P->n1;
-    lp.phi= -1.0*pj_phi2(P->ctx, exp((LC-P->c)/P->n1),P->e);
-    /*fprintf(stderr,"inv:\nL      =%16.13f\nsinC   =%16.13f\nLC     =%16.13f\nXY(%16.4f,%16.4f)=LP(%16.13f,%16.13f)\n",L,sinC,LC,((xy.x/P->ra)+P->x0)/P->to_meter,((xy.y/P->ra)+P->y0)/P->to_meter,lp.lam+P->lam0,lp.phi);*/
-	return (lp);
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(gstmerc)
-    P->lamc= P->lam0;
-    P->n1= sqrt(1.0+P->es*pow(cos(P->phi0),4.0)/(1.0-P->es));
-    P->phic= asin(sin(P->phi0)/P->n1);
-    P->c=       log(pj_tsfn(-1.0*P->phic,0.0,0.0))
-         -P->n1*log(pj_tsfn(-1.0*P->phi0,-1.0*sin(P->phi0),P->e));
-    P->n2= P->k0*P->a*sqrt(1.0-P->es)/(1.0-P->es*sin(P->phi0)*sin(P->phi0));
-    P->XS= 0;/* -P->x0 */
-    P->YS= -1.0*P->n2*P->phic;/* -P->y0 */
-	P->inv= s_inverse;
-	P->fwd= s_forward;
-    /*fprintf(stderr,"a  (m) =%16.4f\ne      =%16.13f\nl0(rad)=%16.13f\np0(rad)=%16.13f\nk0     =%16.4f\nX0  (m)=%16.4f\nY0  (m)=%16.4f\n\nlC(rad)=%16.13f\npC(rad)=%16.13f\nc      =%16.13f\nn1     =%16.13f\nn2 (m) =%16.4f\nXS (m) =%16.4f\nYS (m) =%16.4f\n", P->a, P->e, P->lam0, P->phi0, P->k0, P->x0, P->y0, P->lamc, P->phic, P->c, P->n1, P->n2, P->XS +P->x0, P->YS + P->y0);*/
-ENDENTRY(P)
+
+
+PJ *PROJECTION(gstmerc) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->lamc = P->lam0;
+    Q->n1 = sqrt(1.0 + P->es * pow(cos(P->phi0), 4.0) / (1.0 - P->es));
+    Q->phic = asin(sin(P->phi0) / Q->n1);
+    Q->c = log(pj_tsfn(-1.0 * Q->phic, 0.0, 0.0))
+         - Q->n1 * log(pj_tsfn(-1.0 * P->phi0, -1.0 * sin(P->phi0), P->e));
+    Q->n2 = P->k0 * P->a * sqrt(1.0 - P->es) / (1.0 - P->es * sin(P->phi0) * sin(P->phi0));
+    Q->XS = 0;
+    Q->YS = -1.0 * Q->n2 * Q->phic;
+
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_gstmerc_selftest (void) {return 0;}
+#else
+
+int pj_gstmerc_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=gstmerc   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+
+    XY s_fwd_expect[] = {
+        { 223413.46640632182,  111769.14504058557},
+        { 223413.46640632182, -111769.14504058668},
+        {-223413.46640632302,  111769.14504058557},
+        {-223413.46640632302, -111769.14504058668},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931097109673,  0.0008952465544509083},
+        { 0.0017904931097109673, -0.0008952465544509083},
+        {-0.0017904931097109673,  0.0008952465544509083},
+        {-0.0017904931097109673, -0.0008952465544509083},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
+
+#endif
diff --git a/src/PJ_hammer.c b/src/PJ_hammer.c
index 31e7a127..3118640b 100644
--- a/src/PJ_hammer.c
+++ b/src/PJ_hammer.c
@@ -1,43 +1,129 @@
-#define PROJ_PARMS__ \
-	double w; \
-	double m, rm;
 #define PJ_LIB__
-#define EPS 1.0e-10
-# include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(hammer, "Hammer & Eckert-Greifendorff")
-	"\n\tMisc Sph, \n\tW= M=";
-FORWARD(s_forward); /* spheroid */
-	double cosphi, d;
-
-	d = sqrt(2./(1. + (cosphi = cos(lp.phi)) * cos(lp.lam *= P->w)));
-	xy.x = P->m * d * cosphi * sin(lp.lam);
-	xy.y = P->rm * d * sin(lp.phi);
-	return (xy);
+    "\n\tMisc Sph, \n\tW= M=";
+
+#define EPS 1.0e-10
+
+struct pj_opaque {
+    double w; \
+    double m, rm;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double cosphi, d;
+
+    d = sqrt(2./(1. + (cosphi = cos(lp.phi)) * cos(lp.lam *= Q->w)));
+    xy.x = Q->m * d * cosphi * sin(lp.lam);
+    xy.y = Q->rm * d * sin(lp.phi);
+    return xy;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double z;
+
+    z = sqrt(1. - 0.25*Q->w*Q->w*xy.x*xy.x - 0.25*xy.y*xy.y);
+    if (fabs(2.*z*z-1.) < EPS) {
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_errno = -14;
+    } else {
+        lp.lam = aatan2(Q->w * xy.x * z,2. * z * z - 1)/Q->w;
+        lp.phi = aasin(P->ctx,z * 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);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-INVERSE(s_inverse); /* spheroid */
-	double z;
-	z = sqrt(1. - 0.25*P->w*P->w*xy.x*xy.x - 0.25*xy.y*xy.y);
-	if (fabs(2.*z*z-1.) < EPS) {
-		lp.lam = HUGE_VAL;
-		lp.phi = HUGE_VAL;
-		pj_errno = -14;
-	} else {
-		lp.lam = aatan2(P->w * xy.x * z,2. * z * z - 1)/P->w;
-		lp.phi = aasin(P->ctx,z * xy.y);
-	}
-	return (lp);
+
+
+PJ *PROJECTION(hammer) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    if (pj_param(P->ctx, P->params, "tW").i) {
+        if ((Q->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) E_ERROR(-27);
+    } else
+        Q->w = .5;
+    if (pj_param(P->ctx, P->params, "tM").i) {
+        if ((Q->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) E_ERROR(-27);
+    } else
+        Q->m = 1.;
+
+    Q->rm = 1. / Q->m;
+    Q->m /= Q->w;
+
+    P->es = 0.;
+    P->fwd = s_forward;
+    P->inv = s_inverse;
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_hammer_selftest (void) {return 0;}
+#else
+
+int pj_hammer_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=hammer   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223373.78870324057,  111703.90739776699},
+        { 223373.78870324057, -111703.90739776699},
+        {-223373.78870324057,  111703.90739776699},
+        {-223373.78870324057, -111703.90739776699},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.001790493109965961,  0.00089524655487369749},
+        { 0.001790493109965961, -0.00089524655487369749},
+        {-0.001790493109965961,  0.00089524655487369749},
+        {-0.001790493109965961, -0.00089524655487369749},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(hammer)
-	if (pj_param(P->ctx, P->params, "tW").i) {
-		if ((P->w = fabs(pj_param(P->ctx, P->params, "dW").f)) <= 0.) E_ERROR(-27);
-	} else
-		P->w = .5;
-	if (pj_param(P->ctx, P->params, "tM").i) {
-		if ((P->m = fabs(pj_param(P->ctx, P->params, "dM").f)) <= 0.) E_ERROR(-27);
-	} else
-		P->m = 1.;
-	P->rm = 1. / P->m;
-	P->m /= P->w;
-	P->es = 0.; P->fwd = s_forward; P->inv = s_inverse;
-ENDENTRY(P)
+
+
+#endif
diff --git a/src/PJ_hatano.c b/src/PJ_hatano.c
index 7516ba6e..ca97849b 100644
--- a/src/PJ_hatano.c
+++ b/src/PJ_hatano.c
@@ -1,51 +1,134 @@
-#define PJ_LIB__
-#include	<projects.h>
-PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph.";
-#define NITER	20
-#define EPS	1e-7
-#define ONETOL 1.000001
-#define CN	2.67595
-#define CS	2.43763
-#define RCN	0.37369906014686373063
-#define RCS	0.41023453108141924738
-#define FYCN	1.75859
-#define FYCS	1.93052
-#define RYCN	0.56863737426006061674
-#define RYCS	0.51799515156538134803
-#define FXC	0.85
-#define RXC	1.17647058823529411764
-FORWARD(s_forward); /* spheroid */
-	double th1, c;
-	int i;
-	(void) P;
-
-	c = sin(lp.phi) * (lp.phi < 0. ? CS : CN);
-	for (i = NITER; i; --i) {
-		lp.phi -= th1 = (lp.phi + sin(lp.phi) - c) / (1. + cos(lp.phi));
-		if (fabs(th1) < EPS) break;
-	}
-	xy.x = FXC * lp.lam * cos(lp.phi *= .5);
-	xy.y = sin(lp.phi) * (lp.phi < 0. ? FYCS : FYCN);
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	double th;
-
-	th = xy.y * ( xy.y < 0. ? RYCS : RYCN);
-	if (fabs(th) > 1.)
-		if (fabs(th) > ONETOL)	I_ERROR
-		else			th = th > 0. ? HALFPI : - HALFPI;
-	else
-		th = asin(th);
-	lp.lam = RXC * xy.x / cos(th);
-	th += th;
-	lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN);
-	if (fabs(lp.phi) > 1.)
-		if (fabs(lp.phi) > ONETOL)	I_ERROR
-		else			lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
-	else
-		lp.phi = asin(lp.phi);
-	return (lp);
-}
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(hatano) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
+#define PJ_LIB__

+#include <projects.h>

+

+PROJ_HEAD(hatano, "Hatano Asymmetrical Equal Area") "\n\tPCyl, Sph.";

+

+#define NITER   20

+#define EPS 1e-7

+#define ONETOL 1.000001

+#define CN  2.67595

+#define CS  2.43763

+#define RCN 0.37369906014686373063

+#define RCS 0.41023453108141924738

+#define FYCN    1.75859

+#define FYCS    1.93052

+#define RYCN    0.56863737426006061674

+#define RYCS    0.51799515156538134803

+#define FXC 0.85

+#define RXC 1.17647058823529411764

+

+

+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */

+    XY xy = {0.0,0.0};

+    double th1, c;

+    int i;

+    (void) P;

+

+    c = sin(lp.phi) * (lp.phi < 0. ? CS : CN);

+    for (i = NITER; i; --i) {

+        lp.phi -= th1 = (lp.phi + sin(lp.phi) - c) / (1. + cos(lp.phi));

+        if (fabs(th1) < EPS) break;

+    }

+    xy.x = FXC * lp.lam * cos(lp.phi *= .5);

+    xy.y = sin(lp.phi) * (lp.phi < 0. ? FYCS : FYCN);

+

+    return xy;

+}

+

+

+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */

+    LP lp = {0.0,0.0};

+    double th;

+

+    th = xy.y * ( xy.y < 0. ? RYCS : RYCN);

+    if (fabs(th) > 1.) {

+        if (fabs(th) > ONETOL) {

+            I_ERROR;

+        } else {

+            th = th > 0. ? HALFPI : - HALFPI;

+        }

+    } else {

+        th = asin(th);

+    }

+

+    lp.lam = RXC * xy.x / cos(th);

+    th += th;

+    lp.phi = (th + sin(th)) * (xy.y < 0. ? RCS : RCN);

+    if (fabs(lp.phi) > 1.) {

+        if (fabs(lp.phi) > ONETOL) {

+           I_ERROR;

+        } else {

+            lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;

+        }

+    } else {

+        lp.phi = asin(lp.phi);

+    }

+

+    return (lp);

+}

+

+

+static void *freeup_new (PJ *P) {                       /* Destructor */

+    if (0==P)

+        return 0;

+

+    return pj_dealloc(P);

+}

+

+static void freeup (PJ *P) {

+    freeup_new (P);

+    return;

+}

+

+

+PJ *PROJECTION(hatano) {

+    P->es = 0.;

+    P->inv = s_inverse;

+    P->fwd = s_forward;

+

+    return P;

+}

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_hatano_selftest (void) {return 0;}

+#else

+

+int pj_hatano_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

+    char s_args[] = {"+proj=hatano   +a=6400000    +lat_1=0.5 +lat_2=2"};

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+

+    XY s_fwd_expect[] = {

+        { 189878.87894652804,  131409.8024406255
},

+        { 189881.08195244463, -131409.14227607418
},

+        {-189878.87894652804,  131409.8024406255
},

+        {-189881.08195244463, -131409.14227607418
},

+    };

+

+    XY inv_in[] = {

+        { 200, 100},

+        { 200,-100},

+        {-200, 100},

+        {-200,-100}

+    };

+

+    LP s_inv_expect[] = {

+        { 0.0021064624821817597,  0.00076095689425791926
},

+        { 0.0021064624821676096, -0.00076095777439265377
},

+        {-0.0021064624821817597,  0.00076095689425791926
},

+        {-0.0021064624821676096, -0.00076095777439265377
},

+    };

+

+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);

+}

+

+

+#endif

diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index 12d57ff2..ae5aa5be 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -2,10 +2,10 @@
  * Project: PROJ.4
  * Purpose: Implementation of the HEALPix and rHEALPix projections.
  *          For background see <http://code.scenzgrid.org/index.php/p/scenzgrid-py/source/tree/master/docs/rhealpix_dggs.pdf>.
- * Authors: Alex Raichev (raichev@cs.auckland.ac.nz) 
+ * Authors: Alex Raichev (raichev@cs.auckland.ac.nz)
  *          Michael Speth (spethm@landcareresearch.co.nz)
- * Notes:   Raichev implemented these projections in Python and 
- *          Speth translated them into C here.   
+ * Notes:   Raichev implemented these projections in Python and
+ *          Speth translated them into C here.
  ******************************************************************************
  * Copyright (c) 2001, Thomas Flemming, tf@ttqv.com
  *
@@ -28,38 +28,46 @@
  * CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER DEALINGS IN THE
  * SOFTWARE.
  *****************************************************************************/
-# define PROJ_PARMS__ \
-    int north_square; \
-    int south_square; \
-    double qp; \
-    double *apa;
-# define PJ_LIB__    
-# include	<projects.h>
+# define PJ_LIB__
+# include <projects.h>
+
 PROJ_HEAD(healpix, "HEALPix") "\n\tSph., Ellps.";
 PROJ_HEAD(rhealpix, "rHEALPix") "\n\tSph., Ellps.\n\tnorth_square= south_square=";
-# include	<stdio.h>
+
+# include <stdio.h>
 /* Matrix for counterclockwise rotation by pi/2: */
-# define R1 {{ 0,-1},{ 1, 0}}	
+# define R1 {{ 0,-1},{ 1, 0}}
 /* Matrix for counterclockwise rotation by pi: */
-# define R2 {{-1, 0},{ 0,-1}}	
+# define R2 {{-1, 0},{ 0,-1}}
 /* Matrix for counterclockwise rotation by 3*pi/2:  */
-# define R3 {{ 0, 1},{-1, 0}}	
+# define R3 {{ 0, 1},{-1, 0}}
 /* Identity matrix */
 # define IDENT {{1, 0},{0, 1}}
 /* IDENT, R1, R2, R3, R1 inverse, R2 inverse, R3 inverse:*/
 # define ROT {IDENT, R1, R2, R3, R3, R2, R1}
 /* Fuzz to handle rounding errors: */
 # define EPS 1e-15
+
+struct pj_opaque {
+    int north_square; \
+    int south_square; \
+    double qp; \
+    double *apa;
+};
+
 typedef struct {
     int cn; /* An integer 0--3 indicating the position of the polar cap. */
     double x, y;  /* Coordinates of the pole point (point of most extreme latitude on the polar caps). */
     enum Region {north, south, equatorial} region;
 } CapMap;
+
 typedef struct {
     double x, y;
 } Point;
+
 double rot[7][2][2] = ROT;
 
+
 /**
  * Returns the sign of the double.
  * @param v the parameter whose sign is returned.
@@ -68,31 +76,35 @@ double rot[7][2][2] = ROT;
 double pj_sign (double v) {
     return v > 0 ? 1 : (v < 0 ? -1 : 0);
 }
+
+
 /**
  * Return the index of the matrix in ROT.
  * @param index ranges from -3 to 3.
  */
 static int get_rotate_index(int index) {
     switch(index) {
-	case 0:
-	    return 0;
-	case 1:
-	    return 1;
-	case 2: 
-	    return 2;
-	case 3: 
-	    return 3;
-	case -1:
-	    return 4;
-	case -2:
-	    return 5;
-	case -3:
-	    return 6;
+    case 0:
+        return 0;
+    case 1:
+        return 1;
+    case 2:
+        return 2;
+    case 3:
+        return 3;
+    case -1:
+        return 4;
+    case -2:
+        return 5;
+    case -3:
+        return 6;
     }
     return 0;
 }
+
+
 /**
- * Return 1 if point (testx, testy) lies in the interior of the polygon 
+ * Return 1 if point (testx, testy) lies in the interior of the polygon
  * determined by the vertices in vert, and return 0 otherwise.
  * See http://paulbourke.net/geometry/polygonmesh/ for more details.
  * @param nvert the number of vertices in the polygon.
@@ -103,6 +115,7 @@ static int pnpoly(int nvert, double vert[][2], double testx, double testy) {
     int counter = 0;
     double xinters;
     Point p1, p2;
+
     /* Check for boundrary cases */
     for (i = 0; i < nvert; i++) {
         if (testx == vert[i][0] && testy == vert[i][1]) {
@@ -135,76 +148,84 @@ static int pnpoly(int nvert, double vert[][2], double testx, double testy) {
     }
     return c;
 }
+
+
 /**
  * Return 1 if (x, y) lies in (the interior or boundary of) the image of the
- * HEALPix projection (in case proj=0) or in the image the rHEALPix projection 
+ * HEALPix projection (in case proj=0) or in the image the rHEALPix projection
  * (in case proj=1), and return 0 otherwise.
- * @param north_square the position of the north polar square (rHEALPix only)       
- * @param south_square the position of the south polar square (rHEALPix only)  
+ * @param north_square the position of the north polar square (rHEALPix only)
+ * @param south_square the position of the south polar square (rHEALPix only)
  **/
 int in_image(double x, double y, int proj, int north_square, int south_square) {
     if (proj == 0) {
-    	double healpixVertsJit[][2] = {
-    	    {-1.0*PI- EPS, PI/4.0},
-    	    {-3.0*PI/4.0, PI/2.0 + EPS},
-    	    {-1.0*PI/2.0, PI/4.0 + EPS},
-    	    {-1.0*PI/4.0, PI/2.0 + EPS},
-    	    {0.0, PI/4.0 + EPS},
-    	    {PI/4.0, PI/2.0 + EPS},
-    	    {PI/2.0, PI/4.0 + EPS},
-    	    {3.0*PI/4.0, PI/2.0 + EPS},
-    	    {PI+ EPS, PI/4.0},
-    	    {PI+ EPS, -1.0*PI/4.0},
-    	    {3.0*PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {PI/2.0, -1.0*PI/4.0 - EPS},
-    	    {PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {0.0, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {-1.0*PI/2.0, -1.0*PI/4.0 - EPS},
-    	    {-3.0*PI/4.0, -1.0*PI/2.0 - EPS},
-    	    {-1.0*PI - EPS, -1.0*PI/4.0}
-    	};
-    	return pnpoly((int)sizeof(healpixVertsJit)/
-    	              sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
+        double healpixVertsJit[][2] = {
+            {-1.0*PI- EPS, PI/4.0},
+            {-3.0*PI/4.0, PI/2.0 + EPS},
+            {-1.0*PI/2.0, PI/4.0 + EPS},
+            {-1.0*PI/4.0, PI/2.0 + EPS},
+            {0.0, PI/4.0 + EPS},
+            {PI/4.0, PI/2.0 + EPS},
+            {PI/2.0, PI/4.0 + EPS},
+            {3.0*PI/4.0, PI/2.0 + EPS},
+            {PI+ EPS, PI/4.0},
+            {PI+ EPS, -1.0*PI/4.0},
+            {3.0*PI/4.0, -1.0*PI/2.0 - EPS},
+            {PI/2.0, -1.0*PI/4.0 - EPS},
+            {PI/4.0, -1.0*PI/2.0 - EPS},
+            {0.0, -1.0*PI/4.0 - EPS},
+            {-1.0*PI/4.0, -1.0*PI/2.0 - EPS},
+            {-1.0*PI/2.0, -1.0*PI/4.0 - EPS},
+            {-3.0*PI/4.0, -1.0*PI/2.0 - EPS},
+            {-1.0*PI - EPS, -1.0*PI/4.0}
+        };
+        return pnpoly((int)sizeof(healpixVertsJit)/
+                      sizeof(healpixVertsJit[0]), healpixVertsJit, x, y);
     } else {
-    	double rhealpixVertsJit[][2] = {
-    	    {-1.0*PI - EPS, PI/4.0 + EPS},
-    	    {-1.0*PI + north_square*PI/2.0- EPS, PI/4.0 + EPS},
-    	    {-1.0*PI + north_square*PI/2.0- EPS, 3*PI/4.0 + EPS},
-    	    {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, 3*PI/4.0 + EPS},
-    	    {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, PI/4.0 + EPS},
-    	    {PI + EPS, PI/4.0 + EPS},
-    	    {PI + EPS, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -3.0*PI/4.0 - EPS},
-    	    {-1.0*PI + south_square*PI/2.0 - EPS, -3.0*PI/4.0 - EPS},
-    	    {-1.0*PI + south_square*PI/2.0 - EPS, -1.0*PI/4.0 - EPS},
-    	    {-1.0*PI - EPS, -1.0*PI/4.0 - EPS}};
-    	return pnpoly((int)sizeof(rhealpixVertsJit)/
-    	              sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
+        double rhealpixVertsJit[][2] = {
+            {-1.0*PI - EPS, PI/4.0 + EPS},
+            {-1.0*PI + north_square*PI/2.0- EPS, PI/4.0 + EPS},
+            {-1.0*PI + north_square*PI/2.0- EPS, 3*PI/4.0 + EPS},
+            {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, 3*PI/4.0 + EPS},
+            {-1.0*PI + (north_square + 1.0)*PI/2.0 + EPS, PI/4.0 + EPS},
+            {PI + EPS, PI/4.0 + EPS},
+            {PI + EPS, -1.0*PI/4.0 - EPS},
+            {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -1.0*PI/4.0 - EPS},
+            {-1.0*PI + (south_square + 1.0)*PI/2.0 + EPS, -3.0*PI/4.0 - EPS},
+            {-1.0*PI + south_square*PI/2.0 - EPS, -3.0*PI/4.0 - EPS},
+            {-1.0*PI + south_square*PI/2.0 - EPS, -1.0*PI/4.0 - EPS},
+            {-1.0*PI - EPS, -1.0*PI/4.0 - EPS}};
+        return pnpoly((int)sizeof(rhealpixVertsJit)/
+                      sizeof(rhealpixVertsJit[0]), rhealpixVertsJit, x, y);
     }
 }
+
+
 /**
  * Return the authalic latitude of latitude alpha (if inverse=0) or
- * return the approximate latitude of authalic latitude alpha (if inverse=1).    
- * P contains the relavent ellipsoid parameters. 
+ * return the approximate latitude of authalic latitude alpha (if inverse=1).
+ * P contains the relavent ellipsoid parameters.
  **/
 double auth_lat(PJ *P, double alpha, int inverse) {
+    struct pj_opaque *Q = P->opaque;
     if (inverse == 0) {
         /* Authalic latitude. */
         double q = pj_qsfn(sin(alpha), P->e, 1.0 - P->es);
-        double qp = P->qp; 
-	    double ratio = q/qp;
-	    if (fabsl(ratio) > 1) {
-	        /* Rounding error. */
-	        ratio = pj_sign(ratio);
-	    }
-	    return asin(ratio);
+        double qp = Q->qp;
+        double ratio = q/qp;
+
+        if (fabsl(ratio) > 1) {
+            /* Rounding error. */
+            ratio = pj_sign(ratio);
+        }
+        return asin(ratio);
     } else {
         /* Approximation to inverse authalic latitude. */
-        return pj_authlat(alpha, P->apa);
+        return pj_authlat(alpha, Q->apa);
     }
 }
+
+
 /**
  * Return the HEALPix projection of the longitude-latitude point lp on
  * the unit sphere.
@@ -214,51 +235,57 @@ XY healpix_sphere(LP lp) {
     double phi = lp.phi;
     double phi0 = asin(2.0/3.0);
     XY xy;
+
     /* equatorial region */
     if ( fabsl(phi) <= phi0) {
-	    xy.x = lam;
-	    xy.y = 3.0*PI/8.0*sin(phi);
+        xy.x = lam;
+        xy.y = 3.0*PI/8.0*sin(phi);
     } else {
-	    double lamc;
-	    double sigma = sqrt(3.0*(1 - fabsl(sin(phi))));
-	    double cn = floor(2*lam / PI + 2);
-	    if (cn >= 4) {
-	        cn = 3;
-	    }
-	    lamc = -3*PI/4 + (PI/2)*cn;
-	    xy.x = lamc + (lam - lamc)*sigma;
-	    xy.y = pj_sign(phi)*PI/4*(2 - sigma);
+        double lamc;
+        double sigma = sqrt(3.0*(1 - fabsl(sin(phi))));
+        double cn = floor(2*lam / PI + 2);
+        if (cn >= 4) {
+            cn = 3;
+        }
+        lamc = -3*PI/4 + (PI/2)*cn;
+        xy.x = lamc + (lam - lamc)*sigma;
+        xy.y = pj_sign(phi)*PI/4*(2 - sigma);
     }
     return xy;
 }
+
+
 /**
- * Return the inverse of healpix_sphere(). 
+ * Return the inverse of healpix_sphere().
 **/
 LP healpix_sphere_inverse(XY xy) {
-    LP lp; 
+    LP lp;
     double x = xy.x;
     double y = xy.y;
     double y0 = PI/4.0;
+
     /* Equatorial region. */
     if (fabsl(y) <= y0) {
-	    lp.lam = x;
-	    lp.phi = asin(8.0*y/(3.0*PI));	
+        lp.lam = x;
+        lp.phi = asin(8.0*y/(3.0*PI));
     } else if (fabsl(y) < PI/2.0) {
-	    double cn = floor(2.0*x/PI + 2.0);
+        double cn = floor(2.0*x/PI + 2.0);
         double xc, tau;
-	    if (cn >= 4) {
-	        cn = 3;
-	    }
-	    xc = -3.0*PI/4.0 + (PI/2.0)*cn;
-	    tau = 2.0 - 4.0*fabsl(y)/PI;
-	    lp.lam = xc + (x - xc)/tau;	
-	    lp.phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
+        if (cn >= 4) {
+            cn = 3;
+        }
+        xc = -3.0*PI/4.0 + (PI/2.0)*cn;
+        tau = 2.0 - 4.0*fabsl(y)/PI;
+        lp.lam = xc + (x - xc)/tau;
+        lp.phi = pj_sign(y)*asin(1.0 - pow(tau , 2.0)/3.0);
     } else {
-	    lp.lam = -1.0*PI;
-	    lp.phi = pj_sign(y)*PI/2.0;
+        lp.lam = -1.0*PI;
+        lp.phi = pj_sign(y)*PI/2.0;
     }
     return (lp);
 }
+
+
 /**
  * Return the vector sum a + b, where a and b are 2-dimensional vectors.
  * @param ret holds a + b.
@@ -266,21 +293,25 @@ LP healpix_sphere_inverse(XY xy) {
 static void vector_add(double a[2], double b[2], double *ret) {
     int i;
     for(i = 0; i < 2; i++) {
-	    ret[i] = a[i] + b[i];
+        ret[i] = a[i] + b[i];
     }
 }
+
+
 /**
  * Return the vector difference a - b, where a and b are 2-dimensional vectors.
  * @param ret holds a - b.
  **/
 static void vector_sub(double a[2], double b[2], double*ret) {
     int i;
-    for(i = 0; i < 2; i++) { 
-	    ret[i] = a[i] - b[i];
+    for(i = 0; i < 2; i++) {
+        ret[i] = a[i] - b[i];
     }
 }
+
+
 /**
- * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and 
+ * Return the 2 x 1 matrix product a*b, where a is a 2 x 2 matrix and
  * b is a 2 x 1 matrix.
  * @param ret holds a*b.
  **/
@@ -288,105 +319,110 @@ static void dot_product(double a[2][2], double b[2], double *ret) {
     int i, j;
     int length = 2;
     for(i = 0; i < length; i++) {
-	    ret[i] = 0;
-	    for(j = 0; j < length; j++) {
-	        ret[i] += a[i][j]*b[j];
-	    }
+        ret[i] = 0;
+        for(j = 0; j < length; j++) {
+            ret[i] += a[i][j]*b[j];
+        }
     }
 }
+
+
 /**
- * Return the number of the polar cap, the pole point coordinates, and 
+ * Return the number of the polar cap, the pole point coordinates, and
  * the region that (x, y) lies in.
- * If inverse=0, then assume (x,y) lies in the image of the HEALPix 
+ * If inverse=0, then assume (x,y) lies in the image of the HEALPix
  * projection of the unit sphere.
- * If inverse=1, then assume (x,y) lies in the image of the 
+ * If inverse=1, then assume (x,y) lies in the image of the
  * (north_square, south_square)-rHEALPix projection of the unit sphere.
  **/
 static CapMap get_cap(double x, double y, int north_square, int south_square,
                       int inverse) {
     CapMap capmap;
     double c;
+
     capmap.x = x;
     capmap.y = y;
     if (inverse == 0) {
-    	if (y > PI/4.0) {
-    	    capmap.region = north;
-    	    c = PI/2.0; 
-    	} else if (y < -1*PI/4.0) {
-    	    capmap.region = south;
-    	    c = -1*PI/2.0;
-    	} else {
-    	    capmap.region = equatorial;
-    	    capmap.cn = 0;
-    	    return capmap;
-    	}
-    	/* polar region */
-    	if (x < -1*PI/2.0) {
-    	    capmap.cn = 0;
-    	    capmap.x = (-1*3.0*PI/4.0);
-    	    capmap.y = c;
-    	} else if (x >= -1*PI/2.0 && x < 0) {
-    	    capmap.cn = 1;
-    	    capmap.x = -1*PI/4.0;
-    	    capmap.y = c;
-    	} else if (x >= 0 && x < PI/2.0) {
-    	    capmap.cn = 2;
-    	    capmap.x = PI/4.0;
-    	    capmap.y = c;
-    	} else {
-    	    capmap.cn = 3;
-    	    capmap.x = 3.0*PI/4.0;
-    	    capmap.y = c;
-    	}
-    	return capmap;
+        if (y > PI/4.0) {
+            capmap.region = north;
+            c = PI/2.0;
+        } else if (y < -1*PI/4.0) {
+            capmap.region = south;
+            c = -1*PI/2.0;
+        } else {
+            capmap.region = equatorial;
+            capmap.cn = 0;
+            return capmap;
+        }
+        /* polar region */
+        if (x < -1*PI/2.0) {
+            capmap.cn = 0;
+            capmap.x = (-1*3.0*PI/4.0);
+            capmap.y = c;
+        } else if (x >= -1*PI/2.0 && x < 0) {
+            capmap.cn = 1;
+            capmap.x = -1*PI/4.0;
+            capmap.y = c;
+        } else if (x >= 0 && x < PI/2.0) {
+            capmap.cn = 2;
+            capmap.x = PI/4.0;
+            capmap.y = c;
+        } else {
+            capmap.cn = 3;
+            capmap.x = 3.0*PI/4.0;
+            capmap.y = c;
+        }
+        return capmap;
     } else {
-    	double eps;
-    	if (y > PI/4.0) {
-    	    capmap.region = north;
-    	    capmap.x = (-3.0*PI/4.0 + north_square*PI/2.0); 
-    	    capmap.y = PI/2.0;
-    	    x = x - north_square*PI/2.0;
-    	} else if (y < -1*PI/4.0) {
-    	    capmap.region = south;
-    	    capmap.x = (-3.0*PI/4.0 + south_square*PI/2); 
-    	    capmap.y = -1*PI/2.0;
-    	    x = x - south_square*PI/2.0;
-    	} else {
-    	    capmap.region = equatorial;
-    	    capmap.cn = 0;
-    	    return capmap;
-    	}
-    	/* Polar Region, find the HEALPix polar cap number that 
-    	   x, y moves to when rHEALPix polar square is disassembled. */
-    	eps = 1e-15; /* Kludge.  Fuzz to avoid some rounding errors. */
-    	if (capmap.region == north) {
-    	    if (y >= -1*x - PI/4.0 - eps && y < x + 5.0*PI/4.0 - eps) {
-    		    capmap.cn = (north_square + 1) % 4;
-    	    } else if (y > -1*x -1*PI/4.0 + eps && y >= x + 5.0*PI/4.0 - eps) {
-    		    capmap.cn = (north_square + 2) % 4;
-    	    } else if (y <= -1*x -1*PI/4.0 + eps && y > x + 5.0*PI/4.0 + eps) {
-    		    capmap.cn = (north_square + 3) % 4;
-    	    } else {
-    		    capmap.cn = north_square;
-    	    }
-    	} else if (capmap.region == south) {
-    	    if (y <= x + PI/4.0 + eps && y > -1*x - 5.0*PI/4 + eps) {
-    		    capmap.cn = (south_square + 1) % 4;
-    	    } else if (y < x + PI/4.0 - eps && y <= -1*x - 5.0*PI/4.0 + eps) {
-    		    capmap.cn = (south_square + 2) % 4;
-    	    } else if (y >= x + PI/4.0 - eps && y < -1*x - 5.0*PI/4.0 - eps) {
-    		    capmap.cn = (south_square + 3) % 4;
-    	    } else {
-    		    capmap.cn = south_square;
-    	    }
-    	}
-    	return capmap;
+        double eps;
+        if (y > PI/4.0) {
+            capmap.region = north;
+            capmap.x = (-3.0*PI/4.0 + north_square*PI/2.0);
+            capmap.y = PI/2.0;
+            x = x - north_square*PI/2.0;
+        } else if (y < -1*PI/4.0) {
+            capmap.region = south;
+            capmap.x = (-3.0*PI/4.0 + south_square*PI/2);
+            capmap.y = -1*PI/2.0;
+            x = x - south_square*PI/2.0;
+        } else {
+            capmap.region = equatorial;
+            capmap.cn = 0;
+            return capmap;
+        }
+        /* Polar Region, find the HEALPix polar cap number that
+           x, y moves to when rHEALPix polar square is disassembled. */
+        eps = 1e-15; /* Kludge.  Fuzz to avoid some rounding errors. */
+        if (capmap.region == north) {
+            if (y >= -1*x - PI/4.0 - eps && y < x + 5.0*PI/4.0 - eps) {
+                capmap.cn = (north_square + 1) % 4;
+            } else if (y > -1*x -1*PI/4.0 + eps && y >= x + 5.0*PI/4.0 - eps) {
+                capmap.cn = (north_square + 2) % 4;
+            } else if (y <= -1*x -1*PI/4.0 + eps && y > x + 5.0*PI/4.0 + eps) {
+                capmap.cn = (north_square + 3) % 4;
+            } else {
+                capmap.cn = north_square;
+            }
+        } else if (capmap.region == south) {
+            if (y <= x + PI/4.0 + eps && y > -1*x - 5.0*PI/4 + eps) {
+                capmap.cn = (south_square + 1) % 4;
+            } else if (y < x + PI/4.0 - eps && y <= -1*x - 5.0*PI/4.0 + eps) {
+                capmap.cn = (south_square + 2) % 4;
+            } else if (y >= x + PI/4.0 - eps && y < -1*x - 5.0*PI/4.0 - eps) {
+                capmap.cn = (south_square + 3) % 4;
+            } else {
+                capmap.cn = south_square;
+            }
+        }
+        return capmap;
     }
 }
+
+
 /**
- * Rearrange point (x, y) in the HEALPix projection by 
+ * Rearrange point (x, y) in the HEALPix projection by
  * combining the polar caps into two polar squares.
- * Put the north polar square in position north_square and 
+ * Put the north polar square in position north_square and
  * the south polar square in position south_square.
  * If inverse=1, then uncombine the polar caps.
  * @param north_square integer between 0 and 3.
@@ -400,173 +436,352 @@ static XY combine_caps(double x, double y, int north_square, int south_square,
     double vector[2];
     double v_min_c[2];
     double ret_dot[2];
+
     CapMap capmap = get_cap(x, y, north_square, south_square, inverse);
     if (capmap.region == equatorial) {
-	    xy.x = capmap.x;
-	    xy.y = capmap.y;
-	    return xy;
+        xy.x = capmap.x;
+        xy.y = capmap.y;
+        return xy;
     }
     v[0] = x;
     v[1] = y;
     if (inverse == 0) {
-        /* Rotate (x, y) about its polar cap tip and then translate it to  
+        /* Rotate (x, y) about its polar cap tip and then translate it to
            north_square or south_square. */
-    	int pole = 0;
-    	double (*tmpRot)[2];
-    	double c[2] = {capmap.x, capmap.y};
-    	if (capmap.region == north) {
-    	    pole = north_square;
-    	    a[0] =  (-3.0*PI/4.0 + pole*PI/2);
-    	    a[1] =  (PI/2.0 + pole*0);
-    	    tmpRot = rot[get_rotate_index(capmap.cn - pole)];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	} else {
-    	    pole = south_square;
-    	    a[0] =  (-3.0*PI/4.0 + pole*PI/2);
-    	    a[1] =  (PI/-2.0 + pole*0);
-    	    tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	}
-    	xy.x = vector[0];
-    	xy.y = vector[1];
-    	return xy;
+        int pole = 0;
+        double (*tmpRot)[2];
+        double c[2] = {capmap.x, capmap.y};
+        if (capmap.region == north) {
+            pole = north_square;
+            a[0] =  (-3.0*PI/4.0 + pole*PI/2);
+            a[1] =  (PI/2.0 + pole*0);
+            tmpRot = rot[get_rotate_index(capmap.cn - pole)];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        } else {
+            pole = south_square;
+            a[0] =  (-3.0*PI/4.0 + pole*PI/2);
+            a[1] =  (PI/-2.0 + pole*0);
+            tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        }
+        xy.x = vector[0];
+        xy.y = vector[1];
+        return xy;
     } else {
         /* Inverse function.
          Unrotate (x, y) and then translate it back. */
-    	int pole = 0;
-    	double (*tmpRot)[2];
-    	double c[2] = {capmap.x, capmap.y};
-    	/* disassemble */
-    	if (capmap.region == north) {
-    	    pole = north_square;
-    	    a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
-    	    a[1] =  (PI/2.0 + capmap.cn*0);
-    	    tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	} else {
-    	    pole = south_square;
-    	    a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
-    	    a[1] =  (PI/-2.0 + capmap.cn*0);
-    	    tmpRot = rot[get_rotate_index(capmap.cn - pole)];
-    	    vector_sub(v, c, v_min_c);	
-    	    dot_product(tmpRot, v_min_c, ret_dot);
-    	    vector_add(ret_dot, a, vector);
-    	}
-    	xy.x = vector[0];
-    	xy.y = vector[1];
-    	return xy;
+        int pole = 0;
+        double (*tmpRot)[2];
+        double c[2] = {capmap.x, capmap.y};
+        /* disassemble */
+        if (capmap.region == north) {
+            pole = north_square;
+            a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
+            a[1] =  (PI/2.0 + capmap.cn*0);
+            tmpRot = rot[get_rotate_index(-1*(capmap.cn - pole))];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        } else {
+            pole = south_square;
+            a[0] =  (-3.0*PI/4.0 + capmap.cn*PI/2);
+            a[1] =  (PI/-2.0 + capmap.cn*0);
+            tmpRot = rot[get_rotate_index(capmap.cn - pole)];
+            vector_sub(v, c, v_min_c);
+            dot_product(tmpRot, v_min_c, ret_dot);
+            vector_add(ret_dot, a, vector);
+        }
+        xy.x = vector[0];
+        xy.y = vector[1];
+        return xy;
     }
 }
-FORWARD(s_healpix_forward); /* sphere  */
+
+
+static XY s_healpix_forward(LP lp, PJ *P) { /* sphere  */
     (void) P;
-    (void) xy;
     return healpix_sphere(lp);
 }
-FORWARD(e_healpix_forward); /* ellipsoid  */
-    (void) xy;
+
+
+static XY e_healpix_forward(LP lp, PJ *P) { /* ellipsoid  */
     lp.phi = auth_lat(P, lp.phi, 0);
     return healpix_sphere(lp);
 }
-INVERSE(s_healpix_inverse); /* sphere */
+
+
+static LP s_healpix_inverse(XY xy, PJ *P) { /* sphere */
+    LP lp = {0.0,0.0};
+
     /* Check whether (x, y) lies in the HEALPix image */
     if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
     return healpix_sphere_inverse(xy);
 }
-INVERSE(e_healpix_inverse); /* ellipsoid */
+
+
+static LP e_healpix_inverse(XY xy, PJ *P) { /* ellipsoid */
+    LP lp = {0.0,0.0};
+
     /* Check whether (x, y) lies in the HEALPix image. */
     if (in_image(xy.x, xy.y, 0, 0, 0) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
     lp = healpix_sphere_inverse(xy);
     lp.phi = auth_lat(P, lp.phi, 1);
-    return (lp);
+    return lp;
 }
-FORWARD(s_rhealpix_forward); /* sphere */
-    xy = healpix_sphere(lp);
-    return combine_caps(xy.x, xy.y, P->north_square, P->south_square, 0);
+
+
+static XY s_rhealpix_forward(LP lp, PJ *P) { /* sphere */
+    struct pj_opaque *Q = P->opaque;
+
+    XY xy = healpix_sphere(lp);
+    return combine_caps(xy.x, xy.y, Q->north_square, Q->south_square, 0);
 }
-FORWARD(e_rhealpix_forward); /* ellipsoid */
+
+
+static XY e_rhealpix_forward(LP lp, PJ *P) { /* ellipsoid */
+    struct pj_opaque *Q = P->opaque;
+
     lp.phi = auth_lat(P, lp.phi, 0);
-    xy = healpix_sphere(lp);
-    return combine_caps(xy.x, xy.y, P->north_square, P->south_square, 0);
+    XY xy = healpix_sphere(lp);
+    return combine_caps(xy.x, xy.y, Q->north_square, Q->south_square, 0);
 }
-INVERSE(s_rhealpix_inverse); /* sphere */
+
+
+static LP s_rhealpix_inverse(XY xy, PJ *P) { /* sphere */
+    struct pj_opaque *Q = P->opaque;
+    LP lp = {0.0,0.0};
+
     /* Check whether (x, y) lies in the rHEALPix image. */
-    if (in_image(xy.x, xy.y, 1, P->north_square, P->south_square) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+    if (in_image(xy.x, xy.y, 1, Q->north_square, Q->south_square) == 0) {
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
-    xy = combine_caps(xy.x, xy.y, P->north_square, P->south_square, 1);
+    xy = combine_caps(xy.x, xy.y, Q->north_square, Q->south_square, 1);
     return healpix_sphere_inverse(xy);
 }
-INVERSE(e_rhealpix_inverse); /* ellipsoid */
+
+
+static LP e_rhealpix_inverse(XY xy, PJ *P) { /* ellipsoid */
+    struct pj_opaque *Q = P->opaque;
+    LP lp = {0.0,0.0};
+
     /* Check whether (x, y) lies in the rHEALPix image. */
-    if (in_image(xy.x, xy.y, 1, P->north_square, P->south_square) == 0) {
-	    lp.lam = HUGE_VAL;
-	    lp.phi = HUGE_VAL;
-	    pj_ctx_set_errno(P->ctx, -15);
-	    return lp;
+    if (in_image(xy.x, xy.y, 1, Q->north_square, Q->south_square) == 0) {
+        lp.lam = HUGE_VAL;
+        lp.phi = HUGE_VAL;
+        pj_ctx_set_errno(P->ctx, -15);
+        return lp;
     }
-    xy = combine_caps(xy.x, xy.y, P->north_square, P->south_square, 1);
+    xy = combine_caps(xy.x, xy.y, Q->north_square, Q->south_square, 1);
     lp = healpix_sphere_inverse(xy);
     lp.phi = auth_lat(P, lp.phi, 1);
     return lp;
 }
-FREEUP;
-	if (P) {
-		if (P->apa)
-			pj_dalloc(P->apa);
-		pj_dalloc(P);
-	}
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    if (P->opaque->apa)
+        pj_dealloc(P->opaque->apa);
+
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-ENTRY1(healpix, apa)
+
+
+PJ *PROJECTION(healpix) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
     if (P->es) {
-        P->apa = pj_authset(P->es); /* For auth_lat(). */
-        P->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
-    	P->a = P->a*sqrt(0.5*P->qp); /* Set P->a to authalic radius. */
+        Q->apa = pj_authset(P->es); /* For auth_lat(). */
+        Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
+        P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
         P->ra = 1.0/P->a;
-    	P->fwd = e_healpix_forward;
-    	P->inv = e_healpix_inverse; 
+        P->fwd = e_healpix_forward;
+        P->inv = e_healpix_inverse;
     } else {
-    	P->fwd = s_healpix_forward;
-    	P->inv = s_healpix_inverse; 
+        P->fwd = s_healpix_forward;
+        P->inv = s_healpix_inverse;
     }
-ENDENTRY(P)
-ENTRY1(rhealpix, apa)
-    P->north_square = pj_param(P->ctx, P->params,"inorth_square").i;
-    P->south_square = pj_param(P->ctx, P->params,"isouth_square").i;
+
+    return P;
+}
+
+
+PJ *PROJECTION(rhealpix) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->north_square = pj_param(P->ctx, P->params,"inorth_square").i;
+    Q->south_square = pj_param(P->ctx, P->params,"isouth_square").i;
+
     /* Check for valid north_square and south_square inputs. */
-    if (P->north_square < 0 || P->north_square > 3) {
-	    E_ERROR(-47);
+    if (Q->north_square < 0 || Q->north_square > 3) {
+        E_ERROR(-47);
     }
-    if (P->south_square < 0 || P->south_square > 3) {
-	    E_ERROR(-47);
+    if (Q->south_square < 0 || Q->south_square > 3) {
+        E_ERROR(-47);
     }
     if (P->es) {
-        P->apa = pj_authset(P->es); /* For auth_lat(). */
-        P->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
-	    P->a = P->a*sqrt(0.5*P->qp); /* Set P->a to authalic radius. */
+        Q->apa = pj_authset(P->es); /* For auth_lat(). */
+        Q->qp = pj_qsfn(1.0, P->e, P->one_es); /* For auth_lat(). */
+        P->a = P->a*sqrt(0.5*Q->qp); /* Set P->a to authalic radius. */
         P->ra = 1.0/P->a;
-	    P->fwd = e_rhealpix_forward;
-	    P->inv = e_rhealpix_inverse; 
+        P->fwd = e_rhealpix_forward;
+        P->inv = e_rhealpix_inverse;
     } else {
-	    P->fwd = s_rhealpix_forward;
-	    P->inv = s_rhealpix_inverse; 
+        P->fwd = s_rhealpix_forward;
+        P->inv = s_rhealpix_inverse;
     }
-ENDENTRY(P)
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_healpix_selftest (void) {return 0;}
+#else
+
+int pj_healpix_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=healpix   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=healpix   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222390.10394923863,  130406.58866448226},
+        { 222390.10394923863, -130406.58866448054},
+        {-222390.10394923863,  130406.58866448226},
+        {-222390.10394923863, -130406.58866448054},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223402.14425527418,  131588.04444199943},
+        { 223402.14425527418, -131588.04444199943},
+        {-223402.14425527418,  131588.04444199943},
+        {-223402.14425527418, -131588.04444199943},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017986411845524453,  0.00076679453057823619},
+        { 0.0017986411845524453, -0.00076679453057823619},
+        {-0.0017986411845524453,  0.00076679453057823619},
+        {-0.0017986411845524453, -0.00076679453057823619},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931097838226,  0.00075990887733981202},
+        { 0.0017904931097838226, -0.00075990887733981202},
+        {-0.0017904931097838226,  0.00075990887733981202},
+        {-0.0017904931097838226, -0.00075990887733981202},
+    };
+
+    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);
+}
+
+
+#endif
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_rhealpix_selftest (void) {return 0;}
+#else
+
+int pj_rhealpix_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=rhealpix   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=rhealpix   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222390.10394923863,  130406.58866448226},
+        { 222390.10394923863, -130406.58866448054},
+        {-222390.10394923863,  130406.58866448226},
+        {-222390.10394923863, -130406.58866448054},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223402.14425527418,  131588.04444199943},
+        { 223402.14425527418, -131588.04444199943},
+        {-223402.14425527418,  131588.04444199943},
+        {-223402.14425527418, -131588.04444199943},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017986411845524453,  0.00076679453057823619},
+        { 0.0017986411845524453, -0.00076679453057823619},
+        {-0.0017986411845524453,  0.00076679453057823619},
+        {-0.0017986411845524453, -0.00076679453057823619},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931097838226,  0.00075990887733981202},
+        { 0.0017904931097838226, -0.00075990887733981202},
+        {-0.0017904931097838226,  0.00075990887733981202},
+        {-0.0017904931097838226, -0.00075990887733981202},
+    };
+
+    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);
+}
+
+
+#endif
diff --git a/src/PJ_igh.c b/src/PJ_igh.c
index 4155c856..f82f4f3f 100644
--- a/src/PJ_igh.c
+++ b/src/PJ_igh.c
@@ -1,12 +1,9 @@
-#define PROJ_PARMS__ \
-        struct PJconsts* pj[12]; \
-        double dy0;
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(igh, "Interrupted Goode Homolosine") "\n\tPCyl, Sph.";
-	C_NAMESPACE PJ
-*pj_sinu(PJ *), *pj_moll(PJ *);
 
+C_NAMESPACE PJ *pj_sinu(PJ *), *pj_moll(PJ *);
 static const double d4044118 = (40 + 44/60. + 11.8/3600.) * DEG_TO_RAD; // 40d 44' 11.8" [degrees]
 
 static const double d10  =  10 * DEG_TO_RAD;
@@ -24,106 +21,137 @@ static const double d180 = 180 * DEG_TO_RAD;
 
 static const double EPSLN = 1.e-10; // allow a little 'slack' on zone edge positions
 
-FORWARD(s_forward); /* spheroid */
-        int z;
-        if (lp.phi >=  d4044118) {          // 1|2
-          z = (lp.lam <= -d40 ? 1: 2);
-        }
-        else if (lp.phi >=  0) {            // 3|4
-          z = (lp.lam <= -d40 ? 3: 4);
-        }
-        else if (lp.phi >= -d4044118) {     // 5|6|7|8
-               if (lp.lam <= -d100) z =  5; // 5
-          else if (lp.lam <=  -d20) z =  6; // 6
-          else if (lp.lam <=   d80) z =  7; // 7
-          else z = 8;                       // 8
-        }
-        else {                              // 9|10|11|12
-               if (lp.lam <= -d100) z =  9; // 9
-          else if (lp.lam <=  -d20) z = 10; // 10
-          else if (lp.lam <=   d80) z = 11; // 11
-          else z = 12;                      // 12
-        }
+struct pj_opaque {
+    struct PJconsts* pj[12]; \
+    double dy0;
+};
 
-        lp.lam -= P->pj[z-1]->lam0;
-        xy = P->pj[z-1]->fwd(lp, P->pj[z-1]);
-        xy.x += P->pj[z-1]->x0;
-        xy.y += P->pj[z-1]->y0;
 
-        return (xy);
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    int z;
+
+    if (lp.phi >=  d4044118) {          // 1|2
+      z = (lp.lam <= -d40 ? 1: 2);
+    }
+    else if (lp.phi >=  0) {            // 3|4
+      z = (lp.lam <= -d40 ? 3: 4);
+    }
+    else if (lp.phi >= -d4044118) {     // 5|6|7|8
+           if (lp.lam <= -d100) z =  5; // 5
+      else if (lp.lam <=  -d20) z =  6; // 6
+      else if (lp.lam <=   d80) z =  7; // 7
+      else z = 8;                       // 8
+    }
+    else {                              // 9|10|11|12
+           if (lp.lam <= -d100) z =  9; // 9
+      else if (lp.lam <=  -d20) z = 10; // 10
+      else if (lp.lam <=   d80) z = 11; // 11
+      else z = 12;                      // 12
+    }
+
+    lp.lam -= Q->pj[z-1]->lam0;
+    xy = Q->pj[z-1]->fwd(lp, Q->pj[z-1]);
+    xy.x += Q->pj[z-1]->x0;
+    xy.y += Q->pj[z-1]->y0;
+
+    return xy;
 }
-INVERSE(s_inverse); /* spheroid */
-        const double y90 = P->dy0 + sqrt(2); // lt=90 corresponds to y=y0+sqrt(2)
-
-        int z = 0;
-        if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) // 0
-          z = 0;
-        else if (xy.y >=  d4044118)       // 1|2
-          z = (xy.x <= -d40? 1: 2);
-        else if (xy.y >=  0)              // 3|4
-          z = (xy.x <= -d40? 3: 4);
-        else if (xy.y >= -d4044118) {     // 5|6|7|8
-               if (xy.x <= -d100) z =  5; // 5
-          else if (xy.x <=  -d20) z =  6; // 6
-          else if (xy.x <=   d80) z =  7; // 7
-          else z = 8;                     // 8
-        }
-        else {                            // 9|10|11|12
-               if (xy.x <= -d100) z =  9; // 9
-          else if (xy.x <=  -d20) z = 10; // 10
-          else if (xy.x <=   d80) z = 11; // 11
-          else z = 12;                    // 12
-        }
 
-        if (z)
-        {
-          int ok = 0;
-
-          xy.x -= P->pj[z-1]->x0;
-          xy.y -= P->pj[z-1]->y0;
-          lp = P->pj[z-1]->inv(xy, P->pj[z-1]);
-          lp.lam += P->pj[z-1]->lam0;
-
-          switch (z) {
-            case  1: ok = (lp.lam >= -d180-EPSLN && lp.lam <=  -d40+EPSLN) ||
-                         ((lp.lam >=  -d40-EPSLN && lp.lam <=  -d10+EPSLN) &&
-                          (lp.phi >=   d60-EPSLN && lp.phi <=   d90+EPSLN)); break;
-            case  2: ok = (lp.lam >=  -d40-EPSLN && lp.lam <=  d180+EPSLN) ||
-                         ((lp.lam >= -d180-EPSLN && lp.lam <= -d160+EPSLN) &&
-                          (lp.phi >=   d50-EPSLN && lp.phi <=   d90+EPSLN)) ||
-                         ((lp.lam >=  -d50-EPSLN && lp.lam <=  -d40+EPSLN) &&
-                          (lp.phi >=   d60-EPSLN && lp.phi <=   d90+EPSLN)); break;
-            case  3: ok = (lp.lam >= -d180-EPSLN && lp.lam <=  -d40+EPSLN); break;
-            case  4: ok = (lp.lam >=  -d40-EPSLN && lp.lam <=  d180+EPSLN); break;
-            case  5: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break;
-            case  6: ok = (lp.lam >= -d100-EPSLN && lp.lam <=  -d20+EPSLN); break;
-            case  7: ok = (lp.lam >=  -d20-EPSLN && lp.lam <=   d80+EPSLN); break;
-            case  8: ok = (lp.lam >=   d80-EPSLN && lp.lam <=  d180+EPSLN); break;
-            case  9: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break;
-            case 10: ok = (lp.lam >= -d100-EPSLN && lp.lam <=  -d20+EPSLN); break;
-            case 11: ok = (lp.lam >=  -d20-EPSLN && lp.lam <=   d80+EPSLN); break;
-            case 12: ok = (lp.lam >=   d80-EPSLN && lp.lam <=  d180+EPSLN); break;
-          }
-
-          z = (!ok? 0: z); // projectable?
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    const double y90 = Q->dy0 + sqrt(2); // lt=90 corresponds to y=y0+sqrt(2)
+
+    int z = 0;
+    if (xy.y > y90+EPSLN || xy.y < -y90+EPSLN) // 0
+      z = 0;
+    else if (xy.y >=  d4044118)       // 1|2
+      z = (xy.x <= -d40? 1: 2);
+    else if (xy.y >=  0)              // 3|4
+      z = (xy.x <= -d40? 3: 4);
+    else if (xy.y >= -d4044118) {     // 5|6|7|8
+           if (xy.x <= -d100) z =  5; // 5
+      else if (xy.x <=  -d20) z =  6; // 6
+      else if (xy.x <=   d80) z =  7; // 7
+      else z = 8;                     // 8
+    }
+    else {                            // 9|10|11|12
+           if (xy.x <= -d100) z =  9; // 9
+      else if (xy.x <=  -d20) z = 10; // 10
+      else if (xy.x <=   d80) z = 11; // 11
+      else z = 12;                    // 12
+    }
+
+    if (z) {
+        int ok = 0;
+
+        xy.x -= Q->pj[z-1]->x0;
+        xy.y -= Q->pj[z-1]->y0;
+        lp = Q->pj[z-1]->inv(xy, Q->pj[z-1]);
+        lp.lam += Q->pj[z-1]->lam0;
+
+        switch (z) {
+        case  1: ok = (lp.lam >= -d180-EPSLN && lp.lam <=  -d40+EPSLN) ||
+                     ((lp.lam >=  -d40-EPSLN && lp.lam <=  -d10+EPSLN) &&
+                      (lp.phi >=   d60-EPSLN && lp.phi <=   d90+EPSLN)); break;
+        case  2: ok = (lp.lam >=  -d40-EPSLN && lp.lam <=  d180+EPSLN) ||
+                     ((lp.lam >= -d180-EPSLN && lp.lam <= -d160+EPSLN) &&
+                      (lp.phi >=   d50-EPSLN && lp.phi <=   d90+EPSLN)) ||
+                     ((lp.lam >=  -d50-EPSLN && lp.lam <=  -d40+EPSLN) &&
+                      (lp.phi >=   d60-EPSLN && lp.phi <=   d90+EPSLN)); break;
+        case  3: ok = (lp.lam >= -d180-EPSLN && lp.lam <=  -d40+EPSLN); break;
+        case  4: ok = (lp.lam >=  -d40-EPSLN && lp.lam <=  d180+EPSLN); break;
+        case  5: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break;
+        case  6: ok = (lp.lam >= -d100-EPSLN && lp.lam <=  -d20+EPSLN); break;
+        case  7: ok = (lp.lam >=  -d20-EPSLN && lp.lam <=   d80+EPSLN); break;
+        case  8: ok = (lp.lam >=   d80-EPSLN && lp.lam <=  d180+EPSLN); break;
+        case  9: ok = (lp.lam >= -d180-EPSLN && lp.lam <= -d100+EPSLN); break;
+        case 10: ok = (lp.lam >= -d100-EPSLN && lp.lam <=  -d20+EPSLN); break;
+        case 11: ok = (lp.lam >=  -d20-EPSLN && lp.lam <=   d80+EPSLN); break;
+        case 12: ok = (lp.lam >=   d80-EPSLN && lp.lam <=  d180+EPSLN); break;
         }
-     // if (!z) pj_errno = -15; // invalid x or y
-        if (!z) lp.lam = HUGE_VAL;
-        if (!z) lp.phi = HUGE_VAL;
-        return (lp);
+      z = (!ok? 0: z); // projectable?
+    }
+
+    if (!z) lp.lam = HUGE_VAL;
+    if (!z) lp.phi = HUGE_VAL;
+
+    return lp;
 }
-FREEUP;
-        if (P) {
-                int i;
-                for (i = 0; i < 12; ++i)
-		{
-			if (P->pj[i]) 
-				(*(P->pj[i]->pfree))(P->pj[i]); 
-                }
-                pj_dalloc(P);
-        }
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    int i;
+    for (i = 0; i < 12; ++i) {
+        if (P->opaque->pj[i])
+            pj_dealloc(P->opaque->pj[i]);
+    }
+
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-ENTRY0(igh)
+
+
+PJ *PROJECTION(igh) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
 /*
   Zones:
 
@@ -145,46 +173,90 @@ ENTRY0(igh)
 */
 
 #define SETUP(n, proj, x_0, y_0, lon_0) \
-    if (!(P->pj[n-1] = pj_##proj(0))) E_ERROR_0; \
-    if (!(P->pj[n-1] = pj_##proj(P->pj[n-1]))) E_ERROR_0; \
-    P->pj[n-1]->x0 = x_0; \
-    P->pj[n-1]->y0 = y_0; \
-    P->pj[n-1]->lam0 = lon_0; 
+    if (!(Q->pj[n-1] = pj_##proj(0))) E_ERROR_0; \
+    if (!(Q->pj[n-1] = pj_##proj(Q->pj[n-1]))) E_ERROR_0; \
+    Q->pj[n-1]->x0 = x_0; \
+    Q->pj[n-1]->y0 = y_0; \
+    Q->pj[n-1]->lam0 = lon_0;
 
-        LP lp = { 0, d4044118 };
-        XY xy1;
-        XY xy3;
+    LP lp = { 0, d4044118 };
+    XY xy1;
+    XY xy3;
 
-        // sinusoidal zones
-        SETUP(3, sinu, -d100, 0, -d100);
-        SETUP(4, sinu,   d30, 0,   d30);
-        SETUP(5, sinu, -d160, 0, -d160);
-        SETUP(6, sinu,  -d60, 0,  -d60);
-        SETUP(7, sinu,   d20, 0,   d20);
-        SETUP(8, sinu,  d140, 0,  d140);
+    // sinusoidal zones
+    SETUP(3, sinu, -d100, 0, -d100);
+    SETUP(4, sinu,   d30, 0,   d30);
+    SETUP(5, sinu, -d160, 0, -d160);
+    SETUP(6, sinu,  -d60, 0,  -d60);
+    SETUP(7, sinu,   d20, 0,   d20);
+    SETUP(8, sinu,  d140, 0,  d140);
 
-        // mollweide zones
-        SETUP(1, moll, -d100, 0, -d100);
+    // mollweide zones
+    SETUP(1, moll, -d100, 0, -d100);
 
-        // y0 ?
-        xy1 = P->pj[0]->fwd(lp, P->pj[0]); // zone 1
-        xy3 = P->pj[2]->fwd(lp, P->pj[2]); // zone 3
-        // y0 + xy1.y = xy3.y for lt = 40d44'11.8"
-        P->dy0 = xy3.y - xy1.y;
+    // y0 ?
+    xy1 = Q->pj[0]->fwd(lp, Q->pj[0]); // zone 1
+    xy3 = Q->pj[2]->fwd(lp, Q->pj[2]); // zone 3
+    // y0 + xy1.y = xy3.y for lt = 40d44'11.8"
+    Q->dy0 = xy3.y - xy1.y;
 
-        P->pj[0]->y0 = P->dy0;
+    Q->pj[0]->y0 = Q->dy0;
 
-        // mollweide zones (cont'd)
-        SETUP( 2, moll,   d30,  P->dy0,   d30);
-        SETUP( 9, moll, -d160, -P->dy0, -d160);
-        SETUP(10, moll,  -d60, -P->dy0,  -d60);
-        SETUP(11, moll,   d20, -P->dy0,   d20);
-        SETUP(12, moll,  d140, -P->dy0,  d140);
+    // mollweide zones (cont'd)
+    SETUP( 2, moll,   d30,  Q->dy0,   d30);
+    SETUP( 9, moll, -d160, -Q->dy0, -d160);
+    SETUP(10, moll,  -d60, -Q->dy0,  -d60);
+    SETUP(11, moll,   d20, -Q->dy0,   d20);
+    SETUP(12, moll,  d140, -Q->dy0,  d140);
 
-        P->inv = s_inverse;
-        P->fwd = s_forward;
-        P->es = 0.;
-ENDENTRY(P)
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    P->es = 0.;
+
+    return P;
+}
 
 
+#ifdef PJ_OMIT_SELFTEST
+int pj_igh_selftest (void) {return 0;}
+#else
+
+int pj_igh_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=igh   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223878.49745627123,  111701.07212763709},
+        { 223708.37131305804, -111701.07212763709},
+        {-222857.74059699223,  111701.07212763709},
+        {-223027.86674020503, -111701.07212763709},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.001790489447892545,   0.00089524655489191132},
+        { 0.0017904906685957927, -0.00089524655489191132},
+        {-0.001790496772112032,   0.00089524655489191132},
+        {-0.0017904955514087843, -0.00089524655489191132},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, inv_in, 0, s_inv_expect);
+}
+
 
+#endif
diff --git a/src/PJ_imw_p.c b/src/PJ_imw_p.c
index ae411116..10a06a44 100644
--- a/src/PJ_imw_p.c
+++ b/src/PJ_imw_p.c
@@ -1,151 +1,239 @@
-#define PROJ_PARMS__ \
-	double	P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
-	double	phi_1, phi_2, lam_1; \
-	double	*en; \
-	int	mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(imw_p, "International Map of the World Polyconic")
-	"\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
+    "\n\tMod. Polyconic, Ell\n\tlat_1= and lat_2= [lon_1=]";
+
 #define TOL 1e-10
 #define EPS 1e-10
-	static int
-phi12(PJ *P, double *del, double *sig) {
-	int err = 0;
-
-	if (!pj_param(P->ctx, P->params, "tlat_1").i ||
-		!pj_param(P->ctx, P->params, "tlat_2").i) {
-		err = -41;
-	} else {
-		P->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
-		P->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
-		*del = 0.5 * (P->phi_2 - P->phi_1);
-		*sig = 0.5 * (P->phi_2 + P->phi_1);
-		err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
-	}
-	return err;
+
+struct pj_opaque {
+    double  P, Pp, Q, Qp, R_1, R_2, sphi_1, sphi_2, C2; \
+    double  phi_1, phi_2, lam_1; \
+    double  *en; \
+    int mode; /* = 0, phi_1 and phi_2 != 0, = 1, phi_1 = 0, = -1 phi_2 = 0 */
+};
+
+
+static int phi12(PJ *P, double *del, double *sig) {
+    struct pj_opaque *Q = P->opaque;
+    int err = 0;
+
+    if (!pj_param(P->ctx, P->params, "tlat_1").i ||
+        !pj_param(P->ctx, P->params, "tlat_2").i) {
+        err = -41;
+    } else {
+        Q->phi_1 = pj_param(P->ctx, P->params, "rlat_1").f;
+        Q->phi_2 = pj_param(P->ctx, P->params, "rlat_2").f;
+        *del = 0.5 * (Q->phi_2 - Q->phi_1);
+        *sig = 0.5 * (Q->phi_2 + Q->phi_1);
+        err = (fabs(*del) < EPS || fabs(*sig) < EPS) ? -42 : 0;
+    }
+    return err;
 }
-	static XY
-loc_for(LP lp, PJ *P, double *yc) {
-	XY xy;
-
-	if (! lp.phi) {
-		xy.x = lp.lam;
-		xy.y = 0.;
-	} else {
-		double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
-
-		sp = sin(lp.phi);
-		m = pj_mlfn(lp.phi, sp, cos(lp.phi), P->en);
-		xa = P->Pp + P->Qp * m;
-		ya = P->P + P->Q * m;
-		R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
-		C = sqrt(R * R - xa * xa);
-		if (lp.phi < 0.) C = - C;
-		C += ya - R;
-		if (P->mode < 0) {
-			xb = lp.lam;
-			yb = P->C2;
-		} else {
-			t = lp.lam * P->sphi_2;
-			xb = P->R_2 * sin(t);
-			yb = P->C2 + P->R_2 * (1. - cos(t));
-		}
-		if (P->mode > 0) {
-			xc = lp.lam;
-			*yc = 0.;
-		} else {
-			t = lp.lam * P->sphi_1;
-			xc = P->R_1 * sin(t);
-			*yc = P->R_1 * (1. - cos(t));
-		}
-		D = (xb - xc)/(yb - *yc);
-		B = xc + D * (C + R - *yc);
-		xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
-		if (lp.phi > 0)
-			xy.x = - xy.x;
-		xy.x = (B + xy.x) / (1. + D * D);
-		xy.y = sqrt(R * R - xy.x * xy.x);
-		if (lp.phi > 0)
-			xy.y = - xy.y;
-		xy.y += C + R;
-	}
-	return (xy);
+
+
+static XY loc_for(LP lp, PJ *P, double *yc) {
+    struct pj_opaque *Q = P->opaque;
+    XY xy;
+
+    if (! lp.phi) {
+        xy.x = lp.lam;
+        xy.y = 0.;
+    } else {
+        double xa, ya, xb, yb, xc, D, B, m, sp, t, R, C;
+
+        sp = sin(lp.phi);
+        m = pj_mlfn(lp.phi, sp, cos(lp.phi), Q->en);
+        xa = Q->Pp + Q->Qp * m;
+        ya = Q->P + Q->Q * m;
+        R = 1. / (tan(lp.phi) * sqrt(1. - P->es * sp * sp));
+        C = sqrt(R * R - xa * xa);
+        if (lp.phi < 0.) C = - C;
+        C += ya - R;
+        if (Q->mode < 0) {
+            xb = lp.lam;
+            yb = Q->C2;
+        } else {
+            t = lp.lam * Q->sphi_2;
+            xb = Q->R_2 * sin(t);
+            yb = Q->C2 + Q->R_2 * (1. - cos(t));
+        }
+        if (Q->mode > 0) {
+            xc = lp.lam;
+            *yc = 0.;
+        } else {
+            t = lp.lam * Q->sphi_1;
+            xc = Q->R_1 * sin(t);
+            *yc = Q->R_1 * (1. - cos(t));
+        }
+        D = (xb - xc)/(yb - *yc);
+        B = xc + D * (C + R - *yc);
+        xy.x = D * sqrt(R * R * (1 + D * D) - B * B);
+        if (lp.phi > 0)
+            xy.x = - xy.x;
+        xy.x = (B + xy.x) / (1. + D * D);
+        xy.y = sqrt(R * R - xy.x * xy.x);
+        if (lp.phi > 0)
+            xy.y = - xy.y;
+        xy.y += C + R;
+    }
+    return xy;
 }
-FORWARD(e_forward); /* ellipsoid */
-	double yc;
-	xy = loc_for(lp, P, &yc);
-	return (xy);
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    double yc;
+
+    xy = loc_for(lp, P, &yc);
+    return (xy);
 }
-INVERSE(e_inverse); /* ellipsoid */
-	XY t;
-	double yc;
-
-	lp.phi = P->phi_2;
-	lp.lam = xy.x / cos(lp.phi);
-	do {
-		t = loc_for(lp, P, &yc);
-		lp.phi = ((lp.phi - P->phi_1) * (xy.y - yc) / (t.y - yc)) + P->phi_1;
-		lp.lam = lp.lam * xy.x / t.x;
-	} while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
-	return (lp);
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    XY t;
+    double yc;
+
+    lp.phi = Q->phi_2;
+    lp.lam = xy.x / cos(lp.phi);
+    do {
+        t = loc_for(lp, P, &yc);
+        lp.phi = ((lp.phi - Q->phi_1) * (xy.y - yc) / (t.y - yc)) + Q->phi_1;
+        lp.lam = lp.lam * xy.x / t.x;
+    } while (fabs(t.x - xy.x) > TOL || fabs(t.y - xy.y) > TOL);
+
+    return lp;
 }
-	static void
-xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
-	double F;
-
-	*sp = sin(phi);
-	*R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
-	F = P->lam_1 * *sp;
-	*y = *R * (1 - cos(F));
-	*x = *R * sin(F);
+
+
+static void xy(PJ *P, double phi, double *x, double *y, double *sp, double *R) {
+    double F;
+
+    *sp = sin(phi);
+    *R = 1./(tan(phi) * sqrt(1. - P->es * *sp * *sp ));
+    F = P->opaque->lam_1 * *sp;
+    *y = *R * (1 - cos(F));
+    *x = *R * sin(F);
 }
-FREEUP; if (P) { if (P->en) pj_dalloc(P->en); pj_dalloc(P); } }
-ENTRY1(imw_p, en)
-	double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
-	int i;
-
-	if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-	if( (i = phi12(P, &del, &sig)) != 0)
-		E_ERROR(i);
-	if (P->phi_2 < P->phi_1) { /* make sure P->phi_1 most southerly */
-		del = P->phi_1;
-		P->phi_1 = P->phi_2;
-		P->phi_2 = del;
-	}
-	if (pj_param(P->ctx, P->params, "tlon_1").i)
-		P->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
-	else { /* use predefined based upon latitude */
-		sig = fabs(sig * RAD_TO_DEG);
-		if (sig <= 60)		sig = 2.;
-		else if (sig <= 76) sig = 4.;
-		else				sig = 8.;
-		P->lam_1 = sig * DEG_TO_RAD;
-	}
-	P->mode = 0;
-	if (P->phi_1) xy(P, P->phi_1, &x1, &y1, &P->sphi_1, &P->R_1);
-	else {
-		P->mode = 1;
-		y1 = 0.;
-		x1 = P->lam_1;
-	}
-	if (P->phi_2) xy(P, P->phi_2, &x2, &T2, &P->sphi_2, &P->R_2);
-	else {
-		P->mode = -1;
-		T2 = 0.;
-		x2 = P->lam_1;
-	}
-	m1 = pj_mlfn(P->phi_1, P->sphi_1, cos(P->phi_1), P->en);
-	m2 = pj_mlfn(P->phi_2, P->sphi_2, cos(P->phi_2), P->en);
-	t = m2 - m1;
-	s = x2 - x1;
-	y2 = sqrt(t * t - s * s) + y1;
-	P->C2 = y2 - T2;
-	t = 1. / t;
-	P->P = (m2 * y1 - m1 * y2) * t;
-	P->Q = (y2 - y1) * t;
-	P->Pp = (m2 * x1 - m1 * x2) * t;
-	P->Qp = (x2 - x1) * t;
-	P->fwd = e_forward;
-	P->inv = e_inverse;
-ENDENTRY(P)
+
+
+static void *freeup_new (PJ *P) {              /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(imw_p) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    double del, sig, s, t, x1, x2, T2, y1, m1, m2, y2;
+    int i;
+
+    if (!(Q->en = pj_enfn(P->es))) E_ERROR_0;
+    if( (i = phi12(P, &del, &sig)) != 0)
+        E_ERROR(i);
+    if (Q->phi_2 < Q->phi_1) { /* make sure P->phi_1 most southerly */
+        del = Q->phi_1;
+        Q->phi_1 = Q->phi_2;
+        Q->phi_2 = del;
+    }
+    if (pj_param(P->ctx, P->params, "tlon_1").i)
+        Q->lam_1 = pj_param(P->ctx, P->params, "rlon_1").f;
+    else { /* use predefined based upon latitude */
+        sig = fabs(sig * RAD_TO_DEG);
+        if (sig <= 60)      sig = 2.;
+        else if (sig <= 76) sig = 4.;
+        else                sig = 8.;
+        Q->lam_1 = sig * DEG_TO_RAD;
+    }
+    Q->mode = 0;
+    if (Q->phi_1) xy(P, Q->phi_1, &x1, &y1, &Q->sphi_1, &Q->R_1);
+    else {
+        Q->mode = 1;
+        y1 = 0.;
+        x1 = Q->lam_1;
+    }
+    if (Q->phi_2) xy(P, Q->phi_2, &x2, &T2, &Q->sphi_2, &Q->R_2);
+    else {
+        Q->mode = -1;
+        T2 = 0.;
+        x2 = Q->lam_1;
+    }
+    m1 = pj_mlfn(Q->phi_1, Q->sphi_1, cos(Q->phi_1), Q->en);
+    m2 = pj_mlfn(Q->phi_2, Q->sphi_2, cos(Q->phi_2), Q->en);
+    t = m2 - m1;
+    s = x2 - x1;
+    y2 = sqrt(t * t - s * s) + y1;
+    Q->C2 = y2 - T2;
+    t = 1. / t;
+    Q->P = (m2 * y1 - m1 * y2) * t;
+    Q->Q = (y2 - y1) * t;
+    Q->Pp = (m2 * x1 - m1 * x2) * t;
+    Q->Qp = (x2 - x1) * t;
+
+    P->fwd = e_forward;
+    P->inv = e_inverse;
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_imw_p_selftest (void) {return 0;}
+#else
+
+int pj_imw_p_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=imw_p   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222588.4411393762,   55321.128653809537},
+        { 222756.90637768712, -165827.58428832365},
+        {-222588.4411393762,   55321.128653809537},
+        {-222756.90637768712, -165827.58428832365},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017966991379592214, 0.50090492361427374},
+        { 0.0017966979081574697, 0.49909507588689922},
+        {-0.0017966991379592214, 0.50090492361427374},
+        {-0.0017966979081574697, 0.49909507588689922},
+    };
+
+    return pj_generic_selftest (e_args, 0, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, 0, inv_in, e_inv_expect, 0);
+}
+
+
+#endif
diff --git a/src/PJ_isea.c b/src/PJ_isea.c
index 7c09189c..f1bc89a7 100644
--- a/src/PJ_isea.c
+++ b/src/PJ_isea.c
@@ -1,1119 +1,1176 @@
-/*
- * This code was entirely written by Nathan Wagner
- * and is in the public domain.
- */
-
-#include <math.h>
-#include <stdio.h>
-#include <stdlib.h>
-#include <float.h>
-
-#ifndef M_PI
-#  define M_PI 3.14159265358979323846
-#endif
-
-/*
- * Proj 4 provides its own entry points into
- * the code, so none of the library functions
- * need to be global
- */
-#define ISEA_STATIC static
-#ifndef ISEA_STATIC
-#define ISEA_STATIC
-#endif
-
-struct hex {
-        int iso;
-        int x, y, z;
-};
-
-/* y *must* be positive down as the xy /iso conversion assumes this */
-ISEA_STATIC
-int hex_xy(struct hex *h) {
-	if (!h->iso) return 1;
-	if (h->x >= 0) {
-		h->y = -h->y - (h->x+1)/2;
-	} else {
-		/* need to round toward -inf, not toward zero, so x-1 */
-		h->y = -h->y - h->x/2;
-	}
-	h->iso = 0;
-
-	return 1;
-}
-
-ISEA_STATIC
-int hex_iso(struct hex *h) {
-	if (h->iso) return 1;
-
-	if (h->x >= 0) {
-		h->y = (-h->y - (h->x+1)/2);
-	} else {
-		/* need to round toward -inf, not toward zero, so x-1 */
-		h->y = (-h->y - (h->x)/2);
-	}
-
-	h->z = -h->x - h->y;
-	h->iso = 1;
-	return 1;
-}
-
-ISEA_STATIC
-int hexbin2(double width, double x, double y,
-                int *i, int *j) {
-	double z, rx, ry, rz;
-	double abs_dx, abs_dy, abs_dz;
-	int ix, iy, iz, s;
-	struct hex h;
-
-	x = x / cos(30 * M_PI / 180.0); /* rotated X coord */
-	y = y - x / 2.0; /* adjustment for rotated X */
-
-	/* adjust for actual hexwidth */
-	x /= width;
-	y /= width;
-
-	z = -x - y;
-
-	ix = rx = floor(x + 0.5);
-	iy = ry = floor(y + 0.5);
-	iz = rz = floor(z + 0.5);
-
-	s = ix + iy + iz;
-
-	if (s) {
-		abs_dx = fabs(rx - x);
-		abs_dy = fabs(ry - y);
-		abs_dz = fabs(rz - z);
-
-		if (abs_dx >= abs_dy && abs_dx >= abs_dz) {
-			ix -= s;
-		} else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {
-			iy -= s;
-		} else {
-			iz -= s;
-		}
-	}
-	h.x = ix;
-	h.y = iy;
-	h.z = iz;
-	h.iso = 1;
-
-	hex_xy(&h);
-	*i = h.x;
-	*j = h.y;
-        return ix * 100 + iy;
-}
-#ifndef ISEA_STATIC
-#define ISEA_STATIC
-#endif
-
-enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };
-enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };
-enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,
-	ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX
-};
-
-struct isea_dgg {
-	int	polyhedron; /* ignored, icosahedron */
-	double	o_lat, o_lon, o_az; /* orientation, radians */
-	int	pole; /* true if standard snyder */
-	int	topology; /* ignored, hexagon */
-	int	aperture; /* valid values depend on partitioning method */
-	int	resolution;
-	double	radius; /* radius of the earth in meters, ignored 1.0 */
-	int	output; /* an isea_address_form */
-	int	triangle; /* triangle of last transformed point */
-	int	quad; /* quad of last transformed point */
-	unsigned long serial;
-};
-
-struct isea_pt {
-	double x, y;
-};
-
-struct isea_geo {
-	double lon, lat;
-};
-
-struct isea_address {
-	int	type; /* enum isea_address_form */
-	int	number;
-	double	x,y; /* or i,j or lon,lat depending on type */
-};
-
-/* ENDINC */
-
-enum snyder_polyhedron {
-	SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,
-	SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,
-	SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,
-	SNYDER_POLY_ICOSAHEDRON
-};
-
-struct snyder_constants {
-	double          g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;
-};
-
-/* TODO put these in radians to avoid a later conversion */
-ISEA_STATIC
-struct snyder_constants constants[] = {
-	{23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},
-	{20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},
-	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-	{0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},
-	{37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},
-};
-
-#define E 52.62263186
-#define F 10.81231696
-
-#define DEG60 1.04719755119659774614
-#define DEG120 2.09439510239319549229
-#define DEG72 1.25663706143591729537
-#define DEG90 1.57079632679489661922
-#define DEG144 2.51327412287183459075
-#define DEG36 0.62831853071795864768
-#define DEG108 1.88495559215387594306
-#define DEG180 M_PI
-/* sqrt(5)/M_PI */
-#define ISEA_SCALE 0.8301572857837594396028083
-
-/* 26.565051177 degrees */
-#define V_LAT 0.46364760899944494524
-
-#define RAD2DEG (180.0/M_PI)
-#define DEG2RAD (M_PI/180.0)
-
-ISEA_STATIC
-struct isea_geo vertex[] = {
-	{0.0, DEG90},
-	{DEG180, V_LAT},
-	{-DEG108, V_LAT},
-	{-DEG36, V_LAT},
-	{DEG36, V_LAT},
-	{DEG108, V_LAT},
-	{-DEG144, -V_LAT},
-	{-DEG72, -V_LAT},
-	{0.0, -V_LAT},
-	{DEG72, -V_LAT},
-	{DEG144, -V_LAT},
-	{0.0, -DEG90}
-};
-
-/* TODO make an isea_pt array of the vertices as well */
-
-static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};
-
-/* 52.62263186 */
-#define E_RAD 0.91843818702186776133
-
-/* 10.81231696 */
-#define F_RAD 0.18871053072122403508
-
-/* triangle Centers */
-struct isea_geo icostriangles[] = {
-	{0.0, 0.0},
-	{-DEG144, E_RAD},
-	{-DEG72, E_RAD},
-	{0.0, E_RAD},
-	{DEG72, E_RAD},
-	{DEG144, E_RAD},
-	{-DEG144, F_RAD},
-	{-DEG72, F_RAD},
-	{0.0, F_RAD},
-	{DEG72, F_RAD},
-	{DEG144, F_RAD},
-	{-DEG108, -F_RAD},
-	{-DEG36, -F_RAD},
-	{DEG36, -F_RAD},
-	{DEG108, -F_RAD},
-	{DEG180, -F_RAD},
-	{-DEG108, -E_RAD},
-	{-DEG36, -E_RAD},
-	{DEG36, -E_RAD},
-	{DEG108, -E_RAD},
-	{DEG180, -E_RAD},
-};
-
-static double
-az_adjustment(int triangle)
-{
-	double          adj;
-
-	struct isea_geo v;
-	struct isea_geo c;
-
-	v = vertex[tri_v1[triangle]];
-	c = icostriangles[triangle];
-
-	/* TODO looks like the adjustment is always either 0 or 180 */
-	/* at least if you pick your vertex carefully */
-	adj = atan2(cos(v.lat) * sin(v.lon - c.lon),
-		    cos(c.lat) * sin(v.lat)
-		    - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));
-	return adj;
-}
-
-/* R tan(g) sin(60) */
-#define TABLE_G 0.6615845383
-
-/* H = 0.25 R tan g = */
-#define TABLE_H 0.1909830056
-
-#define RPRIME 0.91038328153090290025
-
-ISEA_STATIC
-struct isea_pt
-isea_triangle_xy(int triangle)
-{
-	struct isea_pt  c;
-	double Rprime = 0.91038328153090290025;
-
-	triangle = (triangle - 1) % 20;
-
-	c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;
-	if (triangle > 9) {
-		c.x += TABLE_G;
-	}
-	switch (triangle / 5) {
-	case 0:
-		c.y = 5.0 * TABLE_H;
-		break;
-	case 1:
-		c.y = TABLE_H;
-		break;
-	case 2:
-		c.y = -TABLE_H;
-		break;
-	case 3:
-		c.y = -5.0 * TABLE_H;
-		break;
-	default:
-		/* should be impossible */
-		exit(EXIT_FAILURE);
-	};
-	c.x *= Rprime;
-	c.y *= Rprime;
-
-	return c;
-}
-
-/* snyder eq 14 */
-static double
-sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)
-{
-	double          az;
-
-	az = atan2(cos(t_lat) * sin(t_lon - f_lon),
-		   cos(f_lat) * sin(t_lat)
-		   - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)
-		);
-	return az;
-}
-
-/* coord needs to be in radians */
-ISEA_STATIC
-int
-isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)
-{
-	int             i;
-
-	/*
-	 * spherical distance from center of polygon face to any of its
-	 * vertexes on the globe
-	 */
-	double          g;
-
-	/*
-	 * spherical angle between radius vector to center and adjacent edge
-	 * of spherical polygon on the globe
-	 */
-	double          G;
-
-	/*
-	 * plane angle between radius vector to center and adjacent edge of
-	 * plane polygon
-	 */
-	double          theta;
-
-	/* additional variables from snyder */
-	double          q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,
-	                x, y;
-
-	/* variables used to store intermediate results */
-	double          cot_theta, tan_g, az_offset;
-
-	/* how many multiples of 60 degrees we adjust the azimuth */
-	int             Az_adjust_multiples;
-
-	struct snyder_constants c;
-
-	/*
-	 * TODO by locality of reference, start by trying the same triangle
-	 * as last time
-	 */
-
-	/* TODO put these constants in as radians to begin with */
-	c = constants[SNYDER_POLY_ICOSAHEDRON];
-	theta = c.theta * DEG2RAD;
-	g = c.g * DEG2RAD;
-	G = c.G * DEG2RAD;
-
-	for (i = 1; i <= 20; i++) {
-		double          z;
-		struct isea_geo center;
-
-		center = icostriangles[i];
-
-		/* step 1 */
-		z = acos(sin(center.lat) * sin(ll->lat)
-			 + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));
-		/* not on this triangle */
-		if (z > g + 0.000005) { /* TODO DBL_EPSILON */
-			continue;
-		}
-
-		Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);
-
-		/* step 2 */
-
-		/* This calculates "some" vertex coordinate */
-		az_offset = az_adjustment(i);
-
-		Az -= az_offset;
-
-		/* TODO I don't know why we do this.  It's not in snyder */
-		/* maybe because we should have picked a better vertex */
-		if (Az < 0.0) {
-			Az += 2.0 * M_PI;
-		}
-		/*
-		 * adjust Az for the point to fall within the range of 0 to
-		 * 2(90 - theta) or 60 degrees for the hexagon, by
-		 * and therefore 120 degrees for the triangle
-		 * of the icosahedron
-		 * subtracting or adding multiples of 60 degrees to Az and
-		 * recording the amount of adjustment
-		 */
-
-		Az_adjust_multiples = 0;
-		while (Az < 0.0) {
-			Az += DEG120;
-			Az_adjust_multiples--;
-		}
-		while (Az > DEG120 + DBL_EPSILON) {
-			Az -= DEG120;
-			Az_adjust_multiples++;
-		}
-
-		/* step 3 */
-		cot_theta = 1.0 / tan(theta);
-		tan_g = tan(g);	/* TODO this is a constant */
-
-		/* Calculate q from eq 9. */
-		/* TODO cot_theta is cot(30) */
-		q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);
-
-		/* not in this triangle */
-		if (z > q + 0.000005) {
-			continue;
-		}
-		/* step 4 */
-
-		/* Apply equations 5-8 and 10-12 in order */
-
-		/* eq 5 */
-		/* Rprime = 0.9449322893 * R; */
-		/* R' in the paper is for the truncated */
-		Rprime = 0.91038328153090290025;
-
-		/* eq 6 */
-		H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));
-
-		/* eq 7 */
-		/* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */
-		Ag = Az + G + H - DEG180;
-
-		/* eq 8 */
-		Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);
-
-		/* eq 10 */
-		/* cot(theta) = 1.73205080756887729355 */
-		dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);
-
-		/* eq 11 */
-		f = dprime / (2.0 * Rprime * sin(q / 2.0));
-
-		/* eq 12 */
-		rho = 2.0 * Rprime * f * sin(z / 2.0);
-
-		/*
-		 * add back the same 60 degree multiple adjustment from step
-		 * 2 to Azprime
-		 */
-
-		Azprime += DEG120 * Az_adjust_multiples;
-
-		/* calculate rectangular coordinates */
-
-		x = rho * sin(Azprime);
-		y = rho * cos(Azprime);
-
-		/*
-		 * TODO
-		 * translate coordinates to the origin for the particular
-		 * hexagon on the flattened polyhedral map plot
-		 */
-
-		out->x = x;
-		out->y = y;
-
-		return i;
-	}
-
-	/*
-	 * should be impossible, this implies that the coordinate is not on
-	 * any triangle
-	 */
-
-	fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",
-		ll->lon * RAD2DEG, ll->lat * RAD2DEG);
-
-	exit(EXIT_FAILURE);
-
-	/* not reached */
-	return 0;		/* supresses a warning */
-}
-
-/*
- * return the new coordinates of any point in orginal coordinate system.
- * Define a point (newNPold) in orginal coordinate system as the North Pole in
- * new coordinate system, and the great circle connect the original and new
- * North Pole as the lon0 longitude in new coordinate system, given any point
- * in orginal coordinate system, this function return the new coordinates.
- */
-
-#define PRECISION 0.0000000000005
-
-/* formula from Snyder, Map Projections: A working manual, p31 */
-/*
- * old north pole at np in new coordinates
- * could be simplified a bit with fewer intermediates
- *
- * TODO take a result pointer
- */
-ISEA_STATIC
-struct isea_geo
-snyder_ctran(struct isea_geo * np, struct isea_geo * pt)
-{
-	struct isea_geo npt;
-	double          alpha, phi, lambda, lambda0, beta, lambdap, phip;
-	double          sin_phip;
-	double          lp_b;	/* lambda prime minus beta */
-	double          cos_p, sin_a;
-
-	phi = pt->lat;
-	lambda = pt->lon;
-	alpha = np->lat;
-	beta = np->lon;
-	lambda0 = beta;
-
-	cos_p = cos(phi);
-	sin_a = sin(alpha);
-
-	/* mpawm 5-7 */
-	sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);
-
-	/* mpawm 5-8b */
-
-	/* use the two argument form so we end up in the right quadrant */
-	lp_b = atan2(cos_p * sin(lambda - lambda0),
-	   (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));
-
-	lambdap = lp_b + beta;
-
-	/* normalize longitude */
-	/* TODO can we just do a modulus ? */
-	lambdap = fmod(lambdap, 2 * M_PI);
-	while (lambdap > M_PI)
-		lambdap -= 2 * M_PI;
-	while (lambdap < -M_PI)
-		lambdap += 2 * M_PI;
-
-	phip = asin(sin_phip);
-
-	npt.lat = phip;
-	npt.lon = lambdap;
-
-	return npt;
-}
-
-ISEA_STATIC
-struct isea_geo
-isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)
-{
-	struct isea_geo npt;
-
-	np->lon += M_PI;
-	npt = snyder_ctran(np, pt);
-	np->lon -= M_PI;
-
-	npt.lon -= (M_PI - lon0 + np->lon);
-
-	/*
-	 * snyder is down tri 3, isea is along side of tri1 from vertex 0 to
-	 * vertex 1 these are 180 degrees apart
-	 */
-	npt.lon += M_PI;
-	/* normalize longitude */
-	npt.lon = fmod(npt.lon, 2 * M_PI);
-	while (npt.lon > M_PI)
-		npt.lon -= 2 * M_PI;
-	while (npt.lon < -M_PI)
-		npt.lon += 2 * M_PI;
-
-	return npt;
-}
-
-/* in radians */
-#define ISEA_STD_LAT 1.01722196792335072101
-#define ISEA_STD_LON .19634954084936207740
-
-/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */
-
-ISEA_STATIC
-int
-isea_grid_init(struct isea_dgg * g)
-{
-	if (!g)
-		return 0;
-
-	g->polyhedron = 20;
-	g->o_lat = ISEA_STD_LAT;
-	g->o_lon = ISEA_STD_LON;
-	g->o_az = 0.0;
-	g->aperture = 4;
-	g->resolution = 6;
-	g->radius = 1.0;
-	g->topology = 6;
-
-	return 1;
-}
-
-ISEA_STATIC
-int
-isea_orient_isea(struct isea_dgg * g)
-{
-	if (!g)
-		return 0;
-	g->o_lat = ISEA_STD_LAT;
-	g->o_lon = ISEA_STD_LON;
-	g->o_az = 0.0;
-	return 1;
-}
-
-ISEA_STATIC
-int
-isea_orient_pole(struct isea_dgg * g)
-{
-	if (!g)
-		return 0;
-	g->o_lat = M_PI / 2.0;
-	g->o_lon = 0.0;
-	g->o_az = 0;
-	return 1;
-}
-
-ISEA_STATIC
-int
-isea_transform(struct isea_dgg * g, struct isea_geo * in,
-	       struct isea_pt * out)
-{
-	struct isea_geo i, pole;
-	int             tri;
-
-	pole.lat = g->o_lat;
-	pole.lon = g->o_lon;
-
-	i = isea_ctran(&pole, in, g->o_az);
-
-	tri = isea_snyder_forward(&i, out);
-	out->x *= g->radius;
-	out->y *= g->radius;
-	g->triangle = tri;
-
-	return tri;
-}
-
-#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)
-
-ISEA_STATIC
-void
-isea_rotate(struct isea_pt * pt, double degrees)
-{
-	double          rad;
-
-	double          x, y;
-
-	rad = -degrees * M_PI / 180.0;
-	while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;
-	while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;
-
-	x = pt->x * cos(rad) + pt->y * sin(rad);
-	y = -pt->x * sin(rad) + pt->y * cos(rad);
-
-	pt->x = x;
-	pt->y = y;
-}
-
-ISEA_STATIC
-int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {
-	struct isea_pt tc; /* center of triangle */
-
-	if (DOWNTRI(tri)) {
-		isea_rotate(pt, 180.0);
-	}
-	tc = isea_triangle_xy(tri);
-	tc.x *= radius;
-	tc.y *= radius;
-	pt->x += tc.x;
-	pt->y += tc.y;
-
-	return tri;
-}
-
-/* convert projected triangle coords to quad xy coords, return quad number */
-ISEA_STATIC
-int
-isea_ptdd(int tri, struct isea_pt *pt) {
-	int             downtri, quad;
-
-	downtri = (((tri - 1) / 5) % 2 == 1);
-	quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;
-
-	isea_rotate(pt, downtri ? 240.0 : 60.0);
-	if (downtri) {
-		pt->x += 0.5;
-		/* pt->y += cos(30.0 * M_PI / 180.0); */
-		pt->y += .86602540378443864672;
-	}
-	return quad;
-}
-
-ISEA_STATIC
-int
-isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)
-{
-	struct isea_pt  v;
-	double          hexwidth;
-	double          sidelength;	/* in hexes */
-	int             d, i;
-	int             maxcoord;
-	struct hex      h;
-
-	/* This is the number of hexes from apex to base of a triangle */
-	sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;
-
-	/* apex to base is cos(30deg) */
-	hexwidth = cos(M_PI / 6.0) / sidelength;
-
-	/* TODO I think sidelength is always x.5, so
-	 * (int)sidelength * 2 + 1 might be just as good
-	 */
-	maxcoord = (int) (sidelength * 2.0 + 0.5);
-
-	v = *pt;
-	hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
-	h.iso = 0;
-	hex_iso(&h);
-
-	d = h.x - h.z;
-	i = h.x + h.y + h.y;
-	
-	/*
-	 * you want to test for max coords for the next quad in the same
-	 * "row" first to get the case where both are max
-	 */
-	if (quad <= 5) {
-		if (d == 0 && i == maxcoord) {
-			/* north pole */
-			quad = 0;
-			d = 0;
-			i = 0;
-		} else if (i == maxcoord) {
-			/* upper right in next quad */
-			quad += 1;
-			if (quad == 6)
-				quad = 1;
-			i = maxcoord - d;
-			d = 0;
-		} else if (d == maxcoord) {
-			/* lower right in quad to lower right */
-			quad += 5;
-			d = 0;
-		}
-	} else if (quad >= 6) {
-		if (i == 0 && d == maxcoord) {
-			/* south pole */
-			quad = 11;
-			d = 0;
-			i = 0;
-		} else if (d == maxcoord) {
-			/* lower right in next quad */
-			quad += 1;
-			if (quad == 11)
-				quad = 6;
-			d = maxcoord - i;
-			i = 0;
-		} else if (i == maxcoord) {
-			/* upper right in quad to upper right */
-			quad = (quad - 4) % 5;
-			i = 0;
-		}
-	}
-
-	di->x = d;
-	di->y = i;
-
-	g->quad = quad;
-	return quad;
-}
-
-ISEA_STATIC
-int
-isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {
-	struct isea_pt  v;
-	double          hexwidth;
-	int             sidelength;	/* in hexes */
-	struct hex      h;
-
-	if (g->aperture == 3 && g->resolution % 2 != 0) {
-		return isea_dddi_ap3odd(g, quad, pt, di);
-	}
-	/* todo might want to do this as an iterated loop */
-	if (g->aperture >0) {
-		sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
-	} else {
-		sidelength = g->resolution;
-	}
-
-	hexwidth = 1.0 / sidelength;
-
-	v = *pt;
-	isea_rotate(&v, -30.0);
-	hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);
-	h.iso = 0;
-	hex_iso(&h);
-
-	/* we may actually be on another quad */
-	if (quad <= 5) {
-		if (h.x == 0 && h.z == -sidelength) {
-			/* north pole */
-			quad = 0;
-			h.z = 0;
-			h.y = 0;
-			h.x = 0;
-		} else if (h.z == -sidelength) {
-			quad = quad + 1;
-			if (quad == 6)
-				quad = 1;
-			h.y = sidelength - h.x;
-			h.z = h.x - sidelength;
-			h.x = 0;
-		} else if (h.x == sidelength) {
-			quad += 5;
-			h.y = -h.z;
-			h.x = 0;
-		}
-	} else if (quad >= 6) {
-		if (h.z == 0 && h.x == sidelength) {
-			/* south pole */
-			quad = 11;
-			h.x = 0;
-			h.y = 0;
-			h.z = 0;
-		} else if (h.x == sidelength) {
-			quad = quad + 1;
-			if (quad == 11)
-				quad = 6;
-			h.x = h.y + sidelength;
-			h.y = 0;
-			h.z = -h.x;
-		} else if (h.y == -sidelength) {
-			quad -= 4;
-			h.y = 0;
-			h.z = -h.x;
-		}
-	}
-	di->x = h.x;
-	di->y = -h.z;
-
-	g->quad = quad;
-	return quad;
-}
-
-ISEA_STATIC
-int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,
-	       	struct isea_pt *di) {
-	struct isea_pt  v;
-	int             quad;
-
-	v = *pt;
-	quad = isea_ptdd(tri, &v);
-	quad = isea_dddi(g, quad, &v, di);
-	return quad;
-}
-
-/* q2di to seqnum */
-ISEA_STATIC
-int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {
-	int             sidelength;
-	int             sn, height;
-	int             hexes;
-
-	if (quad == 0) {
-		g->serial = 1;
-		return g->serial;
-	}
-	/* hexes in a quad */
-	hexes = (int) (pow(g->aperture, g->resolution) + 0.5);
-	if (quad == 11) {
-		g->serial = 1 + 10 * hexes + 1;
-		return g->serial;
-	}
-	if (g->aperture == 3 && g->resolution % 2 == 1) {
-		height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));
-		sn = ((int) di->x) * height;
-		sn += ((int) di->y) / height;
-		sn += (quad - 1) * hexes;
-		sn += 2;
-	} else {
-		sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
-		sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;
-	}
-
-	g->serial = sn;
-	return sn;
-}
-
-/* TODO just encode the quad in the d or i coordinate
- * quad is 0-11, which can be four bits.
- * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf
- */
-/* convert a q2di to global hex coord */
-ISEA_STATIC
-int isea_hex(struct isea_dgg *g, int tri,
-		struct isea_pt *pt, struct isea_pt *hex) {
-	struct isea_pt v;
-	int sidelength;
-	int d, i, x, y, quad;
-
-	quad = isea_ptdi(g, tri, pt, &v);
-
-	hex->x = ((int)v.x << 4) + quad;
-	hex->y = v.y;
-
-	return 1;
-
-	d = v.x;
-	i = v.y;
-
-	/* Aperture 3 odd resolutions */
-	if (g->aperture == 3 && g->resolution % 2 != 0) {
-		int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);
-
-		d += offset * ((g->quad-1) % 5);
-		i += offset * ((g->quad-1) % 5);
-
-		if (quad == 0) {
-			d = 0;
-			i = offset;
-		} else if (quad == 11) {
-			d = 2 * offset;
-			i = 0;
-		} else if (quad > 5) {
-			d += offset;
-		}
-
-		x = (2*d - i) /3;
-		y = (2*i - d) /3;
-
-		hex->x = x + offset / 3;
-		hex->y = y + 2 * offset / 3;
-		return 1;
-	}
-
-	/* aperture 3 even resolutions and aperture 4 */
-	sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);
-	if (g->quad == 0) {
-		hex->x = 0;
-		hex->y = sidelength;
-	} else if (g->quad == 11) {
-		hex->x = sidelength * 2;
-		hex->y = 0;
-	} else {
-		hex->x = d + sidelength * ((g->quad-1) % 5);
-		if (g->quad > 5) hex->x += sidelength;
-		hex->y = i + sidelength * ((g->quad-1) % 5);
-	}
-
-	return 1;
-}
-
-ISEA_STATIC
-struct isea_pt
-isea_forward(struct isea_dgg *g, struct isea_geo *in)
-{
-	int             tri;
-	struct isea_pt  out, coord;
-
-	tri = isea_transform(g, in, &out);
-
-	if (g->output == ISEA_PLANE) {
-		isea_tri_plane(tri, &out, g->radius);
-		return out;
-	}
-
-	/* convert to isea standard triangle size */
-	out.x = out.x / g->radius * ISEA_SCALE;
-	out.y = out.y / g->radius * ISEA_SCALE;
-	out.x += 0.5;
-	out.y += 2.0 * .14433756729740644112;
-
-	switch (g->output) {
-	case ISEA_PROJTRI:
-		/* nothing to do, already in projected triangle */
-		break;
-	case ISEA_VERTEX2DD:
-		g->quad = isea_ptdd(tri, &out);
-		break;
-	case ISEA_Q2DD:
-		/* Same as above, we just don't print as much */
-		g->quad = isea_ptdd(tri, &out);
-		break;
-	case ISEA_Q2DI:
-		g->quad = isea_ptdi(g, tri, &out, &coord);
-		return coord;
-		break;
-	case ISEA_SEQNUM:
-		isea_ptdi(g, tri, &out, &coord);
-		/* disn will set g->serial */
-		isea_disn(g, g->quad, &coord);
-		return coord;
-		break;
-	case ISEA_HEX:
-		isea_hex(g, tri, &out, &coord);
-		return coord;
-		break;
-	}
-
-	return out;
-}
-/*
- * Proj 4 integration code follows
- */
-
-#define PROJ_PARMS__ \
-	struct isea_dgg dgg;
-
-#define PJ_LIB__
-#include <projects.h>
-
-PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";
-
-FORWARD(s_forward);
-	struct isea_pt out;
-	struct isea_geo in;
-
-	in.lon = lp.lam;
-	in.lat = lp.phi;
-
-	out = isea_forward(&P->dgg, &in);
-	
-	xy.x = out.x;
-	xy.y = out.y;
-
-	return xy;
-}
-FREEUP; if (P) pj_dalloc(P); }
-
-ENTRY0(isea)
-	char *opt;
-
-        P->fwd = s_forward;
-        isea_grid_init(&P->dgg);
-
-        P->dgg.output = ISEA_PLANE;
-/*		P->dgg.radius = P->a; / * otherwise defaults to 1 */
-	/* calling library will scale, I think */
-
-	opt = pj_param(P->ctx,P->params, "sorient").s;
-	if (opt) {
-		if (!strcmp(opt, "isea")) {
-			isea_orient_isea(&P->dgg);
-		} else if (!strcmp(opt, "pole")) {
-			isea_orient_pole(&P->dgg);
-		} else {
-			E_ERROR(-34);
-		}
-	}
-
-	if (pj_param(P->ctx,P->params, "tazi").i) {
-		P->dgg.o_az = pj_param(P->ctx,P->params, "razi").f;
-	}
-
-	if (pj_param(P->ctx,P->params, "tlon_0").i) {
-		P->dgg.o_lon = pj_param(P->ctx,P->params, "rlon_0").f;
-	}
-
-	if (pj_param(P->ctx,P->params, "tlat_0").i) {
-		P->dgg.o_lat = pj_param(P->ctx,P->params, "rlat_0").f;
-	}
-
-	if (pj_param(P->ctx,P->params, "taperture").i) {
-		P->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
-	}
-
-	if (pj_param(P->ctx,P->params, "tresolution").i) {
-		P->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
-	}
-
-	opt = pj_param(P->ctx,P->params, "smode").s;
-	if (opt) {
-		if (!strcmp(opt, "plane")) {
-			P->dgg.output = ISEA_PLANE;
-		} else if (!strcmp(opt, "di")) {
-			P->dgg.output = ISEA_Q2DI;
-		}
-		else if (!strcmp(opt, "dd")) {
-			P->dgg.output = ISEA_Q2DD;
-		}
-		else if (!strcmp(opt, "hex")) {
-			P->dgg.output = ISEA_HEX;
-		}
-		else {
-			/* TODO verify error code.  Possibly eliminate magic */
-			E_ERROR(-34);
-		}
-	}
-
-	if (pj_param(P->ctx,P->params, "trescale").i) {
-		P->dgg.radius = ISEA_SCALE;
-	}
-
-	if (pj_param(P->ctx,P->params, "tresolution").i) {
-		P->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;
-	} else {
-		P->dgg.resolution = 4;
-	}
-
-	if (pj_param(P->ctx,P->params, "taperture").i) {
-		P->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;
-	} else {
-		P->dgg.aperture = 3;
-	}
-
-ENDENTRY(P)
+/*

+ * This code was entirely written by Nathan Wagner

+ * and is in the public domain.

+ */

+

+#include <math.h>

+#include <stdio.h>

+#include <stdlib.h>

+#include <float.h>

+

+#ifndef M_PI

+#  define M_PI 3.14159265358979323846

+#endif

+

+/*

+ * Proj 4 provides its own entry points into

+ * the code, so none of the library functions

+ * need to be global

+ */

+#define ISEA_STATIC static

+#ifndef ISEA_STATIC

+#define ISEA_STATIC

+#endif

+

+struct hex {

+        int iso;

+        int x, y, z;

+};

+

+/* y *must* be positive down as the xy /iso conversion assumes this */

+ISEA_STATIC

+int hex_xy(struct hex *h) {

+    if (!h->iso) return 1;

+    if (h->x >= 0) {

+        h->y = -h->y - (h->x+1)/2;

+    } else {

+        /* need to round toward -inf, not toward zero, so x-1 */

+        h->y = -h->y - h->x/2;

+    }

+    h->iso = 0;

+

+    return 1;

+}

+

+ISEA_STATIC

+int hex_iso(struct hex *h) {

+    if (h->iso) return 1;

+

+    if (h->x >= 0) {

+        h->y = (-h->y - (h->x+1)/2);

+    } else {

+        /* need to round toward -inf, not toward zero, so x-1 */

+        h->y = (-h->y - (h->x)/2);

+    }

+

+    h->z = -h->x - h->y;

+    h->iso = 1;

+    return 1;

+}

+

+ISEA_STATIC

+int hexbin2(double width, double x, double y,

+                int *i, int *j) {

+    double z, rx, ry, rz;

+    double abs_dx, abs_dy, abs_dz;

+    int ix, iy, iz, s;

+    struct hex h;

+

+    x = x / cos(30 * M_PI / 180.0); /* rotated X coord */

+    y = y - x / 2.0; /* adjustment for rotated X */

+

+    /* adjust for actual hexwidth */

+    x /= width;

+    y /= width;

+

+    z = -x - y;

+

+    ix = rx = floor(x + 0.5);

+    iy = ry = floor(y + 0.5);

+    iz = rz = floor(z + 0.5);

+

+    s = ix + iy + iz;

+

+    if (s) {

+        abs_dx = fabs(rx - x);

+        abs_dy = fabs(ry - y);

+        abs_dz = fabs(rz - z);

+

+        if (abs_dx >= abs_dy && abs_dx >= abs_dz) {

+            ix -= s;

+        } else if (abs_dy >= abs_dx && abs_dy >= abs_dz) {

+            iy -= s;

+        } else {

+            iz -= s;

+        }

+    }

+    h.x = ix;

+    h.y = iy;

+    h.z = iz;

+    h.iso = 1;

+

+    hex_xy(&h);

+    *i = h.x;

+    *j = h.y;

+        return ix * 100 + iy;

+}

+#ifndef ISEA_STATIC

+#define ISEA_STATIC

+#endif

+

+enum isea_poly { ISEA_NONE, ISEA_ICOSAHEDRON = 20 };

+enum isea_topology { ISEA_HEXAGON=6, ISEA_TRIANGLE=3, ISEA_DIAMOND=4 };

+enum isea_address_form { ISEA_GEO, ISEA_Q2DI, ISEA_SEQNUM, ISEA_INTERLEAVE,

+    ISEA_PLANE, ISEA_Q2DD, ISEA_PROJTRI, ISEA_VERTEX2DD, ISEA_HEX

+};

+

+struct isea_dgg {

+    int polyhedron; /* ignored, icosahedron */

+    double  o_lat, o_lon, o_az; /* orientation, radians */

+    int pole; /* true if standard snyder */

+    int topology; /* ignored, hexagon */

+    int aperture; /* valid values depend on partitioning method */

+    int resolution;

+    double  radius; /* radius of the earth in meters, ignored 1.0 */

+    int output; /* an isea_address_form */

+    int triangle; /* triangle of last transformed point */

+    int quad; /* quad of last transformed point */

+    unsigned long serial;

+};


+

+struct isea_pt {

+    double x, y;

+};

+

+struct isea_geo {

+    double lon, lat;

+};

+

+struct isea_address {

+    int type; /* enum isea_address_form */

+    int number;

+    double  x,y; /* or i,j or lon,lat depending on type */

+};

+

+/* ENDINC */

+

+enum snyder_polyhedron {

+    SNYDER_POLY_HEXAGON, SNYDER_POLY_PENTAGON,

+    SNYDER_POLY_TETRAHEDRON, SNYDER_POLY_CUBE,

+    SNYDER_POLY_OCTAHEDRON, SNYDER_POLY_DODECAHEDRON,

+    SNYDER_POLY_ICOSAHEDRON

+};

+

+struct snyder_constants {

+    double          g, G, theta, ea_w, ea_a, ea_b, g_w, g_a, g_b;

+};

+

+/* TODO put these in radians to avoid a later conversion */

+ISEA_STATIC

+struct snyder_constants constants[] = {

+    {23.80018260, 62.15458023, 60.0, 3.75, 1.033, 0.968, 5.09, 1.195, 1.0},

+    {20.07675127, 55.69063953, 54.0, 2.65, 1.030, 0.983, 3.59, 1.141, 1.027},

+    {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},

+    {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},

+    {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},

+    {0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0, 0.0},

+    {37.37736814, 36.0, 30.0, 17.27, 1.163, 0.860, 13.14, 1.584, 1.0},

+};

+

+#define E 52.62263186

+#define F 10.81231696

+

+#define DEG60 1.04719755119659774614

+#define DEG120 2.09439510239319549229

+#define DEG72 1.25663706143591729537

+#define DEG90 1.57079632679489661922

+#define DEG144 2.51327412287183459075

+#define DEG36 0.62831853071795864768

+#define DEG108 1.88495559215387594306

+#define DEG180 M_PI

+/* sqrt(5)/M_PI */

+#define ISEA_SCALE 0.8301572857837594396028083

+

+/* 26.565051177 degrees */

+#define V_LAT 0.46364760899944494524

+

+#define RAD2DEG (180.0/M_PI)

+#define DEG2RAD (M_PI/180.0)

+

+ISEA_STATIC

+struct isea_geo vertex[] = {

+    {0.0, DEG90},

+    {DEG180, V_LAT},

+    {-DEG108, V_LAT},

+    {-DEG36, V_LAT},

+    {DEG36, V_LAT},

+    {DEG108, V_LAT},

+    {-DEG144, -V_LAT},

+    {-DEG72, -V_LAT},

+    {0.0, -V_LAT},

+    {DEG72, -V_LAT},

+    {DEG144, -V_LAT},

+    {0.0, -DEG90}

+};

+

+/* TODO make an isea_pt array of the vertices as well */

+

+static int      tri_v1[] = {0, 0, 0, 0, 0, 0, 6, 7, 8, 9, 10, 2, 3, 4, 5, 1, 11, 11, 11, 11, 11};

+

+/* 52.62263186 */

+#define E_RAD 0.91843818702186776133

+

+/* 10.81231696 */

+#define F_RAD 0.18871053072122403508

+

+/* triangle Centers */

+struct isea_geo icostriangles[] = {

+    {0.0, 0.0},

+    {-DEG144, E_RAD},

+    {-DEG72, E_RAD},

+    {0.0, E_RAD},

+    {DEG72, E_RAD},

+    {DEG144, E_RAD},

+    {-DEG144, F_RAD},

+    {-DEG72, F_RAD},

+    {0.0, F_RAD},

+    {DEG72, F_RAD},

+    {DEG144, F_RAD},

+    {-DEG108, -F_RAD},

+    {-DEG36, -F_RAD},

+    {DEG36, -F_RAD},

+    {DEG108, -F_RAD},

+    {DEG180, -F_RAD},

+    {-DEG108, -E_RAD},

+    {-DEG36, -E_RAD},

+    {DEG36, -E_RAD},

+    {DEG108, -E_RAD},

+    {DEG180, -E_RAD},

+};

+

+static double

+az_adjustment(int triangle)

+{

+    double          adj;

+

+    struct isea_geo v;

+    struct isea_geo c;

+

+    v = vertex[tri_v1[triangle]];

+    c = icostriangles[triangle];

+

+    /* TODO looks like the adjustment is always either 0 or 180 */

+    /* at least if you pick your vertex carefully */

+    adj = atan2(cos(v.lat) * sin(v.lon - c.lon),

+            cos(c.lat) * sin(v.lat)

+            - sin(c.lat) * cos(v.lat) * cos(v.lon - c.lon));

+    return adj;

+}

+

+/* R tan(g) sin(60) */

+#define TABLE_G 0.6615845383

+

+/* H = 0.25 R tan g = */

+#define TABLE_H 0.1909830056

+

+#define RPRIME 0.91038328153090290025

+

+ISEA_STATIC

+struct isea_pt

+isea_triangle_xy(int triangle)

+{

+    struct isea_pt  c;

+    double Rprime = 0.91038328153090290025;

+

+    triangle = (triangle - 1) % 20;

+

+    c.x = TABLE_G * ((triangle % 5) - 2) * 2.0;

+    if (triangle > 9) {

+        c.x += TABLE_G;

+    }

+    switch (triangle / 5) {

+    case 0:

+        c.y = 5.0 * TABLE_H;

+        break;

+    case 1:

+        c.y = TABLE_H;

+        break;

+    case 2:

+        c.y = -TABLE_H;

+        break;

+    case 3:

+        c.y = -5.0 * TABLE_H;

+        break;

+    default:

+        /* should be impossible */

+        exit(EXIT_FAILURE);

+    };

+    c.x *= Rprime;

+    c.y *= Rprime;

+

+    return c;

+}

+

+/* snyder eq 14 */

+static double

+sph_azimuth(double f_lon, double f_lat, double t_lon, double t_lat)

+{

+    double          az;

+

+    az = atan2(cos(t_lat) * sin(t_lon - f_lon),

+           cos(f_lat) * sin(t_lat)

+           - sin(f_lat) * cos(t_lat) * cos(t_lon - f_lon)

+        );

+    return az;

+}

+

+/* coord needs to be in radians */

+ISEA_STATIC

+int

+isea_snyder_forward(struct isea_geo * ll, struct isea_pt * out)

+{

+    int             i;

+

+    /*

+     * spherical distance from center of polygon face to any of its

+     * vertexes on the globe

+     */

+    double          g;

+

+    /*

+     * spherical angle between radius vector to center and adjacent edge

+     * of spherical polygon on the globe

+     */

+    double          G;

+

+    /*

+     * plane angle between radius vector to center and adjacent edge of

+     * plane polygon

+     */

+    double          theta;

+

+    /* additional variables from snyder */

+    double          q, Rprime, H, Ag, Azprime, Az, dprime, f, rho,

+                    x, y;

+

+    /* variables used to store intermediate results */

+    double          cot_theta, tan_g, az_offset;

+

+    /* how many multiples of 60 degrees we adjust the azimuth */

+    int             Az_adjust_multiples;

+

+    struct snyder_constants c;

+

+    /*

+     * TODO by locality of reference, start by trying the same triangle

+     * as last time

+     */

+

+    /* TODO put these constants in as radians to begin with */

+    c = constants[SNYDER_POLY_ICOSAHEDRON];

+    theta = c.theta * DEG2RAD;

+    g = c.g * DEG2RAD;

+    G = c.G * DEG2RAD;

+

+    for (i = 1; i <= 20; i++) {

+        double          z;

+        struct isea_geo center;

+

+        center = icostriangles[i];

+

+        /* step 1 */

+        z = acos(sin(center.lat) * sin(ll->lat)

+             + cos(center.lat) * cos(ll->lat) * cos(ll->lon - center.lon));

+        /* not on this triangle */

+        if (z > g + 0.000005) { /* TODO DBL_EPSILON */

+            continue;

+        }

+

+        Az = sph_azimuth(center.lon, center.lat, ll->lon, ll->lat);

+

+        /* step 2 */

+

+        /* This calculates "some" vertex coordinate */

+        az_offset = az_adjustment(i);

+

+        Az -= az_offset;

+

+        /* TODO I don't know why we do this.  It's not in snyder */

+        /* maybe because we should have picked a better vertex */

+        if (Az < 0.0) {

+            Az += 2.0 * M_PI;

+        }

+        /*

+         * adjust Az for the point to fall within the range of 0 to

+         * 2(90 - theta) or 60 degrees for the hexagon, by

+         * and therefore 120 degrees for the triangle

+         * of the icosahedron

+         * subtracting or adding multiples of 60 degrees to Az and

+         * recording the amount of adjustment

+         */

+

+        Az_adjust_multiples = 0;

+        while (Az < 0.0) {

+            Az += DEG120;

+            Az_adjust_multiples--;

+        }

+        while (Az > DEG120 + DBL_EPSILON) {

+            Az -= DEG120;

+            Az_adjust_multiples++;

+        }

+

+        /* step 3 */

+        cot_theta = 1.0 / tan(theta);

+        tan_g = tan(g); /* TODO this is a constant */

+

+        /* Calculate q from eq 9. */

+        /* TODO cot_theta is cot(30) */

+        q = atan2(tan_g, cos(Az) + sin(Az) * cot_theta);

+

+        /* not in this triangle */

+        if (z > q + 0.000005) {

+            continue;

+        }

+        /* step 4 */

+

+        /* Apply equations 5-8 and 10-12 in order */

+

+        /* eq 5 */

+        /* Rprime = 0.9449322893 * R; */

+        /* R' in the paper is for the truncated */

+        Rprime = 0.91038328153090290025;

+

+        /* eq 6 */

+        H = acos(sin(Az) * sin(G) * cos(g) - cos(Az) * cos(G));

+

+        /* eq 7 */

+        /* Ag = (Az + G + H - DEG180) * M_PI * R * R / DEG180; */

+        Ag = Az + G + H - DEG180;

+

+        /* eq 8 */

+        Azprime = atan2(2.0 * Ag, Rprime * Rprime * tan_g * tan_g - 2.0 * Ag * cot_theta);

+

+        /* eq 10 */

+        /* cot(theta) = 1.73205080756887729355 */

+        dprime = Rprime * tan_g / (cos(Azprime) + sin(Azprime) * cot_theta);

+

+        /* eq 11 */

+        f = dprime / (2.0 * Rprime * sin(q / 2.0));

+

+        /* eq 12 */

+        rho = 2.0 * Rprime * f * sin(z / 2.0);

+

+        /*

+         * add back the same 60 degree multiple adjustment from step

+         * 2 to Azprime

+         */

+

+        Azprime += DEG120 * Az_adjust_multiples;

+

+        /* calculate rectangular coordinates */

+

+        x = rho * sin(Azprime);

+        y = rho * cos(Azprime);

+

+        /*

+         * TODO

+         * translate coordinates to the origin for the particular

+         * hexagon on the flattened polyhedral map plot

+         */

+

+        out->x = x;

+        out->y = y;

+

+        return i;

+    }

+

+    /*

+     * should be impossible, this implies that the coordinate is not on

+     * any triangle

+     */

+

+    fprintf(stderr, "impossible transform: %f %f is not on any triangle\n",

+        ll->lon * RAD2DEG, ll->lat * RAD2DEG);

+

+    exit(EXIT_FAILURE);

+

+    /* not reached */

+    return 0;       /* supresses a warning */

+}

+

+/*

+ * return the new coordinates of any point in orginal coordinate system.

+ * Define a point (newNPold) in orginal coordinate system as the North Pole in

+ * new coordinate system, and the great circle connect the original and new

+ * North Pole as the lon0 longitude in new coordinate system, given any point

+ * in orginal coordinate system, this function return the new coordinates.

+ */

+

+#define PRECISION 0.0000000000005

+

+/* formula from Snyder, Map Projections: A working manual, p31 */

+/*

+ * old north pole at np in new coordinates

+ * could be simplified a bit with fewer intermediates

+ *

+ * TODO take a result pointer

+ */

+ISEA_STATIC

+struct isea_geo

+snyder_ctran(struct isea_geo * np, struct isea_geo * pt)

+{

+    struct isea_geo npt;

+    double          alpha, phi, lambda, lambda0, beta, lambdap, phip;

+    double          sin_phip;

+    double          lp_b;   /* lambda prime minus beta */

+    double          cos_p, sin_a;

+

+    phi = pt->lat;

+    lambda = pt->lon;

+    alpha = np->lat;

+    beta = np->lon;

+    lambda0 = beta;

+

+    cos_p = cos(phi);

+    sin_a = sin(alpha);

+

+    /* mpawm 5-7 */

+    sin_phip = sin_a * sin(phi) - cos(alpha) * cos_p * cos(lambda - lambda0);

+

+    /* mpawm 5-8b */

+

+    /* use the two argument form so we end up in the right quadrant */

+    lp_b = atan2(cos_p * sin(lambda - lambda0),

+       (sin_a * cos_p * cos(lambda - lambda0) + cos(alpha) * sin(phi)));

+

+    lambdap = lp_b + beta;

+

+    /* normalize longitude */

+    /* TODO can we just do a modulus ? */

+    lambdap = fmod(lambdap, 2 * M_PI);

+    while (lambdap > M_PI)

+        lambdap -= 2 * M_PI;

+    while (lambdap < -M_PI)

+        lambdap += 2 * M_PI;

+

+    phip = asin(sin_phip);

+

+    npt.lat = phip;

+    npt.lon = lambdap;

+

+    return npt;

+}

+

+ISEA_STATIC

+struct isea_geo

+isea_ctran(struct isea_geo * np, struct isea_geo * pt, double lon0)

+{

+    struct isea_geo npt;

+

+    np->lon += M_PI;

+    npt = snyder_ctran(np, pt);

+    np->lon -= M_PI;

+

+    npt.lon -= (M_PI - lon0 + np->lon);

+

+    /*

+     * snyder is down tri 3, isea is along side of tri1 from vertex 0 to

+     * vertex 1 these are 180 degrees apart

+     */

+    npt.lon += M_PI;

+    /* normalize longitude */

+    npt.lon = fmod(npt.lon, 2 * M_PI);

+    while (npt.lon > M_PI)

+        npt.lon -= 2 * M_PI;

+    while (npt.lon < -M_PI)

+        npt.lon += 2 * M_PI;

+

+    return npt;

+}

+

+/* in radians */

+#define ISEA_STD_LAT 1.01722196792335072101

+#define ISEA_STD_LON .19634954084936207740

+

+/* fuller's at 5.2454 west, 2.3009 N, adjacent at 7.46658 deg */

+

+ISEA_STATIC

+int

+isea_grid_init(struct isea_dgg * g)

+{

+    if (!g)

+        return 0;

+

+    g->polyhedron = 20;

+    g->o_lat = ISEA_STD_LAT;

+    g->o_lon = ISEA_STD_LON;

+    g->o_az = 0.0;

+    g->aperture = 4;

+    g->resolution = 6;

+    g->radius = 1.0;

+    g->topology = 6;

+

+    return 1;

+}

+

+ISEA_STATIC

+int

+isea_orient_isea(struct isea_dgg * g)

+{

+    if (!g)

+        return 0;

+    g->o_lat = ISEA_STD_LAT;

+    g->o_lon = ISEA_STD_LON;

+    g->o_az = 0.0;

+    return 1;

+}

+

+ISEA_STATIC

+int

+isea_orient_pole(struct isea_dgg * g)

+{

+    if (!g)

+        return 0;

+    g->o_lat = M_PI / 2.0;

+    g->o_lon = 0.0;

+    g->o_az = 0;

+    return 1;

+}

+

+ISEA_STATIC

+int

+isea_transform(struct isea_dgg * g, struct isea_geo * in,

+           struct isea_pt * out)

+{

+    struct isea_geo i, pole;

+    int             tri;

+

+    pole.lat = g->o_lat;

+    pole.lon = g->o_lon;

+

+    i = isea_ctran(&pole, in, g->o_az);

+

+    tri = isea_snyder_forward(&i, out);

+    out->x *= g->radius;

+    out->y *= g->radius;

+    g->triangle = tri;

+

+    return tri;

+}

+

+#define DOWNTRI(tri) (((tri - 1) / 5) % 2 == 1)

+

+ISEA_STATIC

+void

+isea_rotate(struct isea_pt * pt, double degrees)

+{

+    double          rad;

+

+    double          x, y;

+

+    rad = -degrees * M_PI / 180.0;

+    while (rad >= 2.0 * M_PI) rad -= 2.0 * M_PI;

+    while (rad <= -2.0 * M_PI) rad += 2.0 * M_PI;

+

+    x = pt->x * cos(rad) + pt->y * sin(rad);

+    y = -pt->x * sin(rad) + pt->y * cos(rad);

+

+    pt->x = x;

+    pt->y = y;

+}

+

+ISEA_STATIC

+int isea_tri_plane(int tri, struct isea_pt *pt, double radius) {

+    struct isea_pt tc; /* center of triangle */

+

+    if (DOWNTRI(tri)) {

+        isea_rotate(pt, 180.0);

+    }

+    tc = isea_triangle_xy(tri);

+    tc.x *= radius;

+    tc.y *= radius;

+    pt->x += tc.x;

+    pt->y += tc.y;

+

+    return tri;

+}

+

+/* convert projected triangle coords to quad xy coords, return quad number */

+ISEA_STATIC

+int

+isea_ptdd(int tri, struct isea_pt *pt) {

+    int             downtri, quad;

+

+    downtri = (((tri - 1) / 5) % 2 == 1);

+    quad = ((tri - 1) % 5) + ((tri - 1) / 10) * 5 + 1;

+

+    isea_rotate(pt, downtri ? 240.0 : 60.0);

+    if (downtri) {

+        pt->x += 0.5;

+        /* pt->y += cos(30.0 * M_PI / 180.0); */

+        pt->y += .86602540378443864672;

+    }

+    return quad;

+}

+

+ISEA_STATIC

+int

+isea_dddi_ap3odd(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di)

+{

+    struct isea_pt  v;

+    double          hexwidth;

+    double          sidelength; /* in hexes */

+    int             d, i;

+    int             maxcoord;

+    struct hex      h;

+

+    /* This is the number of hexes from apex to base of a triangle */

+    sidelength = (pow(2.0, g->resolution) + 1.0) / 2.0;

+

+    /* apex to base is cos(30deg) */

+    hexwidth = cos(M_PI / 6.0) / sidelength;

+

+    /* TODO I think sidelength is always x.5, so

+     * (int)sidelength * 2 + 1 might be just as good

+     */

+    maxcoord = (int) (sidelength * 2.0 + 0.5);

+

+    v = *pt;

+    hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);

+    h.iso = 0;

+    hex_iso(&h);

+

+    d = h.x - h.z;

+    i = h.x + h.y + h.y;

+

+    /*

+     * you want to test for max coords for the next quad in the same

+     * "row" first to get the case where both are max

+     */

+    if (quad <= 5) {

+        if (d == 0 && i == maxcoord) {

+            /* north pole */

+            quad = 0;

+            d = 0;

+            i = 0;

+        } else if (i == maxcoord) {

+            /* upper right in next quad */

+            quad += 1;

+            if (quad == 6)

+                quad = 1;

+            i = maxcoord - d;

+            d = 0;

+        } else if (d == maxcoord) {

+            /* lower right in quad to lower right */

+            quad += 5;

+            d = 0;

+        }

+    } else if (quad >= 6) {

+        if (i == 0 && d == maxcoord) {

+            /* south pole */

+            quad = 11;

+            d = 0;

+            i = 0;

+        } else if (d == maxcoord) {

+            /* lower right in next quad */

+            quad += 1;

+            if (quad == 11)

+                quad = 6;

+            d = maxcoord - i;

+            i = 0;

+        } else if (i == maxcoord) {

+            /* upper right in quad to upper right */

+            quad = (quad - 4) % 5;

+            i = 0;

+        }

+    }

+

+    di->x = d;

+    di->y = i;

+

+    g->quad = quad;

+    return quad;

+}

+

+ISEA_STATIC

+int

+isea_dddi(struct isea_dgg *g, int quad, struct isea_pt *pt, struct isea_pt *di) {

+    struct isea_pt  v;

+    double          hexwidth;

+    int             sidelength; /* in hexes */

+    struct hex      h;

+

+    if (g->aperture == 3 && g->resolution % 2 != 0) {

+        return isea_dddi_ap3odd(g, quad, pt, di);

+    }

+    /* todo might want to do this as an iterated loop */

+    if (g->aperture >0) {

+        sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);

+    } else {

+        sidelength = g->resolution;

+    }

+

+    hexwidth = 1.0 / sidelength;

+

+    v = *pt;

+    isea_rotate(&v, -30.0);

+    hexbin2(hexwidth, v.x, v.y, &h.x, &h.y);

+    h.iso = 0;

+    hex_iso(&h);

+

+    /* we may actually be on another quad */

+    if (quad <= 5) {

+        if (h.x == 0 && h.z == -sidelength) {

+            /* north pole */

+            quad = 0;

+            h.z = 0;

+            h.y = 0;

+            h.x = 0;

+        } else if (h.z == -sidelength) {

+            quad = quad + 1;

+            if (quad == 6)

+                quad = 1;

+            h.y = sidelength - h.x;

+            h.z = h.x - sidelength;

+            h.x = 0;

+        } else if (h.x == sidelength) {

+            quad += 5;

+            h.y = -h.z;

+            h.x = 0;

+        }

+    } else if (quad >= 6) {

+        if (h.z == 0 && h.x == sidelength) {

+            /* south pole */

+            quad = 11;

+            h.x = 0;

+            h.y = 0;

+            h.z = 0;

+        } else if (h.x == sidelength) {

+            quad = quad + 1;

+            if (quad == 11)

+                quad = 6;

+            h.x = h.y + sidelength;

+            h.y = 0;

+            h.z = -h.x;

+        } else if (h.y == -sidelength) {

+            quad -= 4;

+            h.y = 0;

+            h.z = -h.x;

+        }

+    }

+    di->x = h.x;

+    di->y = -h.z;

+

+    g->quad = quad;

+    return quad;

+}

+

+ISEA_STATIC

+int isea_ptdi(struct isea_dgg *g, int tri, struct isea_pt *pt,

+            struct isea_pt *di) {

+    struct isea_pt  v;

+    int             quad;

+

+    v = *pt;

+    quad = isea_ptdd(tri, &v);

+    quad = isea_dddi(g, quad, &v, di);

+    return quad;

+}

+

+/* q2di to seqnum */

+ISEA_STATIC

+int isea_disn(struct isea_dgg *g, int quad, struct isea_pt *di) {

+    int             sidelength;

+    int             sn, height;

+    int             hexes;

+

+    if (quad == 0) {

+        g->serial = 1;

+        return g->serial;

+    }

+    /* hexes in a quad */

+    hexes = (int) (pow(g->aperture, g->resolution) + 0.5);

+    if (quad == 11) {

+        g->serial = 1 + 10 * hexes + 1;

+        return g->serial;

+    }

+    if (g->aperture == 3 && g->resolution % 2 == 1) {

+        height = (int) (pow(g->aperture, (g->resolution - 1) / 2.0));

+        sn = ((int) di->x) * height;

+        sn += ((int) di->y) / height;

+        sn += (quad - 1) * hexes;

+        sn += 2;

+    } else {

+        sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);

+        sn = (quad - 1) * hexes + sidelength * di->x + di->y + 2;

+    }

+

+    g->serial = sn;

+    return sn;

+}

+

+/* TODO just encode the quad in the d or i coordinate

+ * quad is 0-11, which can be four bits.

+ * d' = d << 4 + q, d = d' >> 4, q = d' & 0xf

+ */

+/* convert a q2di to global hex coord */

+ISEA_STATIC

+int isea_hex(struct isea_dgg *g, int tri,

+        struct isea_pt *pt, struct isea_pt *hex) {

+    struct isea_pt v;

+    int sidelength;

+    int d, i, x, y, quad;

+

+    quad = isea_ptdi(g, tri, pt, &v);

+

+    hex->x = ((int)v.x << 4) + quad;

+    hex->y = v.y;

+

+    return 1;

+

+    d = v.x;

+    i = v.y;

+

+    /* Aperture 3 odd resolutions */

+    if (g->aperture == 3 && g->resolution % 2 != 0) {

+        int offset = (int)(pow(3.0, g->resolution - 1) + 0.5);

+

+        d += offset * ((g->quad-1) % 5);

+        i += offset * ((g->quad-1) % 5);

+

+        if (quad == 0) {

+            d = 0;

+            i = offset;

+        } else if (quad == 11) {

+            d = 2 * offset;

+            i = 0;

+        } else if (quad > 5) {

+            d += offset;

+        }

+

+        x = (2*d - i) /3;

+        y = (2*i - d) /3;

+

+        hex->x = x + offset / 3;

+        hex->y = y + 2 * offset / 3;

+        return 1;

+    }

+

+    /* aperture 3 even resolutions and aperture 4 */

+    sidelength = (int) (pow(g->aperture, g->resolution / 2.0) + 0.5);

+    if (g->quad == 0) {

+        hex->x = 0;

+        hex->y = sidelength;

+    } else if (g->quad == 11) {

+        hex->x = sidelength * 2;

+        hex->y = 0;

+    } else {

+        hex->x = d + sidelength * ((g->quad-1) % 5);

+        if (g->quad > 5) hex->x += sidelength;

+        hex->y = i + sidelength * ((g->quad-1) % 5);

+    }

+

+    return 1;

+}

+

+ISEA_STATIC

+struct isea_pt

+isea_forward(struct isea_dgg *g, struct isea_geo *in)

+{

+    int             tri;

+    struct isea_pt  out, coord;

+

+    tri = isea_transform(g, in, &out);

+

+    if (g->output == ISEA_PLANE) {

+        isea_tri_plane(tri, &out, g->radius);

+        return out;

+    }

+

+    /* convert to isea standard triangle size */

+    out.x = out.x / g->radius * ISEA_SCALE;

+    out.y = out.y / g->radius * ISEA_SCALE;

+    out.x += 0.5;

+    out.y += 2.0 * .14433756729740644112;

+

+    switch (g->output) {

+    case ISEA_PROJTRI:

+        /* nothing to do, already in projected triangle */

+        break;

+    case ISEA_VERTEX2DD:

+        g->quad = isea_ptdd(tri, &out);

+        break;

+    case ISEA_Q2DD:

+        /* Same as above, we just don't print as much */

+        g->quad = isea_ptdd(tri, &out);

+        break;

+    case ISEA_Q2DI:

+        g->quad = isea_ptdi(g, tri, &out, &coord);

+        return coord;

+        break;

+    case ISEA_SEQNUM:

+        isea_ptdi(g, tri, &out, &coord);

+        /* disn will set g->serial */

+        isea_disn(g, g->quad, &coord);

+        return coord;

+        break;

+    case ISEA_HEX:

+        isea_hex(g, tri, &out, &coord);

+        return coord;

+        break;

+    }

+

+    return out;

+}

+/*

+ * Proj 4 integration code follows

+ */

+

+#define PJ_LIB__

+#include <projects.h>

+

+PROJ_HEAD(isea, "Icosahedral Snyder Equal Area") "\n\tSph";

+

+struct pj_opaque {

+    struct isea_dgg dgg;

+};

+

+

+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */

+    XY xy = {0.0,0.0};

+    struct pj_opaque *Q = P->opaque;

+    struct isea_pt out;

+    struct isea_geo in;

+

+    in.lon = lp.lam;

+    in.lat = lp.phi;

+

+    out = isea_forward(&Q->dgg, &in);

+

+    xy.x = out.x;

+    xy.y = out.y;

+

+    return xy;

+}

+

+

+static void *freeup_new (PJ *P) {                       /* Destructor */

+    if (0==P)

+        return 0;

+    if (0==P->opaque)

+        return pj_dealloc (P);

+

+    pj_dealloc (P->opaque);

+    return pj_dealloc(P);

+}

+

+static void freeup (PJ *P) {

+    freeup_new (P);

+    return;

+}

+

+

+PJ *PROJECTION(isea) {

+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));

+    if (0==Q)

+        return freeup_new (P);

+    P->opaque = Q;

+

+    char *opt;

+

+    P->fwd = s_forward;

+    isea_grid_init(&Q->dgg);

+

+    Q->dgg.output = ISEA_PLANE;

+/*      P->dgg.radius = P->a; / * otherwise defaults to 1 */

+    /* calling library will scale, I think */

+

+    opt = pj_param(P->ctx,P->params, "sorient").s;

+    if (opt) {

+        if (!strcmp(opt, "isea")) {

+            isea_orient_isea(&Q->dgg);

+        } else if (!strcmp(opt, "pole")) {

+            isea_orient_pole(&Q->dgg);

+        } else {

+            E_ERROR(-34);

+        }

+    }

+

+    if (pj_param(P->ctx,P->params, "tazi").i) {

+        Q->dgg.o_az = pj_param(P->ctx,P->params, "razi").f;

+    }

+

+    if (pj_param(P->ctx,P->params, "tlon_0").i) {

+        Q->dgg.o_lon = pj_param(P->ctx,P->params, "rlon_0").f;

+    }

+

+    if (pj_param(P->ctx,P->params, "tlat_0").i) {

+        Q->dgg.o_lat = pj_param(P->ctx,P->params, "rlat_0").f;

+    }

+

+    if (pj_param(P->ctx,P->params, "taperture").i) {

+        Q->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;

+    }

+

+    if (pj_param(P->ctx,P->params, "tresolution").i) {

+        Q->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;

+    }

+

+    opt = pj_param(P->ctx,P->params, "smode").s;

+    if (opt) {

+        if (!strcmp(opt, "plane")) {

+            Q->dgg.output = ISEA_PLANE;

+        } else if (!strcmp(opt, "di")) {

+            Q->dgg.output = ISEA_Q2DI;

+        }

+        else if (!strcmp(opt, "dd")) {

+            Q->dgg.output = ISEA_Q2DD;

+        }

+        else if (!strcmp(opt, "hex")) {

+            Q->dgg.output = ISEA_HEX;

+        }

+        else {

+            /* TODO verify error code.  Possibly eliminate magic */

+            E_ERROR(-34);

+        }

+    }

+

+    if (pj_param(P->ctx,P->params, "trescale").i) {

+        Q->dgg.radius = ISEA_SCALE;

+    }

+

+    if (pj_param(P->ctx,P->params, "tresolution").i) {

+        Q->dgg.resolution = pj_param(P->ctx,P->params, "iresolution").i;

+    } else {

+        Q->dgg.resolution = 4;

+    }

+

+    if (pj_param(P->ctx,P->params, "taperture").i) {

+        Q->dgg.aperture = pj_param(P->ctx,P->params, "iaperture").i;

+    } else {

+        Q->dgg.aperture = 3;

+    }

+

+    return P;

+}

+

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_isea_selftest (void) {return 0;}

+#else

+

+int pj_isea_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

+    char s_args[] = {"+proj=isea   +a=6400000    +lat_1=0.5 +lat_2=2"};

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+

+    XY s_fwd_expect[] = {

+        {-1097074.9480224741, 3442909.3090371834},

+        {-1097074.9482647954, 3233611.7285857084},

+        {-1575486.3536415542, 3442168.3420281881},

+        {-1575486.353880283,  3234352.6955947056},

+    };

+

+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0);

+}

+

+

+#endif

diff --git a/src/PJ_krovak.c b/src/PJ_krovak.c
index dd250134..bb714189 100644
--- a/src/PJ_krovak.c
+++ b/src/PJ_krovak.c
@@ -29,220 +29,284 @@
  * SOFTWARE.
  *****************************************************************************/
 
-#define PROJ_PARMS__ \
-	double	C_x;
 #define PJ_LIB__
-
 #include <projects.h>
 #include <string.h>
 #include <stdio.h>
 
 PROJ_HEAD(krovak, "Krovak") "\n\tPCyl., Ellps.";
 
+struct pj_opaque {
+    double  C_x;
+};
+
 /**
    NOTES: According to EPSG the full Krovak projection method should have
           the following parameters.  Within PROJ.4 the azimuth, and pseudo
-          standard parallel are hardcoded in the algorithm and can't be 
+          standard parallel are hardcoded in the algorithm and can't be
           altered from outside.  The others all have defaults to match the
           common usage with Krovak projection.
 
   lat_0 = latitude of centre of the projection
-         
+
   lon_0 = longitude of centre of the projection
-  
+
   ** = azimuth (true) of the centre line passing through the centre of the projection
 
   ** = latitude of pseudo standard parallel
-   
+
   k  = scale factor on the pseudo standard parallel
-  
+
   x_0 = False Easting of the centre of the projection at the apex of the cone
-  
+
   y_0 = False Northing of the centre of the projection at the apex of the cone
 
  **/
 
 
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    /* calculate xy from lat/lon */
 
-FORWARD(e_forward); /* ellipsoid */
-/* calculate xy from lat/lon */
+    /* Constants, identical to inverse transform function */
+    double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
+    double gfi, u, fi0, deltav, s, d, eps, ro;
 
-/* Constants, identical to inverse transform function */
-	double s45, s90, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
-	double gfi, u, fi0, deltav, s, d, eps, ro;
 
+    s45 = 0.785398163397448;    /* 45deg */
+    s90 = 2 * s45;
+    fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
 
-	s45 = 0.785398163397448;    /* 45deg */
-	s90 = 2 * s45;
-	fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
-
-   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 
+    /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
       be set to 1 here.
       Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
       e2=0.006674372230614;
-   */
-	a =  1; /* 6377397.155; */
-	/* e2 = P->es;*/       /* 0.006674372230614; */
-	e2 = 0.006674372230614;
-	e = sqrt(e2);
-
-	alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
+    */
+    a =  1; /* 6377397.155; */
+    /* e2 = P->es;*/       /* 0.006674372230614; */
+    e2 = 0.006674372230614;
+    e = sqrt(e2);
 
-	uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
-	u0 = asin(sin(fi0) / alfa);
-	g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
+    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
 
-	k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
+    uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
+    u0 = asin(sin(fi0) / alfa);
+    g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
 
-	k1 = P->k0;
-	n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
-	s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
-	n = sin(s0);
-	ro0 = k1 * n0 / tan(s0);
-	ad = s90 - uq;
+    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
 
-/* Transformation */
+    k1 = P->k0;
+    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
+    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
+    n = sin(s0);
+    ro0 = k1 * n0 / tan(s0);
+    ad = s90 - uq;
 
-	gfi =pow ( ((1. + e * sin(lp.phi)) /
-               (1. - e * sin(lp.phi))) , (alfa * e / 2.));
+    /* Transformation */
+    gfi = pow ( ((1. + e * sin(lp.phi)) /
+                (1. - e * sin(lp.phi))) , (alfa * e / 2.));
 
-	u= 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
+    u = 2. * (atan(k * pow( tan(lp.phi / 2. + s45), alfa) / gfi)-s45);
 
-	deltav = - lp.lam * alfa;
+    deltav = - lp.lam * alfa;
 
-	s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
-	d = asin(cos(u) * sin(deltav) / cos(s));
-	eps = n * d;
-	ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n)   ;
+    s = asin(cos(ad) * sin(u) + sin(ad) * cos(u) * cos(deltav));
+    d = asin(cos(u) * sin(deltav) / cos(s));
+    eps = n * d;
+    ro = ro0 * pow(tan(s0 / 2. + s45) , n) / pow(tan(s / 2. + s45) , n);
 
    /* x and y are reverted! */
-	xy.y = ro * cos(eps) / a;
-	xy.x = ro * sin(eps) / a;
+    xy.y = ro * cos(eps) / a;
+    xy.x = ro * sin(eps) / a;
 
-        if( !pj_param(P->ctx, P->params, "tczech").i )
-	  {
-	    xy.y *= -1.0;
-	    xy.x *= -1.0;
-	  }
+    if( !pj_param(P->ctx, P->params, "tczech").i ) {
+        xy.y *= -1.0;
+        xy.x *= -1.0;
+    }
 
-	return (xy);
+    return xy;
 }
 
 
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    /* calculate lat/lon from xy */
 
-INVERSE(e_inverse); /* ellipsoid */
-	/* calculate lat/lon from xy */
+    /* Constants, identisch wie in der Umkehrfunktion */
+    double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
+    double u, deltav, s, d, eps, ro, fi1, xy0;
+    int ok;
 
-/* Constants, identisch wie in der Umkehrfunktion */
-	double s45, s90, fi0, e2, e, alfa, uq, u0, g, k, k1, n0, ro0, ad, a, s0, n;
-	double u, deltav, s, d, eps, ro, fi1, xy0;
-	int ok;
+    s45 = 0.785398163397448;    /* 45deg */
+    s90 = 2 * s45;
+    fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
 
-	s45 = 0.785398163397448;    /* 45deg */
-	s90 = 2 * s45;
-	fi0 = P->phi0;    /* Latitude of projection centre 49deg 30' */
 
-
-   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must 
+   /* Ellipsoid is used as Parameter in for.c and inv.c, therefore a must
       be set to 1 here.
       Ellipsoid Bessel 1841 a = 6377397.155m 1/f = 299.1528128,
       e2=0.006674372230614;
    */
-	a = 1; /* 6377397.155; */
-	/* e2 = P->es; */      /* 0.006674372230614; */
-	e2 = 0.006674372230614;
-	e = sqrt(e2);
+    a = 1; /* 6377397.155; */
+    /* e2 = P->es; */      /* 0.006674372230614; */
+    e2 = 0.006674372230614;
+    e = sqrt(e2);
 
-	alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
-	uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
-	u0 = asin(sin(fi0) / alfa);
-	g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
+    alfa = sqrt(1. + (e2 * pow(cos(fi0), 4)) / (1. - e2));
+    uq = 1.04216856380474;      /* DU(2, 59, 42, 42.69689) */
+    u0 = asin(sin(fi0) / alfa);
+    g = pow(   (1. + e * sin(fi0)) / (1. - e * sin(fi0)) , alfa * e / 2.  );
 
-	k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
+    k = tan( u0 / 2. + s45) / pow  (tan(fi0 / 2. + s45) , alfa) * g;
 
-	k1 = P->k0;
-	n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
-	s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
-	n = sin(s0);
-	ro0 = k1 * n0 / tan(s0);
-	ad = s90 - uq;
+    k1 = P->k0;
+    n0 = a * sqrt(1. - e2) / (1. - e2 * pow(sin(fi0), 2));
+    s0 = 1.37008346281555;       /* Latitude of pseudo standard parallel 78deg 30'00" N */
+    n = sin(s0);
+    ro0 = k1 * n0 / tan(s0);
+    ad = s90 - uq;
 
 
-/* Transformation */
+    /* Transformation */
    /* revert y, x*/
-	xy0=xy.x;
-	xy.x=xy.y;
-	xy.y=xy0;
-
-        if( !pj_param(P->ctx, P->params, "tczech").i )
-	  {
-	    xy.x *= -1.0;
-	    xy.y *= -1.0;
-	  }
+    xy0 = xy.x;
+    xy.x = xy.y;
+    xy.y = xy0;
 
-	ro = sqrt(xy.x * xy.x + xy.y * xy.y);
-	eps = atan2(xy.y, xy.x);
-	d = eps / sin(s0);
-	s = 2. * (atan(  pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
+    if( !pj_param(P->ctx, P->params, "tczech").i ) {
+        xy.x *= -1.0;
+        xy.y *= -1.0;
+    }
 
-	u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
-	deltav = asin(cos(s) * sin(d) / cos(u));
+    ro = sqrt(xy.x * xy.x + xy.y * xy.y);
+    eps = atan2(xy.y, xy.x);
+    d = eps / sin(s0);
+    s = 2. * (atan(  pow(ro0 / ro, 1. / n) * tan(s0 / 2. + s45)) - s45);
 
-	lp.lam = P->lam0 - deltav / alfa;
+    u = asin(cos(ad) * sin(s) - sin(ad) * cos(s) * cos(d));
+    deltav = asin(cos(s) * sin(d) / cos(u));
 
-/* ITERATION FOR lp.phi */
-   fi1 = u;
+    lp.lam = P->lam0 - deltav / alfa;
 
-   ok = 0;
-   do
-   {
-   	lp.phi = 2. * ( atan( pow( k, -1. / alfa)  *
-                            pow( tan(u / 2. + s45) , 1. / alfa)  *
-                            pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
-                           )  - s45);
+    /* ITERATION FOR lp.phi */
+    fi1 = u;
 
-      if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1;
-      fi1 = lp.phi;
+    ok = 0;
+    do {
+        lp.phi = 2. * ( atan( pow( k, -1. / alfa)  *
+                              pow( tan(u / 2. + s45) , 1. / alfa)  *
+                              pow( (1. + e * sin(fi1)) / (1. - e * sin(fi1)) , e / 2.)
+                            )  - s45);
 
-   }
-   while (ok==0);
+        if (fabs(fi1 - lp.phi) < 0.000000000000001) ok=1;
+        fi1 = lp.phi;
+   } while (ok==0);
 
    lp.lam -= P->lam0;
 
-   return (lp);
+   return lp;
 }
 
-FREEUP; if (P) pj_dalloc(P); }
 
-ENTRY0(krovak)
-	double ts;
-	/* read some Parameters,
-	 * here Latitude Truescale */
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
 
-	ts = pj_param(P->ctx, P->params, "rlat_ts").f;
-	P->C_x = ts;
-	
-	/* we want Bessel as fixed ellipsoid */
-	P->a = 6377397.155;
-	P->e = sqrt(P->es = 0.006674372230614);
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(krovak) {
+    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_krovak;
+    double ts;
+    /* read some Parameters,
+     * here Latitude Truescale */
+
+    ts = pj_param(P->ctx, P->params, "rlat_ts").f;
+    Q->C_x = ts;
+
+    /* we want Bessel as fixed ellipsoid */
+    P->a = 6377397.155;
+    P->e = sqrt(P->es = 0.006674372230614);
 
         /* if latitude of projection center is not set, use 49d30'N */
-	if (!pj_param(P->ctx, P->params, "tlat_0").i)
-            P->phi0 = 0.863937979737193; 
+    if (!pj_param(P->ctx, P->params, "tlat_0").i)
+            P->phi0 = 0.863937979737193;
 
         /* if center long is not set use 42d30'E of Ferro - 17d40' for Ferro */
         /* that will correspond to using longitudes relative to greenwich    */
         /* as input and output, instead of lat/long relative to Ferro */
-	if (!pj_param(P->ctx, P->params, "tlon_0").i)
+    if (!pj_param(P->ctx, P->params, "tlon_0").i)
             P->lam0 = 0.7417649320975901 - 0.308341501185665;
 
         /* if scale not set default to 0.9999 */
-	if (!pj_param(P->ctx, P->params, "tk").i)
+    if (!pj_param(P->ctx, P->params, "tk").i)
             P->k0 = 0.9999;
 
-	/* always the same */
-        P->inv = e_inverse; 
-	P->fwd = e_forward;
+    /* always the same */
+    P->inv = e_inverse;
+    P->fwd = e_forward;
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_krovak_selftest (void) {return 0;}
+#else
+
+int pj_krovak_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=krovak   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {-3196535.2325636409,  -6617878.8675514441},
+        {-3260035.4405521089,  -6898873.6148780314},
+        {-3756305.3288691747,  -6478142.5615715114},
+        {-3831703.6585019818,  -6759107.1701553948},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {24.836218918719162,  59.758403933233858},
+        {24.836315484509566,  59.756888425730189},
+        {24.830447747947495,  59.758403933233858},
+        {24.830351182157091,  59.756888425730189},
+    };
+
+    return pj_generic_selftest (e_args, 0, tolerance_xy, tolerance_lp, 4, 4, fwd_in, e_fwd_expect, 0, inv_in, e_inv_expect, 0);
+}
 
-ENDENTRY(P)
 
+#endif