Added regression tests for projections beginning with a and b

Continuing alphabetically, with a few detours due to the occasional case
of more than one projection in one source file.

8 files changed, 949 insertions, 274 deletions

diff --git a/src/PJ_airy.c b/src/PJ_airy.c
index 73a51e5b..0addf2b8 100644
--- a/src/PJ_airy.c
+++ b/src/PJ_airy.c
@@ -1,7 +1,8 @@
 /******************************************************************************
  * Project:  PROJ.4
  * Purpose:  Implementation of the airy (Airy) projection.
- * Author:   Gerald Evenden
+ * Author:   Gerald Evenden (1995)
+ *           Thomas Knudsen (2016) - revise/add regression tests
  *
  ******************************************************************************
  * Copyright (c) 1995, Gerald Evenden
@@ -25,96 +26,163 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#define PROJ_PARMS__ \
-	double	p_halfpi; \
-	double	sinph0; \
-	double	cosph0; \
-	double	Cb; \
-	int		mode; \
-	int		no_cut;	/* do not cut at hemisphere limit */
 #define PJ_LIB__
 #include <projects.h>
 
 PROJ_HEAD(airy, "Airy") "\n\tMisc Sph, no inv.\n\tno_cut lat_b=";
 
+
+struct pj_opaque {
+	double	p_halfpi;
+	double	sinph0;
+	double	cosph0;
+	double	Cb;
+	int		mode;
+	int		no_cut;	/* do not cut at hemisphere limit */
+};
+
+
 # define EPS 1.e-10
 # define N_POLE	0
 # define S_POLE 1
 # define EQUIT	2
 # define OBLIQ	3
-FORWARD(s_forward); /* spheroid */
+
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double  sinlam, coslam, cosphi, sinphi, t, s, Krho, cosz;
 
 	sinlam = sin(lp.lam);
 	coslam = cos(lp.lam);
