First steps toward simplified macros/internals

The brief version::

In an attempt to make proj.4 code slightly more secure and much easier
to read and maintain, I'm trying to eliminate a few unfortunate design
decisions from the early days of proj.4

The work will be *very* intrusive, especially in the PJ_xxx segment of
the code tree, but great care has been taken to design a process that
can be implemented stepwise and localized, one projection at a time,
then finalized with a relatively small and concentrated work package.

The (very) long version: See the comments in PJ_minimal.c

7 files changed, 519 insertions, 139 deletions

diff --git a/src/PJ_cass.c b/src/PJ_cass.c
index 38fa9db5..5d7aea3c 100644
--- a/src/PJ_cass.c
+++ b/src/PJ_cass.c
@@ -1,79 +1,126 @@
-#define PROJ_PARMS__ \
-	double m0; \
-	double n; \
-	double t; \
-	double a1; \
-	double c; \
-	double r; \
-	double dd; \
-	double d2; \
-	double a2; \
-	double tn; \
-	double *en;
 #define PJ_LIB__
 # include	<projects.h>
 PROJ_HEAD(cass, "Cassini") "\n\tCyl, Sph&Ell";
+
+
 # define EPS10	1e-10
 # define C1	.16666666666666666666
 # define C2	.00833333333333333333
 # define C3	.04166666666666666666
 # define C4	.33333333333333333333
 # define C5	.06666666666666666666
-FORWARD(e_forward); /* ellipsoid */
-	xy.y = pj_mlfn(lp.phi, P->n = sin(lp.phi), P->c = cos(lp.phi), P->en);
-	P->n = 1./sqrt(1. - P->es * P->n * P->n);
-	P->tn = tan(lp.phi); P->t = P->tn * P->tn;
-	P->a1 = lp.lam * P->c;
-	P->c *= P->es * P->c / (1 - P->es);
-	P->a2 = P->a1 * P->a1;
-	xy.x = P->n * P->a1 * (1. - P->a2 * P->t *
-		(C1 - (8. - P->t + 8. * P->c) * P->a2 * C2));
-	xy.y -= P->m0 - P->n * P->tn * P->a2 *
-		(.5 + (5. - P->t + 6. * P->c) * P->a2 * C3);
-	return (xy);
+
+
+struct opaque {
+    /* These are the only opaque members actually initialized */
+	double *en;
+	double m0;
+
+    /* Apparently this group of opaque members should be demoted to local variables in e_forward and e_inverse */
+	double n;
+	double t;
+	double a1;
+	double c;
+	double r;
+	double dd;
+	double d2;
+	double a2;
+	double tn;
+};
+
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    struct opaque *O = P->opaq;
+	xy.y = pj_mlfn(lp.phi, O->n = sin(lp.phi), O->c = cos(lp.phi), O->en);
+	O->n = 1./sqrt(1. - P->es * O->n * O->n);
+	O->tn = tan(lp.phi); O->t = O->tn * O->tn;
+	O->a1 = lp.lam * O->c;
+	O->c *= P->es * O->c / (1 - P->es);
+	O->a2 = O->a1 * O->a1;
+	xy.x = O->n * O->a1 * (1. - O->a2 * O->t *
+		(C1 - (8. - O->t + 8. * O->c) * O->a2 * C2));
+	xy.y -= O->m0 - O->n * O->tn * O->a2 *
+		(.5 + (5. - O->t + 6. * O->c) * O->a2 * C3);
+	return xy;
 }
-FORWARD(s_forward); /* spheroid */
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
 	xy.x = asin(cos(lp.phi) * sin(lp.lam));
 	xy.y = atan2(tan(lp.phi) , cos(lp.lam)) - P->phi0;
-	return (xy);
+	return xy;
 }
-INVERSE(e_inverse); /* ellipsoid */
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    struct opaque *O = P->opaq;
 	double ph1;
 
-	ph1 = pj_inv_mlfn(P->ctx, P->m0 + xy.y, P->es, P->en);
-	P->tn = tan(ph1); P->t = P->tn * P->tn;
-	P->n = sin(ph1);
-	P->r = 1. / (1. - P->es * P->n * P->n);
-	P->n = sqrt(P->r);
-	P->r *= (1. - P->es) * P->n;
-	P->dd = xy.x / P->n;
-	P->d2 = P->dd * P->dd;
-	lp.phi = ph1 - (P->n * P->tn / P->r) * P->d2 *
-		(.5 - (1. + 3. * P->t) * P->d2 * C3);
-	lp.lam = P->dd * (1. + P->t * P->d2 *
-		(-C4 + (1. + 3. * P->t) * P->d2 * C5)) / cos(ph1);
-	return (lp);
+	ph1 = pj_inv_mlfn(P->ctx, O->m0 + xy.y, P->es, O->en);
+	O->tn = tan(ph1); O->t = O->tn * O->tn;
+	O->n = sin(ph1);
+	O->r = 1. / (1. - P->es * O->n * O->n);
+	O->n = sqrt(O->r);
+	O->r *= (1. - P->es) * O->n;
+	O->dd = xy.x / O->n;
+	O->d2 = O->dd * O->dd;
+	lp.phi = ph1 - (O->n * O->tn / O->r) * O->d2 *
+		(.5 - (1. + 3. * O->t) * O->d2 * C3);
+	lp.lam = O->dd * (1. + O->t * O->d2 *
+		(-C4 + (1. + 3. * O->t) * O->d2 * C5)) / cos(ph1);
+	return lp;
 }
