Merge pull request #3 from kbevers/fix-projections-with-c

Projections starting with the letter c refactored

6 files changed, 589 insertions, 227 deletions

diff --git a/src/PJ_aea.c b/src/PJ_aea.c
index bd6c90ac..0822efd6 100644
--- a/src/PJ_aea.c
+++ b/src/PJ_aea.c
@@ -355,12 +355,6 @@ source files
 
 int pj_aeqd_selftest     (void)     {return 10000;}
 int pj_alsk_selftest     (void)     {return 10000;}
-
-int pj_cea_selftest      (void)     {return 10000;}
-int pj_chamb_selftest    (void)     {return 10000;}
-int pj_collg_selftest    (void)     {return 10000;}
-int pj_comill_selftest   (void)     {return 10000;}
-int pj_crast_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;}
diff --git a/src/PJ_cea.c b/src/PJ_cea.c
index 44c0a887..87c0da07 100644
--- a/src/PJ_cea.c
+++ b/src/PJ_cea.c
@@ -1,63 +1,157 @@
-#define PROJ_PARMS__ \
-	double qp; \
-	double *apa;
 #define PJ_LIB__
-# include	<projects.h>
+#include   <projects.h>
+
+struct pj_opaque {
+    double qp;
+    double *apa;
+};
+
 PROJ_HEAD(cea, "Equal Area Cylindrical") "\n\tCyl, Sph&Ell\n\tlat_ts=";