-	switch (P->mode) {
+	switch (Q->mode) {
 	case EQUIT:
 	case OBLIQ:
 		sinphi = sin(lp.phi);
 		cosphi = cos(lp.phi);
 		cosz = cosphi * coslam;
-		if (P->mode == OBLIQ)
-			cosz = P->sinph0 * sinphi + P->cosph0 * cosz;
-		if (!P->no_cut && cosz < -EPS)
+		if (Q->mode == OBLIQ)
+			cosz = Q->sinph0 * sinphi + Q->cosph0 * cosz;
+		if (!Q->no_cut && cosz < -EPS)
 			F_ERROR;
 		if (fabs(s = 1. - cosz) > EPS) {
 			t = 0.5 * (1. + cosz);
-			Krho = -log(t)/s - P->Cb / t;
+			Krho = -log(t)/s - Q->Cb / t;
 		} else
-			Krho = 0.5 - P->Cb;
+			Krho = 0.5 - Q->Cb;
 		xy.x = Krho * cosphi * sinlam;
-		if (P->mode == OBLIQ)
-			xy.y = Krho * (P->cosph0 * sinphi -
-				P->sinph0 * cosphi * coslam);
+		if (Q->mode == OBLIQ)
+			xy.y = Krho * (Q->cosph0 * sinphi -
+				Q->sinph0 * cosphi * coslam);
 		else
 			xy.y = Krho * sinphi;
 		break;
 	case S_POLE:
 	case N_POLE:
-		lp.phi = fabs(P->p_halfpi - lp.phi);
-		if (!P->no_cut && (lp.phi - EPS) > HALFPI)
+		lp.phi = fabs(Q->p_halfpi - lp.phi);
+		if (!Q->no_cut && (lp.phi - EPS) > HALFPI)
 			F_ERROR;
 		if ((lp.phi *= 0.5) > EPS) {
 			t = tan(lp.phi);
-			Krho = -2.*(log(cos(lp.phi)) / t + t * P->Cb);
+			Krho = -2.*(log(cos(lp.phi)) / t + t * Q->Cb);
 			xy.x = Krho * sinlam;
 			xy.y = Krho * coslam;
-			if (P->mode == N_POLE)
+			if (Q->mode == N_POLE)
 				xy.y = -xy.y;
 		} else
 			xy.x = xy.y = 0.;
 	}
-	return (xy);
+	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;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(airy)
+
+
+PJ *PROJECTION(airy) {
 	double beta;
 
-	P->no_cut = pj_param(P->ctx, P->params, "bno_cut").i;
+    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_airy;
+
+	Q->no_cut = pj_param(P->ctx, P->params, "bno_cut").i;
 	beta = 0.5 * (HALFPI - pj_param(P->ctx, P->params, "rlat_b").f);
 	if (fabs(beta) < EPS)
-		P->Cb = -0.5;
+		Q->Cb = -0.5;
 	else {
-		P->Cb = 1./tan(beta);
-		P->Cb *= P->Cb * log(cos(beta));
+		Q->Cb = 1./tan(beta);
+		Q->Cb *= Q->Cb * log(cos(beta));
 	}
+
 	if (fabs(fabs(P->phi0) - HALFPI) < EPS)
 		if (P->phi0 < 0.) {
-			P->p_halfpi = -HALFPI;
-			P->mode = S_POLE;
+			Q->p_halfpi = -HALFPI;
+			Q->mode = S_POLE;
 		} else {
-			P->p_halfpi =  HALFPI;
-			P->mode = N_POLE;
+			Q->p_halfpi =  HALFPI;
+			Q->mode = N_POLE;
 		}
 	else {
 		if (fabs(P->phi0) < EPS)
-			P->mode = EQUIT;
+			Q->mode = EQUIT;
 		else {
-			P->mode = OBLIQ;
-			P->sinph0 = sin(P->phi0);
-			P->cosph0 = cos(P->phi0);
+			Q->mode = OBLIQ;
+			Q->sinph0 = sin(P->phi0);
+			Q->cosph0 = cos(P->phi0);
 		}
 	}
 	P->fwd = s_forward;
 	P->es = 0.;
-ENDENTRY(P)
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_airy_selftest (void) {return 0;}
+#else
+
+int pj_airy_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=airy   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 189109.88690862127,   94583.752387504152},
+        { 189109.88690862127,  -94583.752387504152},
+        {-189109.88690862127,   94583.752387504152},
+        {-189109.88690862127,  -94583.752387504152},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 0, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+
+#endif
diff --git a/src/PJ_aitoff.c b/src/PJ_aitoff.c
index d5350f2c..6c766e5a 100644
--- a/src/PJ_aitoff.c
+++ b/src/PJ_aitoff.c
@@ -2,7 +2,9 @@
  * Project:  PROJ.4
  * Purpose:  Implementation of the aitoff (Aitoff) and wintri (Winkel Tripel)
  *           projections.
- * Author:   Gerald Evenden
+ * Author:   Gerald Evenden (1995)
+ *           Drazen Tutic, Lovro Gradiser (2015) - add inverse
+ *           Thomas Knudsen (2016) - revise/add regression tests
  *
  ******************************************************************************
  * Copyright (c) 1995, Gerald Evenden
@@ -26,12 +28,15 @@
  * DEALINGS IN THE SOFTWARE.
  *****************************************************************************/
 
-#define PROJ_PARMS__ \
-	double	cosphi1; \
-	int		mode;
 #define PJ_LIB__
 #include <projects.h>
 
+
+struct pj_opaque {
+	double	cosphi1;
+	int		mode;
+};
+
 #ifndef M_PI
 #  define M_PI 3.14159265358979323846
 #endif
@@ -42,7 +47,16 @@
 PROJ_HEAD(aitoff, "Aitoff") "\n\tMisc Sph";
 PROJ_HEAD(wintri, "Winkel Tripel") "\n\tMisc Sph\n\tlat_1";
 
+
+
+#if 0
 FORWARD(s_forward); /* spheroid */
+#endif
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double c, d;
 
 	if((d = acos(cos(lp.phi) * cos(c = 0.5 * lp.lam)))) {/* basic Aitoff */
@@ -50,8 +64,8 @@ FORWARD(s_forward); /* spheroid */
 		xy.y *= d * sin(lp.phi);
 	} else
 		xy.x = xy.y = 0.;
-	if (P->mode) { /* Winkel Tripel */
-		xy.x = (xy.x + lp.lam * P->cosphi1) * 0.5;
+	if (Q->mode) { /* Winkel Tripel */
+		xy.x = (xy.x + lp.lam * Q->cosphi1) * 0.5;
 		xy.y = (xy.y + lp.phi) * 0.5;
 	}
 	return (xy);
@@ -60,13 +74,13 @@ FORWARD(s_forward); /* spheroid */
 /***********************************************************************************
 *
 * Inverse functions added by Drazen Tutic and Lovro Gradiser based on paper:
-*           
-* I.Özbug Biklirici and Cengizhan Ipbüker. A General Algorithm for the Inverse 
+*
+* I.Özbug Biklirici and Cengizhan Ipbüker. A General Algorithm for the Inverse
 * Transformation of Map Projections Using Jacobian Matrices. In Proceedings of the
 * Third International Symposium Mathematical & Computational Applications,
 * pages 175{182, Turkey, September 2002.
 *
-* Expected accuracy is defined by EPSILON = 1e-12. Should be appropriate for 
+* Expected accuracy is defined by EPSILON = 1e-12. Should be appropriate for
 * most applications of Aitoff and Winkel Tripel projections.
 *
 * Longitudes of 180W and 180E can be mixed in solution obtained.
@@ -78,19 +92,21 @@ FORWARD(s_forward); /* spheroid */
 *
 ************************************************************************************/
 
-INVERSE(s_inverse); /* sphere */  
-        int iter, MAXITER = 10, round = 0, MAXROUND = 20;
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    int iter, MAXITER = 10, round = 0, MAXROUND = 20;
 	double EPSILON = 1e-12, D, C, f1, f2, f1p, f1l, f2p, f2l, dp, dl, sl, sp, cp, cl, x, y;
 
-	if ((fabs(xy.x) < EPSILON) && (fabs(xy.y) < EPSILON )) { lp.phi = 0.; lp.lam = 0.; return (lp); }
+	if ((fabs(xy.x) < EPSILON) && (fabs(xy.y) < EPSILON )) { lp.phi = 0.; lp.lam = 0.; return lp; }
 
 	/* intial values for Newton-Raphson method */
-	lp.phi = xy.y; lp.lam = xy.x; 
+	lp.phi = xy.y; lp.lam = xy.x;
 	do {
-		iter = 0; 
-		do { 
-			sl = sin(lp.lam * 0.5); cl = cos(lp.lam * 0.5); 
-			sp = sin(lp.phi); cp = cos(lp.phi); 
+		iter = 0;
+		do {
+			sl = sin(lp.lam * 0.5); cl = cos(lp.lam * 0.5);
+			sp = sin(lp.phi); cp = cos(lp.phi);
 			D = cp * cl;
  		      	C = 1. - D * D;
 			D = acos(D) / pow(C, 1.5);
@@ -100,11 +116,11 @@ INVERSE(s_inverse); /* sphere */
 		       	f1l = cp * cp * sl * sl / C + D * cp * cl * sp * sp;
 		    	f2p = sp * sp * cl / C + D * sl * sl * cp;
 		      	f2l = 0.5 * (sp * cp * sl / C - D * sp * cp * cp * sl * cl);
-		      	if (P->mode) { /* Winkel Tripel */
-				f1 = 0.5 * (f1 + lp.lam * P->cosphi1);
+		      	if (Q->mode) { /* Winkel Tripel */
+				f1 = 0.5 * (f1 + lp.lam * Q->cosphi1);
 				f2 = 0.5 * (f2 + lp.phi);
 				f1p *= 0.5;
-				f1l = 0.5 * (f1l + P->cosphi1);
+				f1l = 0.5 * (f1l + Q->cosphi1);
 				f2p = 0.5 * (f2p + 1.);
 				f2l *= 0.5;
 			}
@@ -115,9 +131,9 @@ INVERSE(s_inverse); /* sphere */
 			while (dl < -M_PI) dl += M_PI; /* set to interval [-M_PI, M_PI]  */
 			lp.phi -= dp;	lp.lam -= dl;
 		} while ((fabs(dp) > EPSILON || fabs(dl) > EPSILON) && (iter++ < MAXITER));
-		if (lp.phi > M_PI_2) lp.phi -= 2.*(lp.phi-M_PI_2); /* correct if symmetrical solution for Aitoff */ 
-		if (lp.phi < -M_PI_2) lp.phi -= 2.*(lp.phi+M_PI_2); /* correct if symmetrical solution for Aitoff */ 
-		if ((fabs(fabs(lp.phi) - M_PI_2) < EPSILON) && (!P->mode)) lp.lam = 0.; /* if pole in Aitoff, return longitude of 0 */ 
+		if (lp.phi > M_PI_2) lp.phi -= 2.*(lp.phi-M_PI_2); /* correct if symmetrical solution for Aitoff */
+		if (lp.phi < -M_PI_2) lp.phi -= 2.*(lp.phi+M_PI_2); /* correct if symmetrical solution for Aitoff */
+		if ((fabs(fabs(lp.phi) - M_PI_2) < EPSILON) && (!Q->mode)) lp.lam = 0.; /* if pole in Aitoff, return longitude of 0 */
 
 		/* calculate x,y coordinates with solution obtained */
 		if((D = acos(cos(lp.phi) * cos(C = 0.5 * lp.lam)))) {/* Aitoff */
@@ -125,8 +141,8 @@ INVERSE(s_inverse); /* sphere */
 			y *= D * sin(lp.phi);
 		} else
 			x = y = 0.;
-		if (P->mode) { /* Winkel Tripel */
-			x = (x + lp.lam * P->cosphi1) * 0.5;
+		if (Q->mode) { /* Winkel Tripel */
+			x = (x + lp.lam * Q->cosphi1) * 0.5;
 			y = (y + lp.phi) * 0.5;
 		}
 	/* if too far from given values of x,y, repeat with better approximation of phi,lam */
@@ -134,27 +150,154 @@ INVERSE(s_inverse); /* sphere */
 
 	if (iter == MAXITER && round == MAXROUND) fprintf(stderr, "Warning: Accuracy of 1e-12 not reached. Last increments: dlat=%e and dlon=%e\n", dp, dl);
 
-	return (lp);
+	return lp;
 }
 
-FREEUP; if (P) pj_dalloc(P); }
-	static PJ *
-setup(PJ *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;
+}
+
+static PJ *setup(PJ *P) {
 	P->inv = s_inverse;
 	P->fwd = s_forward;
 	P->es = 0.;
 	return P;
 }
-ENTRY0(aitoff)
-	P->mode = 0;
-ENDENTRY(setup(P))
-ENTRY0(wintri)
-	P->mode = 1;
-	if (pj_param(P->ctx, P->params, "tlat_1").i)
-        {
-		if ((P->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f)) == 0.)
+
+
+PJ *PROJECTION(aitoff) {
+    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_aitoff;
+	Q->mode = 0;
+    return setup(P);
+}
+
+
+PJ *PROJECTION(wintri) {
+    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_wintri;
+	Q->mode = 1;
+	if (pj_param(P->ctx, P->params, "tlat_1").i) {
+		if ((Q->cosphi1 = cos(pj_param(P->ctx, P->params, "rlat_1").f)) == 0.)
 			E_ERROR(-22)
-        }
+    }
 	else /* 50d28' or acos(2/pi) */
-		P->cosphi1 = 0.636619772367581343;
-ENDENTRY(setup(P))
+		Q->cosphi1 = 0.636619772367581343;
+    return setup(P);
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_aitoff_selftest (void) {return 0;}
+#else
+
+int pj_aitoff_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=aitoff   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+
+    XY s_fwd_expect[] = {
+        {223379.45881169615,  111706.74288385305},
+        {223379.45881169615,  -111706.74288385305},
+        {-223379.45881169615,  111706.74288385305},
+        {-223379.45881169615,  -111706.74288385305},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+
+    LP s_inv_expect[] = {
+        {0.0017904931100388164,  0.00089524655491012516},
+        {0.0017904931100388164,  -0.00089524655491012516},
+        {-0.0017904931100388164,  0.00089524655491012516},
+        {-0.0017904931100388164,  -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
+
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_wintri_selftest (void) {return 0;}
+#else
+
+int pj_wintri_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=wintri   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        {223390.80153348515,  111703.90750574505},
+        {223390.80153348515,  -111703.90750574505},
+        {-223390.80153348515,  111703.90750574505},
+        {-223390.80153348515,  -111703.90750574505},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        {0.0017904931099113196,  0.00089524655490101819},
+        {0.0017904931099113196,  -0.00089524655490101819},
+        {-0.0017904931099113196,  0.00089524655490101819},
+        {-0.0017904931099113196,  -0.00089524655490101819},
+    };
+
+    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_august.c b/src/PJ_august.c
index 496fc3c3..4d4ecf6e 100644
--- a/src/PJ_august.c
+++ b/src/PJ_august.c
@@ -1,8 +1,13 @@
 #define PJ_LIB__
 #include	<projects.h>
+
+
 PROJ_HEAD(august, "August Epicycloidal") "\n\tMisc Sph, no inv.";
 #define M 1.333333333333333
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
 	double t, c1, c, x1, x12, y1, y12;
 	(void) P;
 
@@ -15,5 +20,54 @@ FORWARD(s_forward); /* spheroid */
 	xy.y = M * y1 * (3. + 3. * x12 - y12);
 	return (xy);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(august) P->inv = 0; P->fwd = s_forward; P->es = 0.; ENDENTRY(P)
+
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(august) {
+    P->inv = 0;
+    P->fwd = s_forward;
+    P->es = 0.;
+    P->pfree = freeup;
+    P->descr = des_august;
+    return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_august_selftest (void) {return 0;}
+#else
+
+int pj_august_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=august   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        {223404.97818097242,  111722.34028976287},
+        {223404.97818097242,  -111722.34028976287},
+        {-223404.97818097242,  111722.34028976287},
+        {-223404.97818097242,  -111722.34028976287},
+    };
+
+    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_bacon.c b/src/PJ_bacon.c
index dfc559b8..b71fae28 100644
--- a/src/PJ_bacon.c
+++ b/src/PJ_bacon.c
@@ -1,19 +1,27 @@
 # define HLFPI2	2.46740110027233965467
 # define EPS	1e-10
-#define PROJ_PARMS__ \
-	int bacn; \
-	int ortl;
 #define PJ_LIB__
 #include	<projects.h>
+
+
+struct pj_opaque {
+	int bacn;
+	int ortl;
+};
+
 PROJ_HEAD(apian, "Apian Globular I") "\n\tMisc Sph, no inv.";
 PROJ_HEAD(ortel, "Ortelius Oval") "\n\tMisc Sph, no inv.";
 PROJ_HEAD(bacon, "Bacon Globular") "\n\tMisc Sph, no inv.";
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double ax, f;
 
-	xy.y = P->bacn ? HALFPI * sin(lp.phi) : lp.phi;
+	xy.y = Q->bacn ? HALFPI * sin(lp.phi) : lp.phi;
 	if ((ax = fabs(lp.lam)) >= EPS) {
-		if (P->ortl && ax >= HALFPI)
+		if (Q->ortl && ax >= HALFPI)
 			xy.x = sqrt(HLFPI2 - lp.phi * lp.phi + EPS) + ax - HALFPI;
 		else {
 			f = 0.5 * (HLFPI2 / ax + ax);
@@ -24,18 +32,154 @@ FORWARD(s_forward); /* spheroid */
 		xy.x = 0.;
 	return (xy);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(bacon)
-	P->bacn = 1;
-	P->ortl = 0;
-	P->es = 0.; P->fwd = s_forward;
-ENDENTRY(P)
-ENTRY0(apian)
-	P->bacn = P->ortl = 0;
-	P->es = 0.; P->fwd = s_forward;
-ENDENTRY(P)
-ENTRY0(ortel)
-	P->bacn = 0;
-	P->ortl = 1;
-	P->es = 0.; P->fwd = s_forward;
-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(bacon) {
+    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_bacon;
+	Q->bacn = 1;
+	Q->ortl = 0;
+	P->es = 0.;
+    P->fwd = s_forward;
+    return P;
+}
+
+
+PJ *PROJECTION(apian) {
+    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_apian;
+	Q->bacn = Q->ortl = 0;
+	P->es = 0.;
+    P->fwd = s_forward;
+    return P;
+}
+
+
+PJ *PROJECTION(ortel) {
+    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_ortel;
+	Q->bacn = 0;
+	Q->ortl = 1;
+	P->es = 0.;
+    P->fwd = s_forward;
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_bacon_selftest (void) {return 0;}
+#else
+int pj_bacon_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=bacon   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        {223334.13255596498,  175450.72592266591},
+        {223334.13255596498,  -175450.72592266591},
+        {-223334.13255596498,  175450.72592266591},
+        {-223334.13255596498,  -175450.72592266591},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 0, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+#endif
+
+
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_apian_selftest (void) {return 0;}
+#else
+int pj_apian_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=apian   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223374.57735525275,   111701.07212763709},
+        { 223374.57735525275,  -111701.07212763709},
+        {-223374.57735525275,   111701.07212763709},
+        {-223374.57735525275,  -111701.07212763709},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 0, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+#endif
+
+
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_ortel_selftest (void) {return 0;}
+#else
+int pj_ortel_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=ortel   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223374.57735525275,   111701.07212763709},
+        { 223374.57735525275,  -111701.07212763709},
+        {-223374.57735525275,   111701.07212763709},
+        {-223374.57735525275,  -111701.07212763709},
+    };
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 0, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+#endif
+
diff --git a/src/PJ_bipc.c b/src/PJ_bipc.c
index 471e7266..5616b817 100644
--- a/src/PJ_bipc.c
+++ b/src/PJ_bipc.c
@@ -1,132 +1,226 @@
-#define PROJ_PARMS__ \
-	int	noskew;
 #define PJ_LIB__
-# include	<projects.h>
-PROJ_HEAD(bipc, "Bipolar conic of western hemisphere")
-	"\n\tConic Sph.";
-# define EPS	1e-10
-# define EPS10	1e-10
+#include <projects.h>
+
+PROJ_HEAD(bipc, "Bipolar conic of western hemisphere") "\n\tConic Sph.";
+
+# define EPS    1e-10
+# define EPS10  1e-10
 # define ONEEPS 1.000000001
-# define NITER	10
-# define lamB	-.34894976726250681539
-# define n	.63055844881274687180
-# define F	1.89724742567461030582
-# define Azab	.81650043674686363166
-# define Azba	1.82261843856185925133
-# define T	1.27246578267089012270
-# define rhoc	1.20709121521568721927
-# define cAzc	.69691523038678375519
-# define sAzc	.71715351331143607555
-# define C45	.70710678118654752469
-# define S45	.70710678118654752410
-# define C20	.93969262078590838411
-# define S20	-.34202014332566873287
-# define R110	1.91986217719376253360
-# define R104	1.81514242207410275904
-FORWARD(s_forward); /* spheroid */
-	double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
-	int tag;
-
-	cphi = cos(lp.phi);
-	sphi = sin(lp.phi);
-	cdlam = cos(sdlam = lamB - lp.lam);
-	sdlam = sin(sdlam);
-	if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
-		Az = lp.phi < 0. ? PI : 0.;
-		tphi = HUGE_VAL;
-	} else {
-		tphi = sphi / cphi;
-		Az = atan2(sdlam , C45 * (tphi - cdlam));
-	}
-	if( (tag = (Az > Azba)) ) {
-		cdlam = cos(sdlam = lp.lam + R110);
-		sdlam = sin(sdlam);
-		z = S20 * sphi + C20 * cphi * cdlam;
-		if (fabs(z) > 1.) {
-			if (fabs(z) > ONEEPS) F_ERROR
-			else z = z < 0. ? -1. : 1.;
-		} else
-			z = acos(z);
-		if (tphi != HUGE_VAL)
-			Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
-		Av = Azab;
-		xy.y = rhoc;
-	} else {
-		z = S45 * (sphi + cphi * cdlam);
-		if (fabs(z) > 1.) {
-			if (fabs(z) > ONEEPS) F_ERROR
-			else z = z < 0. ? -1. : 1.;
-		} else
-			z = acos(z);
-		Av = Azba;
-		xy.y = -rhoc;
-	}
-	if (z < 0.) F_ERROR;
-	r = F * (t = pow(tan(.5 * z), n));
-	if ((al = .5 * (R104 - z)) < 0.) F_ERROR;
-	al = (t + pow(al, n)) / T;
-	if (fabs(al) > 1.) {
-		if (fabs(al) > ONEEPS) F_ERROR
-		else al = al < 0. ? -1. : 1.;
-	} else
-		al = acos(al);
-	if (fabs(t = n * (Av - Az)) < al)
-		r /= cos(al + (tag ? t : -t));
-	xy.x = r * sin(t);
-	xy.y += (tag ? -r : r) * cos(t);
-	if (P->noskew) {
-		t = xy.x;
-		xy.x = -xy.x * cAzc - xy.y * sAzc; 
-		xy.y = -xy.y * cAzc + t * sAzc; 
-	}
-	return (xy);
+# define NITER  10
+# define lamB   -.34894976726250681539
+# define n  .63055844881274687180
+# define F  1.89724742567461030582
+# define Azab   .81650043674686363166
+# define Azba   1.82261843856185925133
+# define T  1.27246578267089012270
+# define rhoc   1.20709121521568721927
+# define cAzc   .69691523038678375519
+# define sAzc   .71715351331143607555
+# define C45    .70710678118654752469
+# define S45    .70710678118654752410
+# define C20    .93969262078590838411
+# define S20    -.34202014332566873287
+# define R110   1.91986217719376253360
+# define R104   1.81514242207410275904
+
+
+struct pj_opaque {
+    int noskew;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double cphi, sphi, tphi, t, al, Az, z, Av, cdlam, sdlam, r;
+    int tag;
+
+    cphi = cos(lp.phi);
+    sphi = sin(lp.phi);
+    cdlam = cos(sdlam = lamB - lp.lam);
+    sdlam = sin(sdlam);
+    if (fabs(fabs(lp.phi) - HALFPI) < EPS10) {
+        Az = lp.phi < 0. ? PI : 0.;
+        tphi = HUGE_VAL;
+    } else {
+        tphi = sphi / cphi;
+        Az = atan2(sdlam , C45 * (tphi - cdlam));
+    }
+    if( (tag = (Az > Azba)) ) {
+        cdlam = cos(sdlam = lp.lam + R110);
+        sdlam = sin(sdlam);
+        z = S20 * sphi + C20 * cphi * cdlam;
+        if (fabs(z) > 1.) {
+            if (fabs(z) > ONEEPS) F_ERROR
+            else z = z < 0. ? -1. : 1.;
+        } else
+            z = acos(z);
+        if (tphi != HUGE_VAL)
+            Az = atan2(sdlam, (C20 * tphi - S20 * cdlam));
+        Av = Azab;
+        xy.y = rhoc;
+    } else {
+        z = S45 * (sphi + cphi * cdlam);
+        if (fabs(z) > 1.) {
+            if (fabs(z) > ONEEPS) F_ERROR
+            else z = z < 0. ? -1. : 1.;
+        } else
+            z = acos(z);
+        Av = Azba;
+        xy.y = -rhoc;
+    }
+    if (z < 0.) F_ERROR;
+    r = F * (t = pow(tan(.5 * z), n));
+    if ((al = .5 * (R104 - z)) < 0.) F_ERROR;
+    al = (t + pow(al, n)) / T;
+    if (fabs(al) > 1.) {
+        if (fabs(al) > ONEEPS) F_ERROR
+        else al = al < 0. ? -1. : 1.;
+    } else
+        al = acos(al);
+    if (fabs(t = n * (Av - Az)) < al)
+        r /= cos(al + (tag ? t : -t));
+    xy.x = r * sin(t);
+    xy.y += (tag ? -r : r) * cos(t);
+    if (Q->noskew) {
+        t = xy.x;
+        xy.x = -xy.x * cAzc - xy.y * sAzc;
+        xy.y = -xy.y * cAzc + t * sAzc;
+    }
+    return (xy);
 }
-INVERSE(s_inverse); /* spheroid */
-	double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
-	int neg, i;
-
-	if (P->noskew) {
-		t = xy.x;
-		xy.x = -xy.x * cAzc + xy.y * sAzc; 
-		xy.y = -xy.y * cAzc - t * sAzc; 
-	}
-	if( (neg = (xy.x < 0.)) ) {
-		xy.y = rhoc - xy.y;
-		s = S20;
-		c = C20;
-		Av = Azab;
-	} else {
-		xy.y += rhoc;
-		s = S45;
-		c = C45;
-		Av = Azba;
-	}
-	rl = rp = r = hypot(xy.x, xy.y);
-	fAz = fabs(Az = atan2(xy.x, xy.y));
-	for (i = NITER; i ; --i) {
-		z = 2. * atan(pow(r / F,1 / n));
-		al = acos((pow(tan(.5 * z), n) +
-		   pow(tan(.5 * (R104 - z)), n)) / T);
-		if (fAz < al)
-			r = rp * cos(al + (neg ? Az : -Az));
-		if (fabs(rl - r) < EPS)
-			break;
-		rl = r;
-	}
-	if (! i) I_ERROR;
-	Az = Av - Az / n;
-	lp.phi = asin(s * cos(z) + c * sin(z) * cos(Az));
-	lp.lam = atan2(sin(Az), c / tan(z) - s * cos(Az));
-	if (neg)
-		lp.lam -= R110;
-	else
-		lp.lam = lamB - lp.lam;
-	return (lp);
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double t, r, rp, rl, al, z, fAz, Az, s, c, Av;
+    int neg, i;
+
+    if (Q->noskew) {
+        t = xy.x;
+        xy.x = -xy.x * cAzc + xy.y * sAzc;
+        xy.y = -xy.y * cAzc - t * sAzc;
+    }
+    if( (neg = (xy.x < 0.)) ) {
+        xy.y = rhoc - xy.y;
+        s = S20;
+        c = C20;
+        Av = Azab;
+    } else {
+        xy.y += rhoc;
+        s = S45;
+        c = C45;
+        Av = Azba;
+    }
+    rl = rp = r = hypot(xy.x, xy.y);
+    fAz = fabs(Az = atan2(xy.x, xy.y));
+    for (i = NITER; i ; --i) {
+        z = 2. * atan(pow(r / F,1 / n));
+        al = acos((pow(tan(.5 * z), n) +
+           pow(tan(.5 * (R104 - z)), n)) / T);
+        if (fAz < al)
+            r = rp * cos(al + (neg ? Az : -Az));
+        if (fabs(rl - r) < EPS)
+            break;
+        rl = r;
+    }
+    if (! i) I_ERROR;
+    Az = Av - Az / n;
+    lp.phi = asin(s * cos(z) + c * sin(z) * cos(Az));
+    lp.lam = atan2(sin(Az), c / tan(z) - s * cos(Az));
+    if (neg)
+        lp.lam -= R110;
+    else
+        lp.lam = lamB - lp.lam;
+    return (lp);
+}
+
+
+static void *freeup_new (PJ *P) {                        /* Destructor */
+    if (0==P)
+        return 0;
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(bipc)
-	P->noskew = pj_param(P->ctx, P->params, "bns").i;
-	P->inv = s_inverse;
-	P->fwd = s_forward;
-	P->es = 0.;
-ENDENTRY(P)
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(bipc) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->noskew = pj_param(P->ctx, P->params, "bns").i;
+    P->pfree = freeup;
+    P->descr = des_bipc;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    P->es = 0.;
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_bipc_selftest (void) {return 0;}
+#else
+
+int pj_bipc_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=bipc   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=bipc   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {2452160.2177257561,  -14548450.759654747},
+        {2447915.213725341,  -14763427.21279873},
+        {2021695.5229349085,  -14540413.695283702},
+        {2018090.5030046992,  -14755620.651414108},
+    };
+
+    XY s_fwd_expect[] = {
+        {2460565.7409749646,  -14598319.9893308},
+        {2456306.1859352002,  -14814033.339502094},
+        {2028625.4978190989,  -14590255.375482792},
+        {2025008.1205891429,  -14806200.018759441},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {-73.038700284978702,  17.248118466239116},
+        {-73.03730373933017,  17.249414978178777},
+        {-73.03589317304332,  17.245536403008771},
+        {-73.034496627213585,  17.246832895573739},
+    };
+
+    LP s_inv_expect[] = {
+        {-73.038693104942126,  17.248116270440242},
+        {-73.037301330021322,  17.24940835333777},
+        {-73.035895582251086,  17.245543027866539},
+        {-73.034503807150301,  17.246835091521532},
+    };
+
+    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_boggs.c b/src/PJ_boggs.c
index b22a25ff..dc4d781f 100644
--- a/src/PJ_boggs.c
+++ b/src/PJ_boggs.c
@@ -8,7 +8,10 @@ PROJ_HEAD(boggs, "Boggs Eumorphic") "\n\tPCyl., no inv., Sph.";
 # define FXC2	1.11072
 # define FYC	0.49931
 # define FYC2	1.41421356237309504880
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
 	double theta, th1, c;
 	int i;
 	(void) P;
@@ -29,5 +32,49 @@ FORWARD(s_forward); /* spheroid */
 	xy.y = FYC * (lp.phi + FYC2 * sin(theta));
 	return (xy);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(boggs) P->es = 0.; P->fwd = s_forward; ENDENTRY(P)
+
+
+static void *freeup_new (PJ *P) {                        /* Destructor */
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(boggs) {
+    P->es = 0.;
+    P->fwd = s_forward;
+    P->pfree = freeup;
+    P->descr = des_boggs;
+    return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_boggs_selftest (void) {return 0;}
+#else
+int pj_boggs_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=boggs   +a=6400000    +lat_1=0 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 211949.70080818201,   117720.99830541089},
+        { 211949.70080818201,  -117720.99830541089},
+        {-211949.70080818201,   117720.99830541089},
+        {-211949.70080818201,  -117720.99830541089},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, tolerance_lp, 4, 0, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+#endif
diff --git a/src/PJ_bonne.c b/src/PJ_bonne.c
index 4b1e7866..cb17104a 100644
--- a/src/PJ_bonne.c
+++ b/src/PJ_bonne.c
@@ -1,51 +1,70 @@
-#define PROJ_PARMS__ \
-	double phi1; \
-	double cphi1; \
-	double am1; \
-	double m1; \
-	double *en;
 #define PJ_LIB__
 #include	<projects.h>
+
 PROJ_HEAD(bonne, "Bonne (Werner lat_1=90)")
 	"\n\tConic Sph&Ell\n\tlat_1=";
 #define EPS10	1e-10
-FORWARD(e_forward); /* ellipsoid */
+
+struct pj_opaque {
+	double phi1;
+	double cphi1;
+	double am1;
+	double m1;
+	double *en;
+};
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double rh, E, c;
 
-	rh = P->am1 + P->m1 - pj_mlfn(lp.phi, E = sin(lp.phi), c = cos(lp.phi), P->en);
+	rh = Q->am1 + Q->m1 - pj_mlfn(lp.phi, E = sin(lp.phi), c = cos(lp.phi), Q->en);
 	E = c * lp.lam / (rh * sqrt(1. - P->es * E * E));
 	xy.x = rh * sin(E);
-	xy.y = P->am1 - rh * cos(E);
-	return (xy);
+	xy.y = Q->am1 - rh * cos(E);
+	return xy;
 }
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double E, rh;
 
-	rh = P->cphi1 + P->phi1 - lp.phi;
+	rh = Q->cphi1 + Q->phi1 - lp.phi;
 	if (fabs(rh) > EPS10) {
 		xy.x = rh * sin(E = lp.lam * cos(lp.phi) / rh);
-		xy.y = P->cphi1 - rh * cos(E);
+		xy.y = Q->cphi1 - rh * cos(E);
 	} else
 		xy.x = xy.y = 0.;
-	return (xy);
+	return xy;
 }
-INVERSE(s_inverse); /* spheroid */
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double rh;
 
-	rh = hypot(xy.x, xy.y = P->cphi1 - xy.y);
-	lp.phi = P->cphi1 + P->phi1 - rh;
+	rh = hypot(xy.x, xy.y = Q->cphi1 - xy.y);
+	lp.phi = Q->cphi1 + Q->phi1 - rh;
 	if (fabs(lp.phi) > HALFPI) I_ERROR;
 	if (fabs(fabs(lp.phi) - HALFPI) <= EPS10)
 		lp.lam = 0.;
 	else
 		lp.lam = rh * atan2(xy.x, xy.y) / cos(lp.phi);
-	return (lp);
+	return lp;
 }
-INVERSE(e_inverse); /* ellipsoid */
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
 	double s, rh;
 
-	rh = hypot(xy.x, xy.y = P->am1 - xy.y);
-	lp.phi = pj_inv_mlfn(P->ctx, P->am1 + P->m1 - rh, P->es, P->en);
+	rh = hypot(xy.x, xy.y = Q->am1 - xy.y);
+	lp.phi = pj_inv_mlfn(P->ctx, Q->am1 + Q->m1 - rh, P->es, Q->en);
 	if ((s = fabs(lp.phi)) < HALFPI) {
 		s = sin(lp.phi);
 		lp.lam = rh * atan2(xy.x, xy.y) *
@@ -53,33 +72,112 @@ INVERSE(e_inverse); /* ellipsoid */
 	} else if (fabs(s - HALFPI) <= EPS10)
 		lp.lam = 0.;
 	else I_ERROR;
-	return (lp);
+	return lp;
 }
-FREEUP;
-	if (P) {
-		if (P->en)
-			pj_dalloc(P->en);
-		pj_dalloc(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->en);
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-ENTRY1(bonne, en)
+
+
+PJ *PROJECTION(bonne) {
 	double c;
+    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_bonne;
 
-	P->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
-	if (fabs(P->phi1) < EPS10) E_ERROR(-23);
+	Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+	if (fabs(Q->phi1) < EPS10) E_ERROR(-23);
 	if (P->es) {
-		P->en = pj_enfn(P->es);
-		P->m1 = pj_mlfn(P->phi1, P->am1 = sin(P->phi1),
-			c = cos(P->phi1), P->en);
-		P->am1 = c / (sqrt(1. - P->es * P->am1 * P->am1) * P->am1);
+		Q->en = pj_enfn(P->es);
+		Q->m1 = pj_mlfn(Q->phi1, Q->am1 = sin(Q->phi1),
+			c = cos(Q->phi1), Q->en);
+		Q->am1 = c / (sqrt(1. - P->es * Q->am1 * Q->am1) * Q->am1);
 		P->inv = e_inverse;
 		P->fwd = e_forward;
 	} else {
-		if (fabs(P->phi1) + EPS10 >= HALFPI)
-			P->cphi1 = 0.;
+		if (fabs(Q->phi1) + EPS10 >= HALFPI)
+			Q->cphi1 = 0.;
 		else
-			P->cphi1 = 1. / tan(P->phi1);
+			Q->cphi1 = 1. / tan(Q->phi1);
 		P->inv = s_inverse;
 		P->fwd = s_forward;
 	}
-ENDENTRY(P)
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_bonne_selftest (void) {return 0;}
+#else
+int pj_bonne_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=bonne   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=bonne   +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.29609715697,   55321.139565494814},
+        { 222605.29609923941,  -165827.64779905154},
+        {-222605.29609715697,   55321.139565494814},
+        {-222605.29609923941,  -165827.64779905154},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223368.11557252839,   55884.555246393575},
+        { 223368.11557463196,  -167517.59936969393},
+        {-223368.11557252839,   55884.555246393575},
+        {-223368.11557463196,  -167517.59936969393},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017966987691132891,  0.50090436853737497},
+        { 0.0017966982774478867,  0.4990956309655612},
+        {-0.0017966987691132891,  0.50090436853737497},
+        {-0.0017966982774478867,  0.4990956309655612},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017905615332457991,  0.50089524631087834},
+        { 0.0017905610449335603,  0.49910475320072978},
+        {-0.0017905615332457991,  0.50089524631087834},
+        {-0.0017905610449335603,  0.49910475320072978},
+    };
+
+    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_run_selftests.c b/src/pj_run_selftests.c
index 44573b2b..45aa7203 100644
--- a/src/pj_run_selftests.c
+++ b/src/pj_run_selftests.c
@@ -33,6 +33,16 @@
 
 extern int pj_aea_selftest(void);
 extern int pj_leac_selftest(void);