-INVERSE(s_inverse); /* spheroid */
-	lp.phi = asin(sin(P->dd = xy.y + P->phi0) * cos(xy.x));
-	lp.lam = atan2(tan(xy.x), cos(P->dd));
-	return (lp);
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+	lp.phi = asin(sin(P->opaq->dd = xy.y + P->phi0) * cos(xy.x));
+	lp.lam = atan2(tan(xy.x), cos(P->opaq->dd));
+	return lp;
 }
-FREEUP;
-	if (P) {
-		if (P->en)
-			pj_dalloc(P->en);
-		pj_dalloc(P);
-	}
+
+
+static void freeup(PJ *P) {                             /* Destructor */
+	if (0==P)
+        return;
+    if (0==P->opaq) {
+        pj_dealloc (P);
+        return;
+    }
+	pj_dealloc(P->opaq->en);
+    pj_dealloc(P->opaq);
+    pj_dealloc(P);
 }
-ENTRY1(cass, en)
-	if (P->es) {
-		if (!(P->en = pj_enfn(P->es))) E_ERROR_0;
-		P->m0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->en);
-		P->inv = e_inverse;
-		P->fwd = e_forward;
-	} else {
+
+
+PJ *PROJECTION(cass) {
+
+    /* Spheroidal? */
+	if (0==P->es) {
 		P->inv = s_inverse;
 		P->fwd = s_forward;
-	}
-ENDENTRY(P)
+    }
+
+    /* Otherwise it's ellipsoidal */
+    P->opaq = pj_calloc (1, sizeof (struct opaque));
+    if (0==P->opaq) {
+        freeup (P);
+        return 0;
+    }
+
+    P->opaq->en = pj_enfn(P->es);
+	if (0==P->opaq->en) {
+        freeup (P);
+        return 0;
+    }
+
+	P->opaq->m0 = pj_mlfn(P->phi0, sin(P->phi0), cos(P->phi0), P->opaq->en);
+	P->inv = e_inverse;
+	P->fwd = e_forward;
+
+    return P;
+}
diff --git a/src/PJ_merc.c b/src/PJ_merc.c
index 4b991c1c..89913e50 100644
--- a/src/PJ_merc.c
+++ b/src/PJ_merc.c
@@ -1,31 +1,53 @@
 #define PJ_LIB__
 #include	<projects.h>
 PROJ_HEAD(merc, "Mercator") "\n\tCyl, Sph&Ell\n\tlat_ts=";
+
 #define EPS10 1.e-10
-FORWARD(e_forward); /* ellipsoid */
-	if (fabs(fabs(lp.phi) - HALFPI) <= EPS10) F_ERROR;
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+	if (fabs(fabs(lp.phi) - HALFPI) <= EPS10)
+        F_ERROR;
 	xy.x = P->k0 * lp.lam;
 	xy.y = - P->k0 * log(pj_tsfn(lp.phi, sin(lp.phi), P->e));
-	return (xy);
+	return xy;
 }
-FORWARD(s_forward); /* spheroid */
-	if (fabs(fabs(lp.phi) - HALFPI) <= EPS10) F_ERROR;
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+	if (fabs(fabs(lp.phi) - HALFPI) <= EPS10)
+        F_ERROR;
 	xy.x = P->k0 * lp.lam;
 	xy.y = P->k0 * log(tan(FORTPI + .5 * lp.phi));
-	return (xy);
+	return xy;
 }
-INVERSE(e_inverse); /* ellipsoid */
-	if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL) I_ERROR;
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+	if ((lp.phi = pj_phi2(P->ctx, exp(- xy.y / P->k0), P->e)) == HUGE_VAL)
+        I_ERROR;
 	lp.lam = xy.x / P->k0;
-	return (lp);
+	return lp;
 }
-INVERSE(s_inverse); /* spheroid */
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
 	lp.phi = HALFPI - 2. * atan(exp(-xy.y / P->k0));
 	lp.lam = xy.x / P->k0;
-	return (lp);
+	return lp;
 }