-# define EPS	1e-10
-FORWARD(e_forward); /* spheroid */
-	xy.x = P->k0 * lp.lam;
-	xy.y = .5 * pj_qsfn(sin(lp.phi), P->e, P->one_es) / P->k0;
-	return (xy);
-}
-FORWARD(s_forward); /* spheroid */
-	xy.x = P->k0 * lp.lam;
-	xy.y = sin(lp.phi) / P->k0;
-	return (xy);
-}
-INVERSE(e_inverse); /* spheroid */
-	lp.phi = pj_authlat(asin( 2. * xy.y * P->k0 / P->qp), P->apa);
-	lp.lam = xy.x / P->k0;
-	return (lp);
-}
-INVERSE(s_inverse); /* spheroid */
-	double t;
-
-	if ((t = fabs(xy.y *= P->k0)) - EPS <= 1.) {
-		if (t >= 1.)
-			lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
-		else
-			lp.phi = asin(xy.y);
-		lp.lam = xy.x / P->k0;
-	} else I_ERROR;
-	return (lp);
-}
-FREEUP;
-	if (P) {
-		if (P->apa)
-			pj_dalloc(P->apa);
-		pj_dalloc(P);
-	}
-}
-ENTRY1(cea, apa)
-	double t = 0.0;
-
-	if (pj_param(P->ctx, P->params, "tlat_ts").i) {
-		P->k0 = cos(t = pj_param(P->ctx, P->params, "rlat_ts").f);
-		if (P->k0 < 0.) {
-			E_ERROR(-24);
-		}
-	}
-	if (P->es) {
-		t = sin(t);
-		P->k0 /= sqrt(1. - P->es * t * t);
-		P->e = sqrt(P->es);
-		if (!(P->apa = pj_authset(P->es))) E_ERROR_0;
-		P->qp = pj_qsfn(1., P->e, P->one_es);
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
-		P->inv = s_inverse;
-		P->fwd = s_forward;
-	}
-ENDENTRY(P)
+# define EPS    1e-10
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    xy.x = P->k0 * lp.lam;
+    xy.y = 0.5 * pj_qsfn(sin(lp.phi), P->e, P->one_es) / P->k0;
+    return xy;
+}
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    xy.x = P->k0 * lp.lam;
+    xy.y = sin(lp.phi) / P->k0;
+    return xy;
+}
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    lp.phi = pj_authlat(asin( 2. * xy.y * P->k0 / Q->qp), Q->apa);
+    lp.lam = xy.x / P->k0;
+    return lp;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    double t;
+
+    if ((t = fabs(xy.y *= P->k0)) - EPS <= 1.) {
+        if (t >= 1.)
+            lp.phi = xy.y < 0. ? -HALFPI : HALFPI;
+        else
+            lp.phi = asin(xy.y);
+        lp.lam = xy.x / P->k0;
+    } else I_ERROR;
+    return (lp);
+}
+
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    if (0==P)
+        return 0;
+    if (0==P->opaque)
+        return pj_dealloc (P);
+
+    pj_dalloc (P->opaque->apa);
+    pj_dealloc (P->opaque);
+    return pj_dealloc (P);
+}
+
+static void freeup (PJ *P) {
+   freeup_new (P);
+    return;
+}
+
+PJ *PROJECTION(cea) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    double t = 0.0;
+
+    if (pj_param(P->ctx, P->params, "tlat_ts").i) {
+        P->k0 = cos(t = pj_param(P->ctx, P->params, "rlat_ts").f);
+        if (P->k0 < 0.) {
+            E_ERROR(-24);
+        }
+    }
+    if (P->es) {
+        t = sin(t);
+        P->k0 /= sqrt(1. - P->es * t * t);
+        P->e = sqrt(P->es);
+        if (!(Q->apa = pj_authset(P->es))) E_ERROR_0;
+        Q->qp = pj_qsfn(1., P->e, P->one_es);
+        P->inv = e_inverse;
+        P->fwd = e_forward;
+    } else {
+        P->inv = s_inverse;
+        P->fwd = s_forward;
+    }
+
+    return P;
+}
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_cea_selftest (void) {return 0;}
+#else
+
+int pj_cea_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=cea   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=cea   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY e_fwd_expect[] = {
+        {222638.981586547132,  110568.812396267356},
+        {222638.981586547132,  -110568.812396265886},
+        {-222638.981586547132,  110568.812396267356},
+        {-222638.981586547132,  -110568.812396265886},
+    };
+
+    XY s_fwd_expect[] = {
+        {223402.144255274179,  111695.401198614476},
+        {223402.144255274179,  -111695.401198614476},
+        {-223402.144255274179,  111695.401198614476},
+        {-223402.144255274179,  -111695.401198614476},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP e_inv_expect[] = {
+        {0.00179663056823904264,  0.000904369476105564289},
+        {0.00179663056823904264,  -0.000904369476105564289},
+        {-0.00179663056823904264,  0.000904369476105564289},
+        {-0.00179663056823904264,  -0.000904369476105564289},
+    };
+
+    LP s_inv_expect[] = {
+        {0.00179049310978382265,  0.000895246554928338998},
+        {0.00179049310978382265,  -0.000895246554928338998},
+        {-0.00179049310978382265,  0.000895246554928338998},
+        {-0.00179049310978382265,  -0.000895246554928338998},
+    };
+
+    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_chamb.c b/src/PJ_chamb.c
index 65f21129..c35061f3 100644
--- a/src/PJ_chamb.c
+++ b/src/PJ_chamb.c
@@ -1,112 +1,179 @@
-typedef struct { double r, Az; } VECT;
-#define PROJ_PARMS__ \
-	struct { /* control point data */ \
-		double phi, lam; \
-		double cosphi, sinphi; \
-		VECT v; \
-		XY	p; \
-		double Az; \
-	} c[3]; \
-	XY p; \
-	double beta_0, beta_1, beta_2;
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
+
+typedef struct { double r, Az; } VECT;
+struct pj_opaque {
+    struct { /* control point data */
+        double phi, lam;
+        double cosphi, sinphi;
+        VECT v;
+        XY  p;
+        double Az;
+    } c[3];
+    XY p;
+    double beta_0, beta_1, beta_2;
+};
+
 PROJ_HEAD(chamb, "Chamberlin Trimetric") "\n\tMisc Sph, no inv."
 "\n\tlat_1= lon_1= lat_2= lon_2= lat_3= lon_3=";
-#include	<stdio.h>
+
+#include <stdio.h>
 #define THIRD 0.333333333333333333
 #define TOL 1e-9
-	static VECT /* distance and azimuth from point 1 to point 2 */
-vect(projCtx ctx, double dphi, double c1, double s1, double c2, double s2, double dlam) {
-	VECT v;
-	double cdl, dp, dl;
-
-	cdl = cos(dlam);
-	if (fabs(dphi) > 1. || fabs(dlam) > 1.)
-		v.r = aacos(ctx, s1 * s2 + c1 * c2 * cdl);
-	else { /* more accurate for smaller distances */
-		dp = sin(.5 * dphi);
-		dl = sin(.5 * dlam);
-		v.r = 2. * aasin(ctx,sqrt(dp * dp + c1 * c2 * dl * dl));
-	}
-	if (fabs(v.r) > TOL)
-		v.Az = atan2(c2 * sin(dlam), c1 * s2 - s1 * c2 * cdl);
-	else
-		v.r = v.Az = 0.;
-	return v;
+
+/* distance and azimuth from point 1 to point 2 */
+static VECT vect(projCtx ctx, double dphi, double c1, double s1, double c2, double s2, double dlam) {
+    VECT v;
+    double cdl, dp, dl;
+
+    cdl = cos(dlam);
+    if (fabs(dphi) > 1. || fabs(dlam) > 1.)
+        v.r = aacos(ctx, s1 * s2 + c1 * c2 * cdl);
+    else { /* more accurate for smaller distances */
+        dp = sin(.5 * dphi);
+        dl = sin(.5 * dlam);
+        v.r = 2. * aasin(ctx,sqrt(dp * dp + c1 * c2 * dl * dl));
+    }
+    if (fabs(v.r) > TOL)
+        v.Az = atan2(c2 * sin(dlam), c1 * s2 - s1 * c2 * cdl);
+    else
+        v.r = v.Az = 0.;
+    return v;
+}
+
+/* law of cosines */
+static double lc(projCtx ctx, double b,double c,double a) {
+    return aacos(ctx, .5 * (b * b + c * c - a * a) / (b * c));
+}
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double sinphi, cosphi, a;
+    VECT v[3];
+    int i, j;
+
+    sinphi = sin(lp.phi);
+    cosphi = cos(lp.phi);
+    for (i = 0; i < 3; ++i) { /* dist/azimiths from control */
+        v[i] = vect(P->ctx, lp.phi - Q->c[i].phi, Q->c[i].cosphi, Q->c[i].sinphi,
+            cosphi, sinphi, lp.lam - Q->c[i].lam);
+        if ( ! v[i].r)
+            break;
+        v[i].Az = adjlon(v[i].Az - Q->c[i].v.Az);
+    }
+    if (i < 3) /* current point at control point */
+        xy = Q->c[i].p;
+    else { /* point mean of intersepts */
+        xy = Q->p;
+        for (i = 0; i < 3; ++i) {
+            j = i == 2 ? 0 : i + 1;
+            a = lc(P->ctx,Q->c[i].v.r, v[i].r, v[j].r);
+            if (v[i].Az < 0.)
+                a = -a;
+            if (! i) { /* coord comp unique to each arc */
+                xy.x += v[i].r * cos(a);
+                xy.y -= v[i].r * sin(a);
+            } else if (i == 1) {
+                a = Q->beta_1 - a;
+                xy.x -= v[i].r * cos(a);
+                xy.y -= v[i].r * sin(a);
+            } else {
+                a = Q->beta_2 - a;
+                xy.x += v[i].r * cos(a);
+                xy.y += v[i].r * sin(a);
+            }
+        }
+        xy.x *= THIRD; /* mean of arc intercepts */
+        xy.y *= THIRD;
+    }
+    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 double /* law of cosines */
-lc(projCtx ctx, double b,double c,double a) {
-	return aacos(ctx, .5 * (b * b + c * c - a * a) / (b * c));
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
 }
-FORWARD(s_forward); /* spheroid */
-	double sinphi, cosphi, a;
-	VECT v[3];
-	int i, j;
-
-	sinphi = sin(lp.phi);
-	cosphi = cos(lp.phi);
-	for (i = 0; i < 3; ++i) { /* dist/azimiths from control */
-		v[i] = vect(P->ctx, lp.phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
-			cosphi, sinphi, lp.lam - P->c[i].lam);
-		if ( ! v[i].r)
-			break;
-		v[i].Az = adjlon(v[i].Az - P->c[i].v.Az);
-	}
-	if (i < 3) /* current point at control point */
-		xy = P->c[i].p;
-	else { /* point mean of intersepts */
-		xy = P->p;
-		for (i = 0; i < 3; ++i) {
-			j = i == 2 ? 0 : i + 1;
-			a = lc(P->ctx,P->c[i].v.r, v[i].r, v[j].r);
-			if (v[i].Az < 0.)
-				a = -a;
-			if (! i) { /* coord comp unique to each arc */
-				xy.x += v[i].r * cos(a);
-				xy.y -= v[i].r * sin(a);
-			} else if (i == 1) {
-				a = P->beta_1 - a;
-				xy.x -= v[i].r * cos(a);
-				xy.y -= v[i].r * sin(a);
-			} else {
-				a = P->beta_2 - a;
-				xy.x += v[i].r * cos(a);
-				xy.y += v[i].r * sin(a);
-			}
-		}
-		xy.x *= THIRD; /* mean of arc intercepts */
-		xy.y *= THIRD;
-	}
-	return xy;
+
+
+PJ *PROJECTION(chamb) {
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return freeup_new (P);
+    P->opaque = Q;
+
+    int i, j;
+    char line[10];
+
+    for (i = 0; i < 3; ++i) { /* get control point locations */
+        (void)sprintf(line, "rlat_%d", i+1);
+        Q->c[i].phi = pj_param(P->ctx, P->params, line).f;
+        (void)sprintf(line, "rlon_%d", i+1);
+        Q->c[i].lam = pj_param(P->ctx, P->params, line).f;
+        Q->c[i].lam = adjlon(Q->c[i].lam - P->lam0);
+        Q->c[i].cosphi = cos(Q->c[i].phi);
+        Q->c[i].sinphi = sin(Q->c[i].phi);
+    }
+    for (i = 0; i < 3; ++i) { /* inter ctl pt. distances and azimuths */
+        j = i == 2 ? 0 : i + 1;
+        Q->c[i].v = vect(P->ctx,Q->c[j].phi - Q->c[i].phi, Q->c[i].cosphi, Q->c[i].sinphi,
+            Q->c[j].cosphi, Q->c[j].sinphi, Q->c[j].lam - Q->c[i].lam);
+        if (! Q->c[i].v.r) E_ERROR(-25);
+        /* co-linearity problem ignored for now */
+    }
+    Q->beta_0 = lc(P->ctx,Q->c[0].v.r, Q->c[2].v.r, Q->c[1].v.r);
+    Q->beta_1 = lc(P->ctx,Q->c[0].v.r, Q->c[1].v.r, Q->c[2].v.r);
+    Q->beta_2 = PI - Q->beta_0;
+    Q->p.y = 2. * (Q->c[0].p.y = Q->c[1].p.y = Q->c[2].v.r * sin(Q->beta_0));
+    Q->c[2].p.y = 0.;
+    Q->c[0].p.x = - (Q->c[1].p.x = 0.5 * Q->c[0].v.r);
+    Q->p.x = Q->c[2].p.x = Q->c[0].p.x + Q->c[2].v.r * cos(Q->beta_0);
+
+    P->es = 0.;
+    P->fwd = s_forward;
+
+    return P;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(chamb)
-	int i, j;
-	char line[10];
-
-	for (i = 0; i < 3; ++i) { /* get control point locations */
-		(void)sprintf(line, "rlat_%d", i+1);
-		P->c[i].phi = pj_param(P->ctx, P->params, line).f;
-		(void)sprintf(line, "rlon_%d", i+1);
-		P->c[i].lam = pj_param(P->ctx, P->params, line).f;
-		P->c[i].lam = adjlon(P->c[i].lam - P->lam0);
-		P->c[i].cosphi = cos(P->c[i].phi);
-		P->c[i].sinphi = sin(P->c[i].phi);
-	}
-	for (i = 0; i < 3; ++i) { /* inter ctl pt. distances and azimuths */
-		j = i == 2 ? 0 : i + 1;
-		P->c[i].v = vect(P->ctx,P->c[j].phi - P->c[i].phi, P->c[i].cosphi, P->c[i].sinphi,
-			P->c[j].cosphi, P->c[j].sinphi, P->c[j].lam - P->c[i].lam);
-		if (! P->c[i].v.r) E_ERROR(-25);
-		/* co-linearity problem ignored for now */
-	}
-	P->beta_0 = lc(P->ctx,P->c[0].v.r, P->c[2].v.r, P->c[1].v.r);
-	P->beta_1 = lc(P->ctx,P->c[0].v.r, P->c[1].v.r, P->c[2].v.r);
-	P->beta_2 = PI - P->beta_0;
-	P->p.y = 2. * (P->c[0].p.y = P->c[1].p.y = P->c[2].v.r * sin(P->beta_0));
-	P->c[2].p.y = 0.;
-	P->c[0].p.x = - (P->c[1].p.x = 0.5 * P->c[0].v.r);
-	P->p.x = P->c[2].p.x = P->c[0].p.x + P->c[2].v.r * cos(P->beta_0);
-	P->es = 0.; P->fwd = s_forward;
-ENDENTRY(P)
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_chamb_selftest (void) {return 0;}
+#else
+
+int pj_chamb_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=chamb   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        {-27864.7795868005815,  -223364.324593274243},
+        {-251312.283053493476,  -223402.145526208304},
+        {-27864.7856491046077,  223364.327328827145},
+        {-251312.289116443484,  223402.142197287147},
+    };
+
+    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_collg.c b/src/PJ_collg.c
index 871dfc97..80029a3a 100644
--- a/src/PJ_collg.c
+++ b/src/PJ_collg.c
@@ -1,10 +1,14 @@
 #define PJ_LIB__
 # include	<projects.h>
+
 PROJ_HEAD(collg, "Collignon") "\n\tPCyl, Sph.";
 #define FXC	1.12837916709551257390
 #define FYC	1.77245385090551602729
 #define ONEEPS	1.0000001
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
 	(void) P;
 	if ((xy.y = 1. - sin(lp.phi)) <= 0.)
 		xy.y = 0.;
@@ -14,7 +18,10 @@ FORWARD(s_forward); /* spheroid */
 	xy.y = FYC * (1. - xy.y);
 	return (xy);
 }
-INVERSE(s_inverse); /* spheroid */
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
 	lp.phi = xy.y / FYC - 1.;
 	if (fabs(lp.phi = 1. - lp.phi * lp.phi) < 1.)
 		lp.phi = asin(lp.phi);
@@ -26,5 +33,66 @@ INVERSE(s_inverse); /* spheroid */
 		lp.lam = xy.x / (FXC * sqrt(lp.lam));
 	return (lp);
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(collg) 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(collg) {
+    P->es = 0.0;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_collg_selftest (void) {return 0;}
+#else
+
+int pj_collg_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=collg   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+    XY s_fwd_expect[] = {
+        {249872.921577929839,  99423.1747884602082},
+        {254272.532301245432,  -98559.3077607425657},
+        {-249872.921577929839,  99423.1747884602082},
+        {-254272.532301245432,  -98559.3077607425657},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        {0.00158679719207879865,  0.00101017310941749921},
+        {0.001586769215623956,  -0.00101018201458258111},
+        {-0.00158679719207879865,  0.00101017310941749921},
+        {-0.001586769215623956,  -0.00101018201458258111},
+    };
+
+    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_comill.c b/src/PJ_comill.c
index ad9914d4..1d88dbd5 100644
--- a/src/PJ_comill.c
+++ b/src/PJ_comill.c
@@ -1,13 +1,14 @@
 /*
-The Compact Miller projection was designed by Tom Patterson, US National 
-Park Service, in 2014. The polynomial equation was developed by Bojan 
-Savric and Bernhard Jenny, College of Earth, Ocean, and Atmospheric 
+The Compact Miller projection was designed by Tom Patterson, US National
+Park Service, in 2014. The polynomial equation was developed by Bojan
+Savric and Bernhard Jenny, College of Earth, Ocean, and Atmospheric
 Sciences, Oregon State University.
 Port to PROJ.4 by Bojan Savric, 4 April 2016
 */
 
 #define PJ_LIB__
-#include	<projects.h>
+#include    <projects.h>
+
 PROJ_HEAD(comill, "Compact Miller") "\n\tCyl., Sph.";
 
 #define K1 0.9902
@@ -19,41 +20,108 @@ PROJ_HEAD(comill, "Compact Miller") "\n\tCyl., Sph.";
 #define EPS 1e-11
 #define MAX_Y (0.6000207669862655 * PI)
 
-FORWARD(s_forward); /* spheroid */
-	double lat_sq;
 
-	lat_sq = lp.phi * lp.phi;
-	xy.x = lp.lam;
-	xy.y = lp.phi * (K1 + lat_sq * (K2 + K3 * lat_sq));
-    return (xy);
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    double lat_sq;
+
+    lat_sq = lp.phi * lp.phi;
+    xy.x = lp.lam;
+    xy.y = lp.phi * (K1 + lat_sq * (K2 + K3 * lat_sq));
+    return xy;
 }
-INVERSE(s_inverse); /* spheroid */
-	double yc, tol, y2, y4, f, fder;
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    double yc, tol, y2, y4, f, fder;
 
     /* make sure y is inside valid range */
-	if (xy.y > MAX_Y) {
-		xy.y = MAX_Y;
-	} else if (xy.y < -MAX_Y) {
-		xy.y = -MAX_Y;
-	}
+    if (xy.y > MAX_Y) {
+        xy.y = MAX_Y;
+    } else if (xy.y < -MAX_Y) {
+        xy.y = -MAX_Y;
+    }
 
     /* latitude */
-	yc = xy.y;
+    yc = xy.y;
     for (;;) { /* Newton-Raphson */
-		y2 = yc * yc;
-		f = (yc * (K1 + y2 * (K2 + K3 * y2))) - xy.y;
-		fder = C1 + y2 * (C2 + C3 * y2);
-		yc -= tol = f / fder;
-		if (fabs(tol) < EPS) {
-			break;
-		}
-	}
-	lp.phi = yc;
+        y2 = yc * yc;
+        f = (yc * (K1 + y2 * (K2 + K3 * y2))) - xy.y;
+        fder = C1 + y2 * (C2 + C3 * y2);
+        yc -= tol = f / fder;
+        if (fabs(tol) < EPS) {
+            break;
+        }
+    }
+    lp.phi = yc;
 
     /* longitude */
-	lp.lam = xy.x;
+    lp.lam = xy.x;
+
+    return lp;
+}
+
+static void *freeup_new (PJ *P) {                       /* Destructor */
+    return pj_dealloc(P);
+}
+
+static void freeup (PJ *P) {
+    freeup_new (P);
+    return;
+}
+
+
+PJ *PROJECTION(comill) {
+    P->es = 0;
+
+    P->inv = s_inverse;
+    P->fwd = s_forward;
 
-	return (lp);
+    return P;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(comill) P->es = 0; P->inv = s_inverse; P->fwd = s_forward; ENDENTRY(P)
+
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_comill_selftest (void) {return 0;}
+#else
+
+int pj_comill_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char s_args[] = {"+proj=comill   +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.144255274179,  110611.859089458536},
+        {223402.144255274179,  -110611.859089458536},
+        {-223402.144255274179,  110611.859089458536},
+        {-223402.144255274179,  -110611.859089458536},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        {0.00179049310978382265,  0.000904106801510605831},
+        {0.00179049310978382265,  -0.000904106801510605831},
+        {-0.00179049310978382265,  0.000904106801510605831},
+        {-0.00179049310978382265,  -0.000904106801510605831},
+    };
+
+    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_crast.c b/src/PJ_crast.c
index 3f251ac6..0c5b706a 100644
--- a/src/PJ_crast.c
+++ b/src/PJ_crast.c
@@ -1,24 +1,95 @@
 #define PJ_LIB__
-# include	<projects.h>
-PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)")
-"\n\tPCyl., Sph.";
-#define XM	0.97720502380583984317
-#define RXM	1.02332670794648848847
-#define YM	3.06998012383946546542
-#define RYM	0.32573500793527994772
-#define THIRD	0.333333333333333333
-FORWARD(s_forward); /* spheroid */
-	(void) P;
-	lp.phi *= THIRD;
-	xy.x = XM * lp.lam * (2. * cos(lp.phi + lp.phi) - 1.);
-	xy.y = YM * sin(lp.phi);
-	return (xy);
+# include   <projects.h>
+
+PROJ_HEAD(crast, "Craster Parabolic (Putnins P4)") "\n\tPCyl., Sph.";
+
+#define XM  0.97720502380583984317
+#define RXM 1.02332670794648848847
+#define YM  3.06998012383946546542
+#define RYM 0.32573500793527994772
+#define THIRD   0.333333333333333333
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    (void) P;
+    lp.phi *= THIRD;
+    xy.x = XM * lp.lam * (2. * cos(lp.phi + lp.phi) - 1.);
+    xy.y = YM * sin(lp.phi);
+    return xy;
 }
-INVERSE(s_inverse); /* spheroid */
-	(void) P;
-	lp.phi = 3. * asin(xy.y * RYM);
-	lp.lam = xy.x * RXM / (2. * cos((lp.phi + lp.phi) * THIRD) - 1);
-	return (lp);
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    (void) P;
+    lp.phi = 3. * asin(xy.y * RYM);
+    lp.lam = xy.x * RXM / (2. * cos((lp.phi + lp.phi) * THIRD) - 1);
+    return lp;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(crast) 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(crast) {
+    P->es = 0.0;
+    P->inv = s_inverse;
+    P->fwd = s_forward;
+
+    return P;
+}
+
+#ifdef PJ_OMIT_SELFTEST
+int pj_crast_selftest (void) {return 0;}
+#else
+
+int pj_crast_selftest (void) {
+    double tolerance_lp = 1e-10;
+    double tolerance_xy = 1e-7;
+
+    char e_args[] = {"+proj=crast   +ellps=GRS80  +lat_1=0.5 +lat_2=2"};
+    char s_args[] = {"+proj=crast   +a=6400000    +lat_1=0.5 +lat_2=2"};
+
+    LP fwd_in[] = {
+        { 2, 1},
+        { 2,-1},
+        {-2, 1},
+        {-2,-1}
+    };
+
+
+    XY s_fwd_expect[] = {
+        {218280.142056780722,  114306.045604279774},
+        {218280.142056780722,  -114306.045604279774},
+        {-218280.142056780722,  114306.045604279774},
+        {-218280.142056780722,  -114306.045604279774},
+    };
+
+    XY inv_in[] = {
+        { 200, 100},
+        { 200,-100},
+        {-200, 100},
+        {-200,-100}
+    };
+
+    LP s_inv_expect[] = {
+        {0.00183225941982580187,  0.00087483943098902331},
+        {0.00183225941982580187,  -0.00087483943098902331},
+        {-0.00183225941982580187,  0.00087483943098902331},
+        {-0.00183225941982580187,  -0.00087483943098902331},
+    };
+
+    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