Merge pull request #4 from kbevers/fix-projs-with-d-e-f

Projections starting with d, e and f converted to new style

13 files changed, 1606 insertions, 531 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index 8cca2386..4db5e386 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -355,25 +355,12 @@ source files
 
 int pj_aeqd_selftest     (void)     {return 10000;}
 int pj_alsk_selftest     (void)     {return 10000;}
-int pj_denoy_selftest    (void)     {return 10000;}
-int pj_eck1_selftest     (void)     {return 10000;}
-int pj_eck2_selftest     (void)     {return 10000;}
-int pj_eck3_selftest     (void)     {return 10000;}
-int pj_eck4_selftest     (void)     {return 10000;}
-int pj_eck5_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_fahey_selftest (void) {return 10000;}
-
-int pj_fouc_s_selftest (void) {return 10000;}
-int pj_gall_selftest (void) {return 10000;}
-int pj_geos_selftest (void) {return 10000;}
 int pj_geocent_selftest (void) {return 10000;}
-int pj_gins8_selftest (void) {return 10000;}
 int pj_gn_sinu_selftest (void) {return 10000;}
-int pj_gnom_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;}
@@ -385,7 +372,6 @@ int pj_igh_selftest (void) {return 10000;}
 int pj_imw_p_selftest (void) {return 10000;}
 int pj_isea_selftest (void) {return 10000;}
 
-int pj_kav7_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;}
@@ -426,7 +412,6 @@ int pj_omerc_selftest (void) {return 10000;}
 int pj_ortho_selftest (void) {return 10000;}
 int pj_patterson_selftest (void) {return 10000;}
 int pj_poly_selftest (void) {return 10000;}
-int pj_putp1_selftest (void) {return 10000;}
 int pj_putp2_selftest (void) {return 10000;}
 int pj_putp3_selftest (void) {return 10000;}
 int pj_putp3p_selftest (void) {return 10000;}
@@ -453,6 +438,5 @@ int pj_vandg3_selftest (void) {return 10000;}
 int pj_vandg4_selftest (void) {return 10000;}
 int pj_wag4_selftest (void) {return 10000;}
 int pj_wag5_selftest (void) {return 10000;}
-int pj_wag6_selftest (void) {return 10000;}
 int pj_weren_selftest (void) {return 10000;}
 #endif
diff --git a/src/PJ_denoy.c b/src/PJ_denoy.c
index b1a6fe80..10005a31 100644
--- a/src/PJ_denoy.c
+++ b/src/PJ_denoy.c
@@ -1,19 +1,69 @@
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
+
 PROJ_HEAD(denoy, "Denoyer Semi-Elliptical") "\n\tPCyl., no inv., Sph.";
-#define C0	0.95
-#define C1	-.08333333333333333333
-#define C3	.00166666666666666666
-#define D1	0.9
-#define D5	0.03
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	xy.y = lp.phi;
-	xy.x = lp.lam;
-	lp.lam = fabs(lp.lam);
-	xy.x *= cos((C0 + lp.lam * (C1 + lp.lam * lp.lam * C3)) *
-		(lp.phi * (D1 + D5 * lp.phi * lp.phi * lp.phi * lp.phi)));
-	return (xy);
+
+#define C0  0.95
+#define C1 -0.08333333333333333333
+#define C3  0.00166666666666666666
+#define D1  0.9
+#define D5  0.03
+
+
+static XY s_forward (LP lp, PJ *P) {            /* Spheroidal, forward */
+    XY xy = {0.0, 0.0};
+    (void) P;
+    xy.y = lp.phi;
+    xy.x = lp.lam;
+    lp.lam = fabs(lp.lam);
+    xy.x *= cos((C0 + lp.lam * (C1 + lp.lam * lp.lam * C3)) *
+            (lp.phi * (D1 + D5 * lp.phi * lp.phi * lp.phi * lp.phi)));
+    return xy;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(denoy) 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(denoy) {
+    P->es = 0.0;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_denoy_selftest (void) {return 0;}
+#else
+
+int pj_denoy_selftest (void) {
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=denoy   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223377.422876954137,  111701.07212763709},
+        { 223377.422876954137, -111701.07212763709},
+        {-223377.422876954137,  111701.07212763709},
+        {-223377.422876954137, -111701.07212763709},
+    };
+
+    return pj_generic_selftest (0, s_args, tolerance_xy, 0, 4, 4, fwd_in, 0, s_fwd_expect, 0, 0, 0);
+}
+
+
+#endif
diff --git a/src/PJ_eck1.c b/src/PJ_eck1.c
index b0b43da5..da2c8685 100644
--- a/src/PJ_eck1.c
+++ b/src/PJ_eck1.c
@@ -1,21 +1,92 @@
-#define PJ_LIB__
-#include	<projects.h>
-PROJ_HEAD(eck1, "Eckert I") "\n\tPCyl., Sph.";
-#define FC	.92131773192356127802
-#define RP	.31830988618379067154
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	xy.x = FC * lp.lam * (1. - RP * fabs(lp.phi));
-	xy.y = FC * lp.phi;
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	(void) P;
-	lp.phi = xy.y / FC;
-	lp.lam = xy.x / (FC * (1. - RP * fabs(lp.phi)));
-	return (lp);
-}
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(eck1)
-	P->es = 0.; P->inv = s_inverse; P->fwd = s_forward;
-ENDENTRY(P)
+#define PJ_LIB__

+#include <projects.h>

+

+PROJ_HEAD(eck1, "Eckert I") "\n\tPCyl., Sph.";

+#define FC  0.92131773192356127802

+#define RP  0.31830988618379067154

+

+

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

+    XY xy = {0.0,0.0};

+    (void) P;

+

+    xy.x = FC * lp.lam * (1. - RP * fabs(lp.phi));

+    xy.y = FC * lp.phi;

+

+    return xy;

+}

+

+

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

+    LP lp = {0.0,0.0};

+    (void) P;

+

+    lp.phi = xy.y / FC;

+    lp.lam = xy.x / (FC * (1. - RP * fabs(lp.phi)));

+

+    return (lp);

+}

+

+

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

+    return pj_dealloc(P);

+}

+

+

+static void freeup (PJ *P) {

+    freeup_new (P);

+    return;

+}

+

+

+PJ *PROJECTION(eck1) {

+    P->es = 0.0;

+    P->inv = s_inverse;

+    P->fwd = s_forward;

+

+    return P
;

+}

+

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_eck1_selftest (void) {return 0;}

+#else

+

+int pj_eck1_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

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

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+    XY s_fwd_expect[] = {

+        { 204680.88820295094,  102912.17842606473},

+        { 204680.88820295094, -102912.17842606473},

+        {-204680.88820295094,  102912.17842606473},

+        {-204680.88820295094, -102912.17842606473},

+    };

+

+    XY inv_in[] = {

+        { 200, 100},

+        { 200,-100},

+        {-200, 100},

+        {-200,-100}

+    };

+

+    LP s_inv_expect[] = {

+        { 0.0019434150820034624,  0.00097170229538813102},

+        { 0.0019434150820034624, -0.00097170229538813102},

+        {-0.0019434150820034624,  0.00097170229538813102},

+        {-0.0019434150820034624, -0.00097170229538813102},

+    };