-FREEUP; if (P) pj_dalloc(P); }
-ENTRY0(merc)
+
+
+static void freeup(PJ *P) {                             /* Destructor */
+    pj_dealloc(P);
+}
+
+
+PJ *PROJECTION(merc) {
 	double phits=0.0;
 	int is_phits;
 
@@ -33,15 +55,20 @@ ENTRY0(merc)
 		phits = fabs(pj_param(P->ctx, P->params, "rlat_ts").f);
 		if (phits >= HALFPI) E_ERROR(-24);
 	}
+
 	if (P->es) { /* ellipsoid */
 		if (is_phits)
 			P->k0 = pj_msfn(sin(phits), cos(phits), P->es);
 		P->inv = e_inverse;
 		P->fwd = e_forward;
-	} else { /* sphere */
+	}
+
+    else { /* sphere */
 		if (is_phits)
 			P->k0 = cos(phits);
 		P->inv = s_inverse;
 		P->fwd = s_forward;
 	}
-ENDENTRY(P)
+
+    return P;
+}
diff --git a/src/PJ_minimal.c b/src/PJ_minimal.c
new file mode 100644
index 00000000..210d25d3
--- /dev/null
+++ b/src/PJ_minimal.c
@@ -0,0 +1,203 @@
+/***********************************************************************
+
+   A minimal example of a new proj.4 projection implementation
+
+   ...and a verbose justification for some highly intrusive code
+   surgery
+
+************************************************************************
+
+**The brief version:**
+
+In an attempt to make proj.4 code slightly more secure and much easier
+to read and maintain, I'm trying to eliminate a few unfortunate design
+decisions from the early days of proj.4
+
+The work will be *very* intrusive, especially in the PJ_xxx segment of
+the code tree, but great care has been taken to design a process that
+can be implemented stepwise and localized, one projection at a time,
+then finalized with a relatively small and concentrated work package.
+
+**The (very) long version:**
+
+Gerald I. Evenden's original design for the proj.4 projection system
+is a beautiful example of software architecture, where a very limited
+set of policy rules leads to a well defined hierarchical structure and
+a high degree of both encapsulation and internal interoperability.
+
+In the proj.4 code, the policy rules are *enforced* by a system of
+preprocessor macros for building the scaffolding for implementation
+of a new projection.
+
+While this system of macros undeniably possesses the property of both
+reducing repetitive code and enforcing policy, unfortunately it also
+possesses two much less desirable properties:
+
+First, while enforcing policy, it also *hides* policy: The "beauty in
+simplicity" of Gerald's design is hidden behind layers of macros,
+whose architectural clarity do not match that of proj.4 in general.
+
+Second (and related), the macros make the source code look like
+something only vaguely related to C, making it hard to read (an effect
+that gets amplified to the tune of syntax highlighters getting confused
+by the macros).
+
+While the policy rule enforcement macros can be eliminated in relatively
+non-intrusive ways, a more fundamental flaw in the proj.4 use of macros
+is found in the PJ_xxx.c files implementing the individual projections:
+The use of internal redefinition of PJ, the fundamental proj data object,
+through the use of the PROJ_PARMS__ macro, makes the sizeof (PJ)
+fundamentally unknown to the calling pj_init function.
+
+This leads to code that is probably not in full conformance with the
+C standard.
+
+It is also a memory management catastrophe waiting to happen.
+
+But first and foremost, it leads to some very clumsy initialization code,
+where pj_init (the constructor function), needs to start the constsruction
+process by asking the PJ_xxx function to do the memory allocation (because
+pj_init does not know the size of the PROJ_PARMS-mangled PJ object being
+instantiated).
+
+Then, after doing some initialization work, pj_init returns control to
+PJ_xxx, asking it to finalize the initialization with the projection
+specific parameters specified by the PROJ_PARMS__ macro.
+
+Behind the scenes, hidden by two layers of macros, what happens is even
+worse, as a lot of the initialization code is duplicated in every PJ_xxx
+file, rather than being centralized in the pj_init function.
+
+**Solution procedure:**
+
+Evidently, the way to eliminate this clumsyness will be to introduce an
+opaque object, that is managed by tne individual PJ_xxx projection code,
+and represented as a simple void-pointer in the PJ object.
+
+This can be done one projection code file at a time, working through the
+code base as time permits (it will take at least a month).
+
+When a PJ_xxx file is on the surgical bench, it will also have its
+ENTRYA/ENTRY0/ENTRY1/ENTRY2/ENDENTRY/etc. etc. macro-guts torn out and
+replaced by the PROJECTION macro (introduced in projects.h).
+
+This leads to code that looks a lot more like real C, and hence is much
+less confusing to both syntax higlighters and humans. It also leads
+to code that, after all projections have been processed, with a final
+sweep over the code base can be brought into the style of the code in
+PJ_minimal.c
+
+In my humble opinion the result wil be a code base that is not only easier
+to maintain, but also more welcoming to new contributors.
+
+And if proj is to expand its strong basis in projections into the fields
+of geodetic transformations and general geometric geodesy, we will need
+to be able to attract quite a few expert geodesist contributors.
+
+And since expert geodesists are not necessarily expert coders, a welcoming
+code base is a real asset (to put the icing on the cake of the already
+welcoming user- and developer community).
+
+Note that the entire process does not touch the algorithmic/mathematical
+parts of the code at all - it is actuallly an attempt to make this part
+stand out more clearly.
+
+---
+
+The attached material is an attempt to show what happens if we remove
+the layers of macros, and introduce a more centralized approach to
+memory allocation and initialization.
+
+Please note, however, that the level of cantralization achieved here
+is not yet fully supported by the proj.4 infrastructure: It is an
+example, intended to show what can be achieved through a smooth,
+gradual and safe refactoring of the existing layered macro system.
+
+In my humble opinion, this version makes the beauty of Gerald's design
+much more evident than the current layered-macro-version.
+
+Thomas Knudsen, thokn@sdfe.dk, 2016-03-31
+
+***********************************************************************/
+
+#define PJ_LIB__
+#include	<projects.h>
+#include <assert.h>
+PROJ_HEAD(minimal, "Minimal example (brief description goes here)");
+
+
+/* Projection specific elements for the PJ object */
+struct opaque {
+	double a;
+	int b;
+};
+
+
+static XY e_forward (LP lp, PJ *P) {          /* Ellipsoidal, forward */
+    XY xy = {0.0,0.0};
+    /* Actual ellipsoidal forward code goes here */
+    xy.y = lp.lam + P->es;
+    xy.x = lp.phi + 42;
+	return xy;
+}
+
+
+static XY s_forward (LP lp, PJ *P) {           /* Spheroidal, forward */
+    XY xy = {0.0,0.0};
+    /* Actual spheroidal forward code goes here */
+    xy.y = lp.lam + P->es;
+    xy.x = lp.phi + 42;
+	return xy;
+}
+
+
+static LP e_inverse (XY xy, PJ *P) {          /* Ellipsoidal, inverse */
+    LP lp = {0.0,0.0};
+    /* Actual ellipsoidal forward code goes here */
+    lp.lam = xy.x - P->es;
+    lp.phi = xy.y - P->opaq->b;
+	return lp;
+}
+
+
+static LP s_inverse (XY xy, PJ *P) {           /* Spheroidal, inverse */
+    LP lp = {0.0,0.0};
+    /* Actual spheroidal forward code goes here */
+    lp.lam = xy.x - P->es;
+    lp.phi = xy.y - P->opaq->b;
+	return lp;
+}
+
+
+static void freeup(PJ *P) {                                    /* Destructor */
+    if (P==0)
+        return;
+    /* Projection specific deallocation goes here */
+    pj_dealloc (P->opaq);
+    pj_dealloc (P);
+    return;
+}
+
+
+PJ *pj_projection_specific_setup_minimal (PJ *P) {
+    pj_prepare (P, des_minimal, freeup, sizeof (struct opaque));
+    if (0==P->opaq) {
+        freeup (P);
+        return 0;
+    }
+
+    P->opaq->a = 42.42;
+    P->opaq->b = 42;
+
+    /* Spheroidal? */
+	if (0==P->es) {
+		P->fwd = s_forward;
+		P->inv = s_inverse;
+    }
+
+    /* Otherwise it's ellipsoidal */
+    P->fwd = e_forward;
+	P->inv = e_inverse;
+
+    return P;
+}
diff --git a/src/pj_init.c b/src/pj_init.c
index 13d469da..ca7b9e4a 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -110,7 +110,7 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
     *sword = 't';
 
     /* loop till we find our target keyword */
-    while (*next_char) 
+    while (*next_char)
     {
         next_char = fill_buffer(state, next_char);
 
@@ -119,9 +119,9 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
             next_char++;
 
         next_char = fill_buffer(state, next_char);
-        
+
         /* for comments, skip past end of line. */
-        if( *next_char == '#' ) 
+        if( *next_char == '#' )
         {
             while( *next_char && *next_char != '\n' )
                 next_char++;
@@ -131,11 +131,11 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
                 next_char++;
             if (*next_char == '\r')
                 next_char++;
-            
-        } 
+
+        }
 
         /* Is this our target? */
-        else if( *next_char == '<' ) 
+        else if( *next_char == '<' )
         {
             /* terminate processing target on the next block definition */
             if (in_target)
@@ -143,7 +143,7 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
 
             next_char++;
             if (strncmp(name, next_char, len) == 0
-                && next_char[len] == '>') 
+                && next_char[len] == '>')
             {
                 /* skip past target word */
                 next_char += len + 1;
@@ -151,14 +151,14 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
                 if(found_def)
                     *found_def = 1;
             }
-            else 
+            else
             {
                 /* skip past end of line */
                 while( *next_char && *next_char != '\n' )
                     next_char++;
             }
         }
-        else if (in_target) 
+        else if (in_target)
         {
             const char *start_of_word = next_char;
             int word_len = 0;
@@ -183,27 +183,27 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
             if (!pj_param(ctx, *start, sword).i) {
                 /* don't default ellipse if datum, ellps or any earth model
                    information is set. */
-                if( strncmp(sword+1,"ellps=",6) != 0 
-                    || (!pj_param(ctx, *start, "tdatum").i 
-                        && !pj_param(ctx, *start, "tellps").i 
-                        && !pj_param(ctx, *start, "ta").i 
-                        && !pj_param(ctx, *start, "tb").i 
-                        && !pj_param(ctx, *start, "trf").i 
+                if( strncmp(sword+1,"ellps=",6) != 0
+                    || (!pj_param(ctx, *start, "tdatum").i
+                        && !pj_param(ctx, *start, "tellps").i
+                        && !pj_param(ctx, *start, "ta").i
+                        && !pj_param(ctx, *start, "tb").i
+                        && !pj_param(ctx, *start, "trf").i
                         && !pj_param(ctx, *start, "tf").i) )
                 {
                     next = next->next = pj_mkparam(sword+1);
                 }
             }
-            
+
         }
-        else 
+        else
         {
             /* skip past word */
             while( *next_char && !isspace(*next_char) )
                 next_char++;
-            
+
         }
-    }        
+    }
 
     if (errno == 25)
         errno = 0;
@@ -229,7 +229,7 @@ get_defaults(projCtx ctx, paralist **start, paralist *next, char *name) {
     if (errno)
         errno = 0; /* don't care if can't open file */
     ctx->last_errno = 0;
-    
+
     return next;
 }
 
@@ -245,11 +245,11 @@ get_init(projCtx ctx, paralist **start, paralist *next, char *name,
     const paralist *orig_next = next;
 
     (void)strncpy(fname, name, MAX_PATH_FILENAME + ID_TAG_MAX + 1);
-	
-    /* 
-    ** Search for file/key pair in cache 
+
+    /*
+    ** Search for file/key pair in cache
     */
-	
+
     init_items = pj_search_initcache( name );
     if( init_items != NULL )
     {
@@ -275,8 +275,8 @@ get_init(projCtx ctx, paralist **start, paralist *next, char *name,
     if (errno == 25)
         errno = 0; /* unknown problem with some sys errno<-25 */
 
-    /* 
-    ** If we seem to have gotten a result, insert it into the 
+    /*
+    ** If we seem to have gotten a result, insert it into the
     ** init file cache.
     */
     if( next != NULL && next != orig_next )
@@ -308,7 +308,7 @@ pj_init_plus_ctx( projCtx ctx, const char *definition )
     char	*defn_copy;
     int		argc = 0, i, blank_count = 0;
     PJ	    *result = NULL;
-    
+
     /* make a copy that we can manipulate */
     defn_copy = (char *) pj_malloc( strlen(definition)+1 );
     strcpy( defn_copy, definition );
@@ -328,13 +328,13 @@ pj_init_plus_ctx( projCtx ctx, const char *definition )
                     defn_copy[i - blank_count] = '\0';
                     blank_count = 0;
                 }
-                
+
                 if( argc+1 == MAX_ARG )
                 {
                     pj_ctx_set_errno( ctx, -44 );
                     goto bum_call;
                 }
-                
+
                 argv[argc++] = defn_copy + i + 1;
             }
             break;
@@ -455,7 +455,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     PIN->rone_es = 1./PIN->one_es;
 
     /* Now that we have ellipse information check for WGS84 datum */
-    if( PIN->datum_type == PJD_3PARAM 
+    if( PIN->datum_type == PJD_3PARAM
         && PIN->datum_params[0] == 0.0
         && PIN->datum_params[1] == 0.0
         && PIN->datum_params[2] == 0.0
@@ -464,7 +464,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     {
         PIN->datum_type = PJD_WGS84;
     }
-        
+
     /* set PIN->geoc coordinate system */
     PIN->geoc = (PIN->es && pj_param(ctx, start, "bgeoc").i);
 
@@ -532,7 +532,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     /* set units */
     s = 0;
-    if ((name = pj_param(ctx, start, "sunits").s) != NULL) { 
+    if ((name = pj_param(ctx, start, "sunits").s) != NULL) {
         for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i) ;
         if (!s) { pj_ctx_set_errno( ctx, -7 ); goto bum_call; }
         s = pj_units[i].to_meter;
@@ -547,7 +547,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     /* set vertical units */
     s = 0;
-    if ((name = pj_param(ctx, start, "svunits").s) != NULL) { 
+    if ((name = pj_param(ctx, start, "svunits").s) != NULL) {
         for (i = 0; (s = pj_units[i].id) && strcmp(name, s) ; ++i) ;
         if (!s) { pj_ctx_set_errno( ctx, -7 ); goto bum_call; }
         s = pj_units[i].to_meter;
@@ -564,7 +564,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     /* prime meridian */
     s = 0;
-    if ((name = pj_param(ctx, start, "spm").s) != NULL) { 
+    if ((name = pj_param(ctx, start, "spm").s) != NULL) {
         const char *value = NULL;
         char *next_str = NULL;
 
@@ -576,8 +576,8 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
                 break;
             }
         }
-            
-        if( value == NULL 
+
+        if( value == NULL
             && (dmstor_ctx(ctx,name,&next_str) != 0.0  || *name == '0')
             && *next_str == '\0' )
             value = name;
@@ -639,3 +639,27 @@ pj_free(PJ *P) {
         P->pfree(P);
     }
 }
+
+
+
+
+
+
+
+
+/************************************************************************/
+/*                              pj_prepare()                            */
+/*                                                                      */
+/*      Helper function for the PJ_xxxx functions providing the         */
+/*      projection specific setup for each projection type.             */
+/*                                                                      */
+/*      Currently not used, but placed here as part of the  material    */
+/*      Demonstrating the idea for a future PJ_xxx architecture         */
+/*      (cf. pj_minimal.c)                                              */
+/*                                                                      */
+/************************************************************************/
+void pj_prepare (PJ *P, const char *description, void (*freeup)(struct PJconsts *), size_t sizeof_struct_opaque) {
+    P->descr = description;
+    P->pfree = freeup;
+    P->opaq  = pj_calloc (1, sizeof_struct_opaque);
+}
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index 80443a2b..aab69e99 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -9,19 +9,64 @@
 pj_malloc(size_t size) {
 /*
 / Currently, pj_malloc is a hack to solve an errno problem.
-/ The problem is described in more details at 
-/ https://bugzilla.redhat.com/bugzilla/show_bug.cgi?id=86420. 
-/ It seems, that pj_init and similar functions incorrectly 
-/ (under debian/glibs-2.3.2) assume that pj_malloc resets 
+/ The problem is described in more details at
+/ https://bugzilla.redhat.com/bugzilla/show_bug.cgi?id=86420.
+/ It seems, that pj_init and similar functions incorrectly
+/ (under debian/glibs-2.3.2) assume that pj_malloc resets
 / errno after success. pj_malloc tries to mimic this.
 */
         int old_errno = errno;
-        void *res = malloc(size);       
+        void *res = malloc(size);
         if ( res && !old_errno )
-                errno = 0;        	 
+                errno = 0;
         return res;
 }
 	void
 pj_dalloc(void *ptr) {
 	free(ptr);
 }
+
+
+/**********************************************************************/
+void *pj_calloc (size_t n, size_t size) {
+/***********************************************************************
+
+pj_calloc is the pj-equivalent of calloc().
+
+It allocates space for an array of <n> elements of size <size>.
+The array is initialized to zeros.
+
+***********************************************************************/
+    void *res = pj_malloc (n*size);
+    if (0==res)
+        return 0;
+    memset (res, 0, n*size);
+    return res;
+}
+
+
+/**********************************************************************/
+void *pj_dealloc (void *ptr) {
+/***********************************************************************
+
+pj_dealloc supports the common use case of "clean up and return a null
+pointer" to signal an error in a multi level allocation:
+
+    struct foo { int bar; int *baz; };
+
+    struct foo *p = pj_calloc (1, sizeof (struct foo));
+    if (0==p)
+        return 0;
+
+    p->baz = pj_calloc (10, sizeof(int));
+    if (0==p->baz)
+        return pj_dealloc (p); // clean up + signal error by 0-return
+
+    return p;  // success
+
+***********************************************************************/
+    if (0==ptr)
+        return 0;
+    pj_dalloc (ptr);
+    return 0;
+}
diff --git a/src/proj_api.h b/src/proj_api.h
index b482c0f5..dd25d61a 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -130,6 +130,8 @@ char *pj_get_def(projPJ, int);
 projPJ pj_latlong_from_proj( projPJ );
 void *pj_malloc(size_t);
 void pj_dalloc(void *);
+void *pj_calloc (size_t n, size_t size);
+void *pj_dealloc (void *ptr);
 char *pj_strerrno(int);
 int *pj_get_errno_ref(void);
 const char *pj_get_release(void);
diff --git a/src/projects.h b/src/projects.h
index 411a3403..3106eb76 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -48,7 +48,7 @@
 #define C_NAMESPACE extern "C"
 #define C_NAMESPACE_VAR extern "C"
 extern "C" {
-#else    
+#else
 #define C_NAMESPACE extern
 #define C_NAMESPACE_VAR
 #endif
@@ -132,8 +132,8 @@ typedef struct {
 
 /* datum_type values */
 #define PJD_UNKNOWN   0
-#define PJD_3PARAM    1   
-#define PJD_7PARAM    2   
+#define PJD_3PARAM    1
+#define PJD_7PARAM    2
 #define PJD_GRIDSHIFT 3
 #define PJD_WGS84     4   /* WGS84 (or anything considered equivelent) */
 
@@ -143,7 +143,7 @@ typedef struct {
 #define PJD_ERR_GRID_AREA           -48
 #define PJD_ERR_CATALOG             -49
 
-#define USE_PROJUV 
+#define USE_PROJUV
 
 typedef struct { double u, v; } projUV;
 typedef struct { double r, i; }	COMPLEX;
@@ -163,7 +163,7 @@ typedef struct { double lam, phi, z; } LPZ;
 
 typedef union { double  f; int  i; char *s; } PROJVALUE;
 struct PJconsts;
-    
+
 struct PJ_LIST {
 	char	*id;		/* projection keyword */
 	struct PJconsts	*(*proj)(struct PJconsts*);/* projection entry point */
@@ -197,14 +197,14 @@ typedef struct {
     double ll_long;      /* lower left corner coordinates (radians) */
     double ll_lat;
     double ur_long;      /* upper right corner coordinates (radians) */
-    double ur_lat; 
+    double ur_lat;
 } PJ_Region;
 
 struct DERIVS {
     double x_l, x_p; /* derivatives of x for lambda-phi */
     double y_l, y_p; /* derivatives of y for lambda-phi */
 };
-    
+
 struct FACTORS {
 	struct DERIVS der;
 	double h, k;	/* meridinal, parallel scales */
@@ -226,12 +226,18 @@ typedef struct ARG_list {
 	/* base projection data structure */
 
 
+#ifdef PJ_LIB__
+    /* we need this forward declaration in order to be able to add a
+       pointer to struct opaque to the typedef struct PJconsts below */
+    struct opaque;
+#endif
+
 typedef struct PJconsts {
     projCtx_t *ctx;
 	XY  (*fwd)(LP, struct PJconsts *);
 	LP  (*inv)(XY, struct PJconsts *);
         XYZ (*fwd3d)(LPZ, struct PJconsts *);
-        LPZ (*inv3d)(XYZ, struct PJconsts *); 
+        LPZ (*inv3d)(XYZ, struct PJconsts *);
 	void (*spc)(LP, struct PJconsts *, struct FACTORS *);
 	void (*pfree)(struct PJconsts *);
 	const char *descr;
@@ -253,7 +259,7 @@ typedef struct PJconsts {
 		x0, y0, /* easting and northing */
 		k0,	/* general scaling factor */
 		to_meter, fr_meter; /* cartesian scaling */
-    
+
         int     datum_type; /* PJD_UNKNOWN/3PARAM/7PARAM/GRIDSHIFT/WGS84 */
         double  datum_params[7];
         struct _pj_gi **gridlist;
@@ -272,9 +278,9 @@ typedef struct PJconsts {
         /* New Datum Shift Grid Catalogs */
         char   *catalog_name;
         struct _PJ_GridCatalog *catalog;
-    
+
         double   datum_date;
-    
+
         struct _pj_gi *last_before_grid;
         PJ_Region     last_before_region;
         double        last_before_date;
@@ -283,6 +289,10 @@ typedef struct PJconsts {
         PJ_Region     last_after_region;
         double        last_after_date;
 
+#ifdef PJ_LIB__
+        struct opaque *opaq;
+#endif
+
 #ifdef PROJ_PARMS__
 PROJ_PARMS__
 #endif /* end of optional extensions */
@@ -297,7 +307,7 @@ extern struct PJ_LIST pj_list[];
 #else
 #define PROJ_HEAD(id, name) \
     struct PJconsts *pj_##id(struct PJconsts*); extern char * const pj_s_##id;
-    
+
 #include "pj_list.h"
 #undef PROJ_HEAD
 #define PROJ_HEAD(id, name) {#id, pj_##id, &pj_s_##id},
@@ -350,7 +360,27 @@ extern struct PJ_PRIME_MERIDIANS pj_prime_meridians[];
 #define INVERSE3D(name) static LPZ name(XYZ xyz, PJ *P) {LPZ lpz = {0.0, 0.0, 0.0}
 #define FREEUP static void freeup(PJ *P) {
 #define SPECIAL(name) static void name(LP lp, PJ *P, struct FACTORS *fac)
+
+
+/* cleaned up alternative to most of the "repetitive projection code" macros */
+#define PROJECTION(name)                                     \
+pj_projection_specific_setup_##name (PJ *P);                 \
+C_NAMESPACE_VAR const char * const pj_s_##name = des_##name; \
+C_NAMESPACE PJ *pj_##name (PJ *P) {                          \
+    if (P)                                                   \
+        return pj_projection_specific_setup_##name (P);      \
+    P = (PJ*) pj_calloc (1, sizeof(PJ));                     \
+    if (0==P)                                                \
+        return 0;                                            \
+    P->pfree = freeup;                                       \
+    P->descr = des_##name;                                   \
+    return P;                                                \
+}                                                            \
+PJ *pj_projection_specific_setup_##name (PJ *P)
+
 #endif
+
+
 #define MAX_TAB_ID 80
 typedef struct { float lam, phi; } FLP;
 typedef struct { int lam, phi; } ILP;
@@ -366,8 +396,8 @@ struct CTABLE {
 typedef struct _pj_gi {
     char *gridname;   /* identifying name of grid, eg "conus" or ntv2_0.gsb */
     char *filename;   /* full path to filename */
-    
-    const char *format; /* format of this grid, ie "ctable", "ntv1", 
+
+    const char *format; /* format of this grid, ie "ctable", "ntv1",
                            "ntv2" or "missing". */
 
     int   grid_offset; /* offset in file, for delayed loading */
@@ -414,6 +444,7 @@ int pj_ell_set(projCtx ctx, paralist *, double *, double *);
 int pj_datum_set(projCtx,paralist *, PJ *);
 int pj_prime_meridian_set(paralist *, PJ *);
 int pj_angular_units_set(paralist *, PJ *);
+void pj_prepare (PJ *P, const char *description, void (*freeup)(struct PJconsts *), size_t sizeof_struct_opaque);
 
 paralist *pj_clone_paralist( const paralist* );
 paralist*pj_search_initcache( const char *filekey );
@@ -439,7 +470,7 @@ struct PW_COEF {/* row coefficient structure */
     int m;		/* number of c coefficients (=0 for none) */
     double *c;	/* power coefficients */
 };
- 
+
 /* Approximation structures and procedures */
 typedef struct {	/* Chebyshev or Power series structure */
 	projUV a, b;		/* power series range for evaluation */
@@ -457,6 +488,7 @@ void **vector2(int, int, int);
 void freev2(void **v, int nrows);
 int bchgen(projUV, projUV, int, int, projUV **, projUV(*)(projUV));
 int bch2bps(projUV, projUV, projUV **, int, int);
+
 /* nadcon related protos */
 LP nad_intr(LP, struct CTABLE *);
 LP nad_cvt(LP, int, struct CTABLE *);
@@ -470,15 +502,15 @@ void nad_free(struct CTABLE *);
 /* higher level handling of datum grid shift files */
 
 int pj_apply_vgridshift( PJ *defn, const char *listname,
-                         PJ_GRIDINFO ***gridlist_p, 
+                         PJ_GRIDINFO ***gridlist_p,
                          int *gridlist_count_p,
-                         int inverse, 
+                         int inverse,
                          long point_count, int point_offset,
                          double *x, double *y, double *z );
-int pj_apply_gridshift_2( PJ *defn, int inverse, 
+int pj_apply_gridshift_2( PJ *defn, int inverse,
                           long point_count, int point_offset,
                           double *x, double *y, double *z );
-int pj_apply_gridshift_3( projCtx ctx, 
+int pj_apply_gridshift_3( projCtx ctx,
                           PJ_GRIDINFO **gridlist, int gridlist_count,
                           int inverse, long point_count, int point_offset,
                           double *x, double *y, double *z );
@@ -493,15 +525,15 @@ void pj_gridinfo_free( projCtx, PJ_GRIDINFO * );
 PJ_GridCatalog *pj_gc_findcatalog( projCtx, const char * );
 PJ_GridCatalog *pj_gc_readcatalog( projCtx, const char * );
 void pj_gc_unloadall( projCtx );
-int pj_gc_apply_gridshift( PJ *defn, int inverse, 
+int pj_gc_apply_gridshift( PJ *defn, int inverse,
                            long point_count, int point_offset,
                            double *x, double *y, double *z );
-int pj_gc_apply_gridshift( PJ *defn, int inverse, 
+int pj_gc_apply_gridshift( PJ *defn, int inverse,
                            long point_count, int point_offset,
                            double *x, double *y, double *z );
 
-PJ_GRIDINFO *pj_gc_findgrid( projCtx ctx, 
-                             PJ_GridCatalog *catalog, int after, 
+PJ_GRIDINFO *pj_gc_findgrid( projCtx ctx,
+                             PJ_GridCatalog *catalog, int after,
                              LP location, double date,
                              PJ_Region *optional_region,
                              double *grid_date );