+extern int pj_airy_selftest(void);
+extern int pj_aitoff_selftest(void);
+extern int pj_wintri_selftest(void);
+extern int pj_august_selftest(void);
+extern int pj_bacon_selftest(void);
+extern int pj_apian_selftest(void);
+extern int pj_ortel_selftest(void);
+extern int pj_boggs_selftest(void);
+extern int pj_bonne_selftest(void);
+extern int pj_bipc_selftest(void);
 
 
 
@@ -48,16 +58,33 @@ static void run_one_test (const char *mnemonic, int (testfunc)(void), int verbos
 }
 
 
+#ifdef PJ_OMIT_SELFTEST
+int pj_run_selftests (int verbosity) {
+    printf ("This version of libproj is not configured for internal regression tests.\n");
+    return 0;
+}
+#else
 int pj_run_selftests (int verbosity) {
     int n_ok = 0, n_ko = 0;
 
     if (verbosity)
        printf ("Running internal regression tests\n");
-    run_one_test ("aea",    pj_aea_selftest,    verbosity, &n_ok, &n_ko);
-    run_one_test ("leac",   pj_leac_selftest,   verbosity, &n_ok, &n_ko);
+    run_one_test ("aea",      pj_aea_selftest,       verbosity, &n_ok, &n_ko);
+    run_one_test ("leac",     pj_leac_selftest,      verbosity, &n_ok, &n_ko);
+    run_one_test ("airy",     pj_airy_selftest,      verbosity, &n_ok, &n_ko);
+    run_one_test ("aitoff",   pj_aitoff_selftest,    verbosity, &n_ok, &n_ko);
+    run_one_test ("wintri",   pj_wintri_selftest,    verbosity, &n_ok, &n_ko);
+    run_one_test ("august",   pj_august_selftest,    verbosity, &n_ok, &n_ko);
+    run_one_test ("bacon",    pj_bacon_selftest,     verbosity, &n_ok, &n_ko);
+    run_one_test ("apian",    pj_apian_selftest,     verbosity, &n_ok, &n_ko);
+    run_one_test ("ortel",    pj_ortel_selftest,     verbosity, &n_ok, &n_ko);
+    run_one_test ("boggs",    pj_boggs_selftest,     verbosity, &n_ok, &n_ko);
+    run_one_test ("bonne",    pj_bonne_selftest,     verbosity, &n_ok, &n_ko);
+    run_one_test ("bipc",     pj_bipc_selftest,     verbosity, &n_ok, &n_ko);
 
     if (0==verbosity)
        printf ("Internal regression tests done. ");
     printf ("Total: %d, Failure: %d, Success: %d\n", n_ok+n_ko, n_ko, n_ok);
     return n_ko;
 }
+#endif