+

+    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_eck2.c b/src/PJ_eck2.c
index 08b65595..6d73b88a 100644
--- a/src/PJ_eck2.c
+++ b/src/PJ_eck2.c
@@ -1,29 +1,104 @@
 #define PJ_LIB__
-# include	<projects.h>
+# include   <projects.h>
+
 PROJ_HEAD(eck2, "Eckert II") "\n\tPCyl. Sph.";
-#define FXC	0.46065886596178063902
-#define FYC	1.44720250911653531871
-#define C13	0.33333333333333333333
-#define ONEEPS	1.0000001
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	xy.x = FXC * lp.lam * (xy.y = sqrt(4. - 3. * sin(fabs(lp.phi))));
-	xy.y = FYC * (2. - xy.y);
-	if ( lp.phi < 0.) xy.y = -xy.y;
-	return (xy);
+
+#define FXC     0.46065886596178063902
+#define FYC     1.44720250911653531871
+#define C13     0.33333333333333333333
+#define ONEEPS  1.0000001
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    (void) P;
+
+    xy.x = FXC * lp.lam * (xy.y = sqrt(4. - 3. * sin(fabs(lp.phi))));
+    xy.y = FYC * (2. - xy.y);
+    if ( lp.phi < 0.) xy.y = -xy.y;
+
+    return (xy);
 }
-INVERSE(s_inverse); /* spheroid */
-	lp.lam = xy.x / (FXC * ( lp.phi = 2. - fabs(xy.y) / FYC) );
-	lp.phi = (4. - lp.phi * lp.phi) * C13;
-	if (fabs(lp.phi) >= 1.) {
-		if (fabs(lp.phi) > ONEEPS)	I_ERROR
-		else
-			lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
-	} else
-		lp.phi = asin(lp.phi);
-	if (xy.y < 0)
-		lp.phi = -lp.phi;
-	return (lp);
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    (void) P;
+
+    lp.lam = xy.x / (FXC * ( lp.phi = 2. - fabs(xy.y) / FYC) );
+    lp.phi = (4. - lp.phi * lp.phi) * C13;
+    if (fabs(lp.phi) >= 1.) {
+        if (fabs(lp.phi) > ONEEPS)  I_ERROR
+        else
+            lp.phi = lp.phi < 0. ? -HALFPI : HALFPI;
+    } else
+        lp.phi = asin(lp.phi);
+    if (xy.y < 0)
+        lp.phi = -lp.phi;
+    return (lp);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(eck2); P->es = 0.; P->inv = s_inverse; 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(eck2) {
+    P->es = 0.;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_eck2_selftest (void) {return 0;}
+#else
+
+int pj_eck2_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=eck2   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 204472.87090796008,  121633.73497524235},
+        { 204472.87090796008, -121633.73497524235},
+        {-204472.87090796008,  121633.73497524235},
+        {-204472.87090796008, -121633.73497524235},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0019434150820034624,  0.00082480429919795412},
+        { 0.0019434150820034624, -0.00082480429919795412},
+        {-0.0019434150820034624,  0.00082480429919795412},
+        {-0.0019434150820034624, -0.00082480429919795412},
+    };
+
+    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_eck3.c b/src/PJ_eck3.c
index d7755f0c..3eb7f8f9 100644
--- a/src/PJ_eck3.c
+++ b/src/PJ_eck3.c
@@ -1,50 +1,295 @@
-#define PROJ_PARMS__ \
-	double C_x, C_y, A, B;
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(eck3, "Eckert III") "\n\tPCyl, Sph.";
 PROJ_HEAD(putp1, "Putnins P1") "\n\tPCyl, Sph.";
 PROJ_HEAD(wag6, "Wagner VI") "\n\tPCyl, Sph.";
 PROJ_HEAD(kav7, "Kavraisky VII") "\n\tPCyl, Sph.";
-FORWARD(s_forward); /* spheroid */
-	xy.y = P->C_y * lp.phi;
-	xy.x = P->C_x * lp.lam * (P->A + asqrt(1. - P->B * lp.phi * lp.phi));
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	lp.phi = xy.y / P->C_y;
-	lp.lam = xy.x / (P->C_x * (P->A + asqrt(1. - P->B * lp.phi * lp.phi)));
-	return (lp);
-}
-FREEUP; if (P) pj_dalloc(P); }
-	static PJ *
-setup(PJ *P) {
-	P->es = 0.;
-	P->inv = s_inverse;
-	P->fwd = s_forward;
-	return P;
-}
-ENTRY0(eck3)
-	P->C_x = .42223820031577120149;
-	P->C_y = .84447640063154240298;
-	P->A = 1.;
-	P->B = 0.4052847345693510857755;
-ENDENTRY(setup(P))
-ENTRY0(kav7)
-	P->C_x = 0.2632401569273184856851;
-	P->C_x = 0.8660254037844;
-	P->C_y = 1.;
-	P->A = 0.;
-	P->B = 0.30396355092701331433;
-ENDENTRY(setup(P))
-ENTRY0(wag6);
-	P->C_x = P->C_y = 0.94745;
-	P->A = 0.;
-	P->B = 0.30396355092701331433;
-ENDENTRY(setup(P))
-ENTRY0(putp1);
-	P->C_x = 1.89490;
-	P->C_y = 0.94745;
-	P->A = -0.5;
-	P->B = 0.30396355092701331433;
-ENDENTRY(setup(P))
+
+struct pj_opaque {
+    double C_x, C_y, A, B;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    xy.y = Q->C_y * lp.phi;
+    xy.x = Q->C_x * lp.lam * (Q->A + asqrt(1. - Q->B * lp.phi * 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;
+
+    lp.phi = xy.y / Q->C_y;
+    lp.lam = xy.x / (Q->C_x * (Q->A + asqrt(1. - Q->B * lp.phi * lp.phi)));
+    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;
+}
+
+
+static PJ *setup(PJ *P) {
+    P->es = 0.;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    return P;
+}
+
+
+PJ *PROJECTION(eck3) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->C_x = 0.42223820031577120149;
+    Q->C_y = 0.84447640063154240298;
+    Q->A = 1.0;
+    Q->B = 0.4052847345693510857755;
+
+    return setup(P);
+}
+
+
+PJ *PROJECTION(kav7) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    /* Defined twice in original code - Using 0.866...,
+     * but leaving the other one here as a safety measure.
+     * Q->C_x = 0.2632401569273184856851; */
+    Q->C_x = 0.8660254037844;
+    Q->C_y = 1.;
+    Q->A = 0.;
+    Q->B = 0.30396355092701331433;
+
+    return setup(P);
+}
+
+
+PJ *PROJECTION(wag6) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->C_x = Q->C_y = 0.94745;
+    Q->A = 0.0;
+    Q->B = 0.30396355092701331433;
+
+    return setup(P);
+}
+
+
+PJ *PROJECTION(putp1) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->C_x = 1.89490;
+    Q->C_y = 0.94745;
+    Q->A = -0.5;
+    Q->B = 0.30396355092701331433;
+
+    return setup(P);
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_eck3_selftest (void) {return 0;}
+#else
+
+int pj_eck3_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=eck3   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 188652.01572153764,  94328.919337031271},
+        { 188652.01572153764, -94328.919337031271},
+        {-188652.01572153764,  94328.919337031271},
+        {-188652.01572153764, -94328.919337031271},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0021202405520236059,  0.0010601202759750307},
+        { 0.0021202405520236059, -0.0010601202759750307},
+        {-0.0021202405520236059,  0.0010601202759750307},
+        {-0.0021202405520236059, -0.0010601202759750307},
+    };
+
+    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_kav7_selftest (void) {return 0;}
+#else
+
+int pj_kav7_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=kav7   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 193462.9749437288,  111701.07212763709},
+        { 193462.9749437288, -111701.07212763709},
+        {-193462.9749437288,  111701.07212763709},
+        {-193462.9749437288, -111701.07212763709}
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0020674833579085268,  0.00089524655489191132},
+        { 0.0020674833579085268, -0.00089524655489191132},
+        {-0.0020674833579085268,  0.00089524655489191132},
+        {-0.0020674833579085268, -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
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_wag6_selftest (void) {return 0;}
+#else
+
+int pj_wag6_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=wag6   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 211652.56216440981,  105831.18078732977},
+        { 211652.56216440981, -105831.18078732977},
+        {-211652.56216440981,  105831.18078732977},
+        {-211652.56216440981, -105831.18078732977}
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0018898022163257513,  0.000944901108123818},
+        { 0.0018898022163257513, -0.000944901108123818},
+        {-0.0018898022163257513,  0.000944901108123818},
+        {-0.0018898022163257513, -0.000944901108123818}
+    };
+
+    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_putp1_selftest (void) {return 0;}
+#else
+
+int pj_putp1_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=putp1   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 211642.76275416015,  105831.18078732977},
+        { 211642.76275416015, -105831.18078732977},
+        {-211642.76275416015,  105831.18078732977},
+        {-211642.76275416015, -105831.18078732977}
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0018898022164038663,  0.000944901108123818},
+        { 0.0018898022164038663, -0.000944901108123818},
+        {-0.0018898022164038663,  0.000944901108123818},
+        {-0.0018898022164038663, -0.000944901108123818}
+    };
+
+    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_eck4.c b/src/PJ_eck4.c
index 0bedbbc9..8a56f019 100644
--- a/src/PJ_eck4.c
+++ b/src/PJ_eck4.c
@@ -1,45 +1,117 @@
-#define PJ_LIB__
-#include	<projects.h>
-PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph.";
-#define C_x	.42223820031577120149
-#define C_y	1.32650042817700232218
-#define RC_y	.75386330736002178205
-#define C_p	3.57079632679489661922
-#define RC_p	.28004957675577868795
-#define EPS	1e-7
-#define NITER	6
-FORWARD(s_forward); /* spheroid */
-	double p, V, s, c;
-	int i;
-	(void) P;
-
-	p = C_p * sin(lp.phi);
-	V = lp.phi * lp.phi;
-	lp.phi *= 0.895168 + V * ( 0.0218849 + V * 0.00826809 );
-	for (i = NITER; i ; --i) {
-		c = cos(lp.phi);
-		s = sin(lp.phi);
-		lp.phi -= V = (lp.phi + s * (c + 2.) - p) /
-			(1. + c * (c + 2.) - s * s);
-		if (fabs(V) < EPS)
-			break;
-	}
-	if (!i) {
-		xy.x = C_x * lp.lam;
-		xy.y = lp.phi < 0. ? -C_y : C_y;
-	} else {
-		xy.x = C_x * lp.lam * (1. + cos(lp.phi));
-		xy.y = C_y * sin(lp.phi);
-	}
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	double c;
-
-	lp.phi = aasin(P->ctx,xy.y / C_y);
-	lp.lam = xy.x / (C_x * (1. + (c = cos(lp.phi))));
-	lp.phi = aasin(P->ctx,(lp.phi + sin(lp.phi) * (c + 2.)) / C_p);
-	return (lp);
-}
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(eck4); P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
+#define PJ_LIB__

+#include    <projects.h>

+

+PROJ_HEAD(eck4, "Eckert IV") "\n\tPCyl, Sph.";

+

+#define C_x .42223820031577120149

+#define C_y 1.32650042817700232218

+#define RC_y    .75386330736002178205

+#define C_p 3.57079632679489661922

+#define RC_p    .28004957675577868795

+#define EPS 1e-7

+#define NITER   6

+

+

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

+    XY xy = {0.0,0.0};

+    double p, V, s, c;

+    int i;

+    (void) P;

+

+    p = C_p * sin(lp.phi);

+    V = lp.phi * lp.phi;

+    lp.phi *= 0.895168 + V * ( 0.0218849 + V * 0.00826809 );

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

+        c = cos(lp.phi);

+        s = sin(lp.phi);

+        lp.phi -= V = (lp.phi + s * (c + 2.) - p) /

+            (1. + c * (c + 2.) - s * s);

+        if (fabs(V) < EPS)

+            break;

+    }

+    if (!i) {

+        xy.x = C_x * lp.lam;

+        xy.y = lp.phi < 0. ? -C_y : C_y;

+    } else {

+        xy.x = C_x * lp.lam * (1. + cos(lp.phi));

+        xy.y = C_y * sin(lp.phi);

+    }

+    return xy;

+}

+

+

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

+    LP lp = {0.0,0.0};

+    double c;

+

+    lp.phi = aasin(P->ctx,xy.y / C_y);

+    lp.lam = xy.x / (C_x * (1. + (c = cos(lp.phi))));

+    lp.phi = aasin(P->ctx,(lp.phi + sin(lp.phi) * (c + 2.)) / C_p);

+    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(eck4) {

+    P->es = 0.0;

+    P->inv = s_inverse;

+    P->fwd = s_forward;

+

+    return P;

+}

+

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_eck4_selftest (void) {return 0;}

+#else

+

+int pj_eck4_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

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

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+

+    XY s_fwd_expect[] = {

+        { 188646.38935641639,  132268.54017406539},

+        { 188646.38935641639, -132268.54017406539},

+        {-188646.38935641639,  132268.54017406539},

+        {-188646.38935641639, -132268.54017406539},

+    };

+

+    XY inv_in[] = {

+        { 200, 100},

+        { 200,-100},

+        {-200, 100},

+        {-200,-100}

+    };

+

+    LP s_inv_expect[] = {

+        { 0.0021202405520236059, 0.00075601458836610643},

+        { 0.0021202405520236059, -0.00075601458836610643},

+        {-0.0021202405520236059, 0.00075601458836610643},

+        {-0.0021202405520236059, -0.00075601458836610643},

+    };

+

+    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_eck5.c b/src/PJ_eck5.c
index 029c18e0..ccfcb88f 100644
--- a/src/PJ_eck5.c
+++ b/src/PJ_eck5.c
@@ -1,20 +1,108 @@
-#define PJ_LIB__
-# include	<projects.h>
-PROJ_HEAD(eck5, "Eckert V") "\n\tPCyl, Sph.";
-#define XF	0.44101277172455148219
-#define RXF	2.26750802723822639137
-#define YF	0.88202554344910296438
-#define RYF	1.13375401361911319568
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	xy.x = XF * (1. + cos(lp.phi)) * lp.lam;
-	xy.y = YF * lp.phi;
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	(void) P;
-	lp.lam = RXF * xy.x / (1. + cos( lp.phi = RYF * xy.y));
-	return (lp);
-}
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(eck5); P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
+#define PJ_LIB__

+#include <projects.h>

+

+PROJ_HEAD(eck5, "Eckert V") "\n\tPCyl, Sph.";

+

+#define XF  0.44101277172455148219

+#define RXF 2.26750802723822639137

+#define YF  0.88202554344910296438

+#define RYF 1.13375401361911319568

+

+

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

+    XY xy = {0.0,0.0};

+    (void) P;

+    xy.x = XF * (1. + cos(lp.phi)) * lp.lam;

+    xy.y = YF * lp.phi;

+

+    return xy;

+}

+

+

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

+    LP lp = {0.0,0.0};

+    (void) P;

+    lp.lam = RXF * xy.x / (1. + cos( lp.phi = RYF * xy.y));

+

+    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(eck5) {

+    P->es = 0.0;

+    P->inv = s_inverse;

+    P->fwd = s_forward;

+

+    return P;

+}

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_eck5_selftest (void) {return 0;}

+#else

+

+int pj_eck5_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

+    char e_args[] = {"+proj=eck5   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};

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

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+

+    XY e_fwd_expect[] = {

+        { 196358.31442683787,  98186.634363414545
},

+        { 196358.31442683787, -98186.634363414545
},

+        {-196358.31442683787,  98186.634363414545
},

+        {-196358.31442683787, -98186.634363414545
},

+    };

+

+    XY s_fwd_expect[] = {

+        { 197031.39213406085,  98523.198847226551
},

+        { 197031.39213406085, -98523.198847226551
},

+        {-197031.39213406085,  98523.198847226551
},

+        {-197031.39213406085, -98523.198847226551
},

+    };

+

+    XY inv_in[] = {

+        { 200, 100},

+        { 200,-100},

+        {-200, 100},

+        {-200,-100}

+    };

+

+    LP e_inv_expect[] = {

+        {0.0020369371178927064,  0.0010184685588659013
},

+        {0.0020369371178927064,  -0.0010184685588659013
},

+        {-0.0020369371178927064,  0.0010184685588659013
},

+        {-0.0020369371178927064,  -0.0010184685588659013
},

+    };

+

+    LP s_inv_expect[] = {

+        {0.002029978749734037,  0.001014989374787388
},

+        {0.002029978749734037,  -0.001014989374787388
},

+        {-0.002029978749734037,  0.001014989374787388
},

+        {-0.002029978749734037,  -0.001014989374787388
},

+    };

+

+    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_fahey.c b/src/PJ_fahey.c
index 007fc906..864225f4 100644
--- a/src/PJ_fahey.c
+++ b/src/PJ_fahey.c
@@ -1,19 +1,96 @@
-#define PJ_LIB__
-# include	<projects.h>
-PROJ_HEAD(fahey, "Fahey") "\n\tPcyl, Sph.";
-#define TOL 1e-6
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	xy.y = 1.819152 * ( xy.x = tan(0.5 * lp.phi) );
-	xy.x = 0.819152 * lp.lam * asqrt(1 - xy.x * xy.x);
-	return (xy);
-}
-INVERSE(s_inverse); /* spheroid */
-	(void) P;
-	lp.phi = 2. * atan(xy.y /= 1.819152);
-	lp.lam = fabs(xy.y = 1. - xy.y * xy.y) < TOL ? 0. :
-		xy.x / (0.819152 * sqrt(xy.y));
-	return (lp);
-}
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(fahey) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
+#define PJ_LIB__

+# include   <projects.h>

+

+PROJ_HEAD(fahey, "Fahey") "\n\tPcyl, Sph.";

+

+#define TOL 1e-6

+

+

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

+    XY xy = {0.0,0.0};

+    (void) P;

+

+    xy.x = tan(0.5 * lp.phi);

+    xy.y = 1.819152 * xy.x;

+    xy.x = 0.819152 * lp.lam * asqrt(1 - xy.x * xy.x);

+    return xy;

+}

+

+

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

+    LP lp = {0.0,0.0};

+    (void) P;

+

+    xy.y /= 1.819152;

+    lp.phi = 2. * atan(xy.y);

+    xy.y = 1. - xy.y * xy.y;

+    lp.lam = fabs(xy.y) < TOL ? 0. : xy.x / (0.819152 * sqrt(xy.y));

+    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(fahey) {

+    P->es = 0.;

+    P->inv = s_inverse;

+    P->fwd = s_forward;

+

+    return P;

+}

+

+#ifdef PJ_OMIT_SELFTEST

+int pj_fahey_selftest (void) {return 0;}

+#else

+

+int pj_fahey_selftest (void) {

+    double tolerance_lp = 1e-10;

+    double tolerance_xy = 1e-7;

+

+    char e_args[] = {"+proj=fahey   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};

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

+

+    LP fwd_in[] = {

+        { 2, 1},

+        { 2,-1},

+        {-2, 1},

+        {-2,-1}

+    };

+

+    XY s_fwd_expect[] = {

+        { 182993.34464912376,  101603.19356988439
},

+        { 182993.34464912376, -101603.19356988439
},

+        {-182993.34464912376,  101603.19356988439
},

+        {-182993.34464912376, -101603.19356988439
},

+    };

+

+    XY inv_in[] = {

+        { 200, 100},

+        { 200,-100},

+        {-200, 100},

+        {-200,-100}

+    };

+

+    LP s_inv_expect[] = {

+        {0.0021857886080359551,  0.00098424601668238403
},

+        {0.0021857886080359551,  -0.00098424601668238403
},

+        {-0.0021857886080359551,  0.00098424601668238403
},

+        {-0.0021857886080359551,  -0.00098424601668238403
},

+    };

+

+    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_fouc_s.c b/src/PJ_fouc_s.c
index b84b3f82..4123497c 100644
--- a/src/PJ_fouc_s.c
+++ b/src/PJ_fouc_s.c
@@ -1,45 +1,126 @@
-#define PROJ_PARMS__ \
-	double n, n1;
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(fouc_s, "Foucaut Sinusoidal") "\n\tPCyl., Sph.";
+
 #define MAX_ITER    10
 #define LOOP_TOL    1e-7
-FORWARD(s_forward); /* spheroid */
-	double t;
 
-	t = cos(lp.phi);
-	xy.x = lp.lam * t / (P->n + P->n1 * t);
-	xy.y = P->n * lp.phi + P->n1 * sin(lp.phi);
-	return (xy);
+struct pj_opaque {
+    double n, n1;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double t;
+
+    t = cos(lp.phi);
+    xy.x = lp.lam * t / (Q->n + Q->n1 * t);
+    xy.y = Q->n * lp.phi + Q->n1 * 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 V;
+    int i;
+
+    if (Q->n) {
+        lp.phi = xy.y;
+        for (i = MAX_ITER; i ; --i) {
+            lp.phi -= V = (Q->n * lp.phi + Q->n1 * sin(lp.phi) - xy.y ) /
+                (Q->n + Q->n1 * cos(lp.phi));
+            if (fabs(V) < LOOP_TOL)
+                break;
+        }
+        if (!i)
+            lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
+    } else
+        lp.phi = aasin(P->ctx,xy.y);
+    V = cos(lp.phi);
+    lp.lam = xy.x * (Q->n + Q->n1 * V) / V;
+    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 V;
-	int i;
-
-	if (P->n) {
-		lp.phi = xy.y;
-		for (i = MAX_ITER; i ; --i) {
-			lp.phi -= V = (P->n * lp.phi + P->n1 * sin(lp.phi) - xy.y ) /
-				(P->n + P->n1 * cos(lp.phi));
-			if (fabs(V) < LOOP_TOL)
-				break;
-		}
-		if (!i)
-			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
-	} else
-		lp.phi = aasin(P->ctx,xy.y);
-	V = cos(lp.phi);
-	lp.lam = xy.x * (P->n + P->n1 * V) / V;
-	return (lp);
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(fouc_s)
-	P->n = pj_param(P->ctx, P->params, "dn").f;
-	if (P->n < 0. || P->n > 1.)
-		E_ERROR(-99)
-	P->n1 = 1. - P->n;
-	P->es = 0;
-	P->inv = s_inverse;
-	P->fwd = s_forward;
-ENDENTRY(P)
+
+
+PJ *PROJECTION(fouc_s) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    Q->n = pj_param(P->ctx, P->params, "dn").f;
+    if (Q->n < 0. || Q->n > 1.)
+        E_ERROR(-99)
+    Q->n1 = 1. - Q->n;
+    P->es = 0;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_fouc_s_selftest (void) {return 0;}
+#else
+
+int pj_fouc_s_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=fouc_s   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223402.14425527424,  111695.40119861449},
+        { 223402.14425527424, -111695.40119861449},
+        {-223402.14425527424,  111695.40119861449},
+        {-223402.14425527424, -111695.40119861449},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931097838226,  0.000895246554928339},
+        { 0.0017904931097838226, -0.000895246554928339},
+        {-0.0017904931097838226,  0.000895246554928339},
+        {-0.0017904931097838226, -0.000895246554928339},
+    };
+
+    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_gall.c b/src/PJ_gall.c
index 1acd7d09..b3d7cc6c 100644
--- a/src/PJ_gall.c
+++ b/src/PJ_gall.c
@@ -1,21 +1,100 @@
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(gall, "Gall (Gall Stereographic)") "\n\tCyl, Sph";
-#define YF	1.70710678118654752440
-#define XF	0.70710678118654752440
-#define RYF	0.58578643762690495119
-#define RXF	1.41421356237309504880
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	xy.x = XF * lp.lam;
-	xy.y = YF * tan(.5 * lp.phi);
-	return (xy);
+
+#define YF  1.70710678118654752440
+#define XF  0.70710678118654752440
+#define RYF 0.58578643762690495119
+#define RXF 1.41421356237309504880
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    (void) P;
+
+    xy.x = XF * lp.lam;
+    xy.y = YF * tan(.5 * lp.phi);
+
+    return xy;
 }
-INVERSE(s_inverse); /* spheroid */
-	(void) P;
-	lp.lam = RXF * xy.x;
-	lp.phi = 2. * atan(xy.y * RYF);
-	return (lp);
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    (void) P;
+
+    lp.lam = RXF * xy.x;
+    lp.phi = 2. * atan(xy.y * RYF);
+
+    return lp;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(gall) P->es = 0.; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
+
+
+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(gall) {
+    P->es = 0.0;
+
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_gall_selftest (void) {return 0;}
+#else
+
+int pj_gall_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=gall   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 157969.17113451968,  95345.249178385886},
+        { 157969.17113451968, -95345.249178385886},
+        {-157969.17113451968,  95345.249178385886},
+        {-157969.17113451968, -95345.249178385886},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0025321396391918614,  0.001048846580346495},
+        { 0.0025321396391918614, -0.001048846580346495},
+        {-0.0025321396391918614,  0.001048846580346495},
+        {-0.0025321396391918614, -0.001048846580346495},
+    };
+
+    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_geos.c b/src/PJ_geos.c
index 09393adf..578608c9 100644
--- a/src/PJ_geos.c
+++ b/src/PJ_geos.c
@@ -3,8 +3,7 @@
 **
 ** Copyright (c) 2004   Gerald I. Evenden
 ** Copyright (c) 2012   Martin Raspaud
-*/
-/*
+**
 ** See also (section 4.4.3.2):
 **   http://www.eumetsat.int/en/area4/msg/news/us_doc/cgms_03_26.pdf
 **
@@ -27,164 +26,274 @@
 ** TORT OR OTHERWISE, ARISING FROM, OUT OF OR IN CONNECTION WITH THE
 ** SOFTWARE OR THE USE OR OTHER DEALINGS IN THE SOFTWARE.
 */
-#define PROJ_PARMS__ \
-	double	h; \
-	double  radius_p; \
-	double  radius_p2; \
-	double  radius_p_inv2; \
-	double  radius_g; \
-	double  radius_g_1; \
-	double  C; \
-        char *  sweep_axis; \
-        int     flip_axis;
+
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
+struct pj_opaque {
+    double h;
+    double radius_p;
+    double radius_p2;
+    double radius_p_inv2;
+    double radius_g;
+    double radius_g_1;
+    double C;
+    char *sweep_axis;
+    int flip_axis;
+};
 
 PROJ_HEAD(geos, "Geostationary Satellite View") "\n\tAzi, Sph&Ell\n\th=";
 
-FORWARD(s_forward); /* spheroid */
-	double Vx, Vy, Vz, tmp;
-
-/* Calculation of the three components of the vector from satellite to
-** position on earth surface (lon,lat).*/
-	tmp = cos(lp.phi);
-	Vx = cos (lp.lam) * tmp;
-	Vy = sin (lp.lam) * tmp;
-	Vz = sin (lp.phi);
-/* Check visibility.*/
-	if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz) < 0.) F_ERROR;
-/* Calculation based on view angles from satellite.*/
-	tmp = P->radius_g - Vx;
-        if(P->flip_axis)
-          {
-            xy.x = P->radius_g_1 * atan(Vy / hypot(Vz, tmp));
-            xy.y = P->radius_g_1 * atan(Vz / tmp);
-          }
-        else
-          {
-            xy.x = P->radius_g_1 * atan(Vy / tmp);
-            xy.y = P->radius_g_1 * atan(Vz / hypot(Vy, tmp));
-          }
-	return (xy);
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double Vx, Vy, Vz, tmp;
+
+    /* Calculation of the three components of the vector from satellite to
+    ** position on earth surface (lon,lat).*/
+    tmp = cos(lp.phi);
+    Vx = cos (lp.lam) * tmp;
+    Vy = sin (lp.lam) * tmp;
+    Vz = sin (lp.phi);
+
+    /* Check visibility*/
+
+
+    /* Calculation based on view angles from satellite.*/
+    tmp = Q->radius_g - Vx;
+
+    if(Q->flip_axis) {
+        xy.x = Q->radius_g_1 * atan(Vy / hypot(Vz, tmp));
+        xy.y = Q->radius_g_1 * atan(Vz / tmp);
+    } else {
+        xy.x = Q->radius_g_1 * atan(Vy / tmp);
+        xy.y = Q->radius_g_1 * atan(Vz / hypot(Vy, tmp));
+    }
+
+    return xy;
 }
-FORWARD(e_forward); /* ellipsoid */
-	double r, Vx, Vy, Vz, tmp;
-
-/* Calculation of geocentric latitude. */
-	lp.phi = atan (P->radius_p2 * tan (lp.phi));
-/* Calculation of the three components of the vector from satellite to
-** position on earth surface (lon,lat).*/
-	r = (P->radius_p) / hypot(P->radius_p * cos (lp.phi), sin (lp.phi));
-	Vx = r * cos (lp.lam) * cos (lp.phi);
-	Vy = r * sin (lp.lam) * cos (lp.phi);
-	Vz = r * sin (lp.phi);
-/* Check visibility. */
-	if (((P->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * P->radius_p_inv2) < 0.)
-		F_ERROR;
-/* Calculation based on view angles from satellite. */
-	tmp = P->radius_g - Vx;
-        if(P->flip_axis)
-          {
-            xy.x = P->radius_g_1 * atan (Vy / hypot (Vz, tmp));
-            xy.y = P->radius_g_1 * atan (Vz / tmp);
-          }
-        else
-          {
-            xy.x = P->radius_g_1 * atan (Vy / tmp);
-            xy.y = P->radius_g_1 * atan (Vz / hypot (Vy, tmp));
-          }
-	return (xy);
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double r, Vx, Vy, Vz, tmp;
+
+    /* Calculation of geocentric latitude. */
+    lp.phi = atan (Q->radius_p2 * tan (lp.phi));
+
+    /* Calculation of the three components of the vector from satellite to
+    ** position on earth surface (lon,lat).*/
+    r = (Q->radius_p) / hypot(Q->radius_p * cos (lp.phi), sin (lp.phi));
+    Vx = r * cos (lp.lam) * cos (lp.phi);
+    Vy = r * sin (lp.lam) * cos (lp.phi);
+    Vz = r * sin (lp.phi);
+
+    /* Check visibility. */
+    if (((Q->radius_g - Vx) * Vx - Vy * Vy - Vz * Vz * Q->radius_p_inv2) < 0.)
+        F_ERROR;
+
+    /* Calculation based on view angles from satellite. */
+    tmp = Q->radius_g - Vx;
+
+    if(Q->flip_axis) {
+        xy.x = Q->radius_g_1 * atan (Vy / hypot (Vz, tmp));
+        xy.y = Q->radius_g_1 * atan (Vz / tmp);
+    } else {
+        xy.x = Q->radius_g_1 * atan (Vy / tmp);
+        xy.y = Q->radius_g_1 * atan (Vz / hypot (Vy, tmp));
+    }
+
+    return xy;
 }
-INVERSE(s_inverse); /* spheroid */
-	double Vx, Vy, Vz, a, b, det, k;
-
-/* Setting three components of vector from satellite to position.*/
-	Vx = -1.0;
-	if(P->flip_axis)
-          {
-            Vz = tan (xy.y / (P->radius_g - 1.0));
-            Vy = tan (xy.x / (P->radius_g - 1.0)) * sqrt (1.0 + Vz * Vz);
-          }
-        else
-          {
-            Vy = tan (xy.x / (P->radius_g - 1.0));
-            Vz = tan (xy.y / (P->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
-          }
-/* Calculation of terms in cubic equation and determinant.*/
-	a   = Vy * Vy + Vz * Vz + Vx * Vx;
-	b   = 2 * P->radius_g * Vx;
-	if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
-/* Calculation of three components of vector from satellite to position.*/
-	k  = (-b - sqrt(det)) / (2 * a);
-	Vx = P->radius_g + k * Vx;
-	Vy *= k;
-	Vz *= k;
-/* Calculation of longitude and latitude.*/
-	lp.lam = atan2 (Vy, Vx);
-	lp.phi = atan (Vz * cos (lp.lam) / Vx);
-	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 Vx, Vy, Vz, a, b, det, k;
+
+    /* Setting three components of vector from satellite to position.*/
+    Vx = -1.0;
+    if(Q->flip_axis) {
+        Vz = tan (xy.y / (Q->radius_g - 1.0));
+        Vy = tan (xy.x / (Q->radius_g - 1.0)) * sqrt (1.0 + Vz * Vz);
+    } else {
+        Vy = tan (xy.x / (Q->radius_g - 1.0));
+        Vz = tan (xy.y / (Q->radius_g - 1.0)) * sqrt (1.0 + Vy * Vy);
+    }
+
+    /* Calculation of terms in cubic equation and determinant.*/
+    a = Vy * Vy + Vz * Vz + Vx * Vx;
+    b = 2 * Q->radius_g * Vx;
+    if ((det = (b * b) - 4 * a * Q->C) < 0.) I_ERROR;
+
+    /* Calculation of three components of vector from satellite to position.*/
+    k  = (-b - sqrt(det)) / (2 * a);
+    Vx = Q->radius_g + k * Vx;
+    Vy *= k;
+    Vz *= k;
+
+    /* Calculation of longitude and latitude.*/
+    lp.lam = atan2 (Vy, Vx);
+    lp.phi = atan (Vz * cos (lp.lam) / Vx);
+
+    return lp;
 }
-INVERSE(e_inverse); /* ellipsoid */
-	double Vx, Vy, Vz, a, b, det, k;
-
-/* Setting three components of vector from satellite to position.*/
-	Vx = -1.0;
-        if(P->flip_axis)
-          {
-            Vz = tan (xy.y / P->radius_g_1);
-            Vy = tan (xy.x / P->radius_g_1) * hypot(1.0, Vz);
-          }
-        else
-          {
-            Vy = tan (xy.x / P->radius_g_1);
-            Vz = tan (xy.y / P->radius_g_1) * hypot(1.0, Vy);
-          }
-/* Calculation of terms in cubic equation and determinant.*/
-	a = Vz / P->radius_p;
-	a   = Vy * Vy + a * a + Vx * Vx;
-	b   = 2 * P->radius_g * Vx;
-	if ((det = (b * b) - 4 * a * P->C) < 0.) I_ERROR;
-/* Calculation of three components of vector from satellite to position.*/
-	k  = (-b - sqrt(det)) / (2. * a);
-	Vx = P->radius_g + k * Vx;
-	Vy *= k;
-	Vz *= k;
-/* Calculation of longitude and latitude.*/
-	lp.lam  = atan2 (Vy, Vx);
-	lp.phi = atan (Vz * cos (lp.lam) / Vx);
-	lp.phi = atan (P->radius_p_inv2 * tan (lp.phi));
-	return (lp);
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double Vx, Vy, Vz, a, b, det, k;
+
+    /* Setting three components of vector from satellite to position.*/
+    Vx = -1.0;
+
+    if(Q->flip_axis) {
+        Vz = tan (xy.y / Q->radius_g_1);
+        Vy = tan (xy.x / Q->radius_g_1) * hypot(1.0, Vz);
+    } else {
+        Vy = tan (xy.x / Q->radius_g_1);
+        Vz = tan (xy.y / Q->radius_g_1) * hypot(1.0, Vy);
+    }
+
+    /* Calculation of terms in cubic equation and determinant.*/
+    a = Vz / Q->radius_p;
+    a   = Vy * Vy + a * a + Vx * Vx;
+    b   = 2 * Q->radius_g * Vx;
+    if ((det = (b * b) - 4 * a * Q->C) < 0.) I_ERROR;
+
+    /* Calculation of three components of vector from satellite to position.*/
+    k  = (-b - sqrt(det)) / (2. * a);
+    Vx = Q->radius_g + k * Vx;
+    Vy *= k;
+    Vz *= k;
+
+    /* Calculation of longitude and latitude.*/
+    lp.lam = atan2 (Vy, Vx);
+    lp.phi = atan (Vz * cos (lp.lam) / Vx);
+    lp.phi = atan (Q->radius_p_inv2 * tan (lp.phi));
+
+    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);
 }
-FREEUP; if (P) free(P); }
-ENTRY0(geos)
-	if ((P->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30);
-	if (P->phi0) E_ERROR(-46);
-        P->sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
-        if (P->sweep_axis == NULL)
-          P->flip_axis = 0;
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(geos) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    if ((Q->h = pj_param(P->ctx, P->params, "dh").f) <= 0.) E_ERROR(-30);
+
+    if (P->phi0) E_ERROR(-46);
+
+    Q->sweep_axis = pj_param(P->ctx, P->params, "ssweep").s;
+    if (Q->sweep_axis == NULL)
+      Q->flip_axis = 0;
+    else {
+        if (Q->sweep_axis[1] != '\0' ||
+            (Q->sweep_axis[0] != 'x' &&
+             Q->sweep_axis[0] != 'y'))
+          E_ERROR(-49);
+        if (Q->sweep_axis[0] == 'x')
+          Q->flip_axis = 1;
         else
-          {
-            if (P->sweep_axis[1] != '\0' ||
-                (P->sweep_axis[0] != 'x' &&
-                 P->sweep_axis[0] != 'y'))
-              E_ERROR(-49);
-            if (P->sweep_axis[0] == 'x')
-              P->flip_axis = 1;
-            else
-              P->flip_axis = 0;
-          }
-        P->radius_g_1 = P->h / P->a;
-	P->radius_g = 1. + P->radius_g_1;
-	P->C  = P->radius_g * P->radius_g - 1.0;
-	if (P->es) {
-		P->radius_p      = sqrt (P->one_es);
-		P->radius_p2     = P->one_es;
-		P->radius_p_inv2 = P->rone_es;
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
-		P->radius_p = P->radius_p2 = P->radius_p_inv2 = 1.0;
-		P->inv = s_inverse;
-		P->fwd = s_forward;
-	}
-ENDENTRY(P)
+          Q->flip_axis = 0;
+    }
+
+    Q->radius_g_1 = Q->h / P->a;
+    Q->radius_g = 1. + Q->radius_g_1;
+    Q->C  = Q->radius_g * Q->radius_g - 1.0;
+    if (P->es) {
+        Q->radius_p      = sqrt (P->one_es);
+        Q->radius_p2     = P->one_es;
+        Q->radius_p_inv2 = P->rone_es;
+        P->inv = e_inverse;
+        P->fwd = e_forward;
+    } else {
+        Q->radius_p = Q->radius_p2 = Q->radius_p_inv2 = 1.0;
+        P->inv = s_inverse;
+        P->fwd = s_forward;
+    }
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_geos_selftest (void) {return 0;}
+#else
+
+int pj_geos_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=geos   +ellps=GRS80  +lat_1=0.5 +lat_2=2 +h=35785831"};
+    char s_args[] = {"+proj=geos   +a=6400000    +lat_1=0.5 +lat_2=2 +h=35785831"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        { 222527.07036580026,  110551.30341332949},
+        { 222527.07036580026, -110551.30341332949},
+        {-222527.07036580026,  110551.30341332949},
+        {-222527.07036580026, -110551.30341332949},
+    };
+
+    XY s_fwd_expect[] = {
+        { 223289.45763579503,  111677.65745653701},
+        { 223289.45763579503,  -111677.65745653701},
+        {-223289.45763579503,  111677.65745653701},
+        {-223289.45763579503,  -111677.65745653701},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        { 0.0017966305689715385,  0.00090436947723267452},
+        { 0.0017966305689715385, -0.00090436947723267452},
+        {-0.0017966305689715385,  0.00090436947723267452},
+        {-0.0017966305689715385, -0.00090436947723267452},
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931105078943,  0.00089524655504237148},
+        { 0.0017904931105078943, -0.00089524655504237148},
+        {-0.0017904931105078943,  0.00089524655504237148},
+        {-0.0017904931105078943, -0.00089524655504237148},
+    };
+
+    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_gins8.c b/src/PJ_gins8.c
index b7f38ad9..3eae7efa 100644
--- a/src/PJ_gins8.c
+++ b/src/PJ_gins8.c
@@ -1,18 +1,75 @@
 #define PJ_LIB__
-# include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(gins8, "Ginsburg VIII (TsNIIGAiK)") "\n\tPCyl, Sph., no inv.";
-#define	Cl 0.000952426
-#define	Cp 0.162388
-#define	C12 0.08333333333333333
-FORWARD(s_forward); /* spheroid */
-	double t = lp.phi * lp.phi;
-	(void) P;
-
-	xy.y = lp.phi * (1. + t * C12);
-	xy.x = lp.lam * (1. - Cp * t);
-	t = lp.lam * lp.lam;
-	xy.x *= (0.87 - Cl * t * t);
-	return (xy);
+
+#define Cl 0.000952426
+#define Cp 0.162388
+#define C12 0.08333333333333333
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    double t = lp.phi * lp.phi;
+    (void) P;
+
+    xy.y = lp.phi * (1. + t * C12);
+    xy.x = lp.lam * (1. - Cp * t);
+    t = lp.lam * lp.lam;
+    xy.x *= (0.87 - Cl * t * t);
+
+    return xy;
+}
+
+
+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(gins8) {
+    P->es = 0.0;
+    P->inv = 0;
+    P->fwd = s_forward;
+
+    return P;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(gins8) P->es = 0.; P->inv = 0; P->fwd = s_forward; ENDENTRY(P)
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_gins8_selftest (void) {return 0;}
+#else
+
+int pj_gins8_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=gins8   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 194350.25093959007,  111703.90763533533},
+        { 194350.25093959007, -111703.90763533533},
+        {-194350.25093959007,  111703.90763533533},
+        {-194350.25093959007, -111703.90763533533},
+    };
+
+    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_gnom.c b/src/PJ_gnom.c
index 11deb86c..151f718e 100644
--- a/src/PJ_gnom.c
+++ b/src/PJ_gnom.c
@@ -1,105 +1,192 @@
-#define PROJ_PARMS__ \
-	double	sinph0; \
-	double	cosph0; \
-	int		mode;
 #define PJ_LIB__
-#include	<projects.h>
+#include <projects.h>
+
 PROJ_HEAD(gnom, "Gnomonic") "\n\tAzi, Sph.";
-#define EPS10 1.e-10
-#define N_POLE	0
+
+#define EPS10  1.e-10
+#define N_POLE 0
 #define S_POLE 1
-#define EQUIT	2
-#define OBLIQ	3
-FORWARD(s_forward); /* spheroid */
-	double  coslam, cosphi, sinphi;
-
-	sinphi = sin(lp.phi);
-	cosphi = cos(lp.phi);
-	coslam = cos(lp.lam);
-	switch (P->mode) {
-	case EQUIT:
-		xy.y = cosphi * coslam;
-		break;
-	case OBLIQ:
-		xy.y = P->sinph0 * sinphi + P->cosph0 * cosphi * coslam;
-		break;
-	case S_POLE:
-		xy.y = - sinphi;
-		break;
-	case N_POLE:
-		xy.y = sinphi;
-		break;
-	}
-	if (xy.y <= EPS10) F_ERROR;
-	xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam);
-	switch (P->mode) {
-	case EQUIT:
-		xy.y *= sinphi;
-		break;
-	case OBLIQ:
-		xy.y *= P->cosph0 * sinphi - P->sinph0 * cosphi * coslam;
-		break;
-	case N_POLE:
-		coslam = - coslam;
-	case S_POLE:
-		xy.y *= cosphi * coslam;
-		break;
-	}
-	return (xy);
+#define EQUIT  2
+#define OBLIQ  3
+
+struct pj_opaque {
+    double  sinph0;
+    double  cosph0;
+    int     mode;
+};
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double  coslam, cosphi, sinphi;
+
+    sinphi = sin(lp.phi);
+    cosphi = cos(lp.phi);
+    coslam = cos(lp.lam);
+
+    switch (Q->mode) {
+        case EQUIT:
+            xy.y = cosphi * coslam;
+            break;
+        case OBLIQ:
+            xy.y = Q->sinph0 * sinphi + Q->cosph0 * cosphi * coslam;
+            break;
+        case S_POLE:
+            xy.y = - sinphi;
+            break;
+        case N_POLE:
+            xy.y = sinphi;
+            break;
+    }
+
+    if (xy.y <= EPS10) F_ERROR;
+
+    xy.x = (xy.y = 1. / xy.y) * cosphi * sin(lp.lam);
+    switch (Q->mode) {
+        case EQUIT:
+            xy.y *= sinphi;
+            break;
+        case OBLIQ:
+            xy.y *= Q->cosph0 * sinphi - Q->sinph0 * cosphi * coslam;
+            break;
+        case N_POLE:
+            coslam = - coslam;
+        case S_POLE:
+            xy.y *= cosphi * coslam;
+            break;
+    }
+    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  rh, cosz, sinz;
+
+    rh = hypot(xy.x, xy.y);
+    sinz = sin(lp.phi = atan(rh));
+    cosz = sqrt(1. - sinz * sinz);
+
+    if (fabs(rh) <= EPS10) {
+        lp.phi = P->phi0;
+        lp.lam = 0.;
+    } else {
+        switch (Q->mode) {
+            case OBLIQ:
+                lp.phi = cosz * Q->sinph0 + xy.y * sinz * Q->cosph0 / rh;
+                if (fabs(lp.phi) >= 1.)
+                    lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
+                else
+                    lp.phi = asin(lp.phi);
+                xy.y = (cosz - Q->sinph0 * sin(lp.phi)) * rh;
+                xy.x *= sinz * Q->cosph0;
+                break;
+            case EQUIT:
+                lp.phi = xy.y * sinz / rh;
+                if (fabs(lp.phi) >= 1.)
+                    lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
+                else
+                    lp.phi = asin(lp.phi);
+                xy.y = cosz * rh;
+                xy.x *= sinz;
+                break;
+            case S_POLE:
+                lp.phi -= HALFPI;
+                break;
+            case N_POLE:
+                lp.phi = HALFPI - lp.phi;
+                xy.y = -xy.y;
+                break;
+        }
+        lp.lam = atan2(xy.x, xy.y);
+    }
+    return lp;
+}
+
+
+static void *freeup_new (PJ *P) {              /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    pj_dealloc (P->opaque);
+    return pj_dealloc(P);
+}
+
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-INVERSE(s_inverse); /* spheroid */
-	double  rh, cosz, sinz;
-
-	rh = hypot(xy.x, xy.y);
-	sinz = sin(lp.phi = atan(rh));
-	cosz = sqrt(1. - sinz * sinz);
-	if (fabs(rh) <= EPS10) {
-		lp.phi = P->phi0;
-		lp.lam = 0.;
-	} else {
-		switch (P->mode) {
-		case OBLIQ:
-			lp.phi = cosz * P->sinph0 + xy.y * sinz * P->cosph0 / rh;
-			if (fabs(lp.phi) >= 1.)
-				lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
-			else
-				lp.phi = asin(lp.phi);
-			xy.y = (cosz - P->sinph0 * sin(lp.phi)) * rh;
-			xy.x *= sinz * P->cosph0;
-			break;
-		case EQUIT:
-			lp.phi = xy.y * sinz / rh;
-			if (fabs(lp.phi) >= 1.)
-				lp.phi = lp.phi > 0. ? HALFPI : - HALFPI;
-			else
-				lp.phi = asin(lp.phi);
-			xy.y = cosz * rh;
-			xy.x *= sinz;
-			break;
-		case S_POLE:
-			lp.phi -= HALFPI;
-			break;
-		case N_POLE:
-			lp.phi = HALFPI - lp.phi;
-			xy.y = -xy.y;
-			break;
-		}
-		lp.lam = atan2(xy.x, xy.y);
-	}
-	return (lp);
+
+
+PJ *PROJECTION(gnom) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    if (fabs(fabs(P->phi0) - HALFPI) < EPS10) {
+        Q->mode = P->phi0 < 0. ? S_POLE : N_POLE;
+    } else if (fabs(P->phi0) < EPS10) {
+        Q->mode = EQUIT;
+    } else {
+        Q->mode = OBLIQ;
+        Q->sinph0 = sin(P->phi0);
+        Q->cosph0 = cos(P->phi0);
+    }
+
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+    P->es = 0.;
+
+    return P;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(gnom)
-	if (fabs(fabs(P->phi0) - HALFPI) < EPS10)
-		P->mode = P->phi0 < 0. ? S_POLE : N_POLE;
-	else if (fabs(P->phi0) < EPS10)
-		P->mode = EQUIT;
-	else {
-		P->mode = OBLIQ;
-		P->sinph0 = sin(P->phi0);
-		P->cosph0 = cos(P->phi0);
-	}
-	P->inv = s_inverse;
-	P->fwd = s_forward;
-	P->es = 0.;
-ENDENTRY(P)
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_gnom_selftest (void) {return 0;}
+#else
+
+int pj_gnom_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=gnom   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        { 223492.92474718543,  111780.50920659291},
+        { 223492.92474718543, -111780.50920659291},
+        {-223492.92474718543,  111780.50920659291},
+        {-223492.92474718543, -111780.50920659291},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        { 0.0017904931092009798,  0.00089524655438192376},
+        { 0.0017904931092009798, -0.00089524655438192376},
+        {-0.0017904931092009798,  0.00089524655438192376},
+        {-0.0017904931092009798, -0.00089524655438192376},
+    };
+
+    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