Replace pj_ell_set with reimplementation... (#675)

* Replace pj_ell_set with reimplementation supporting ellipsoid inheritance

* remove unreachable code from pj_ell_set.c

* Swap steps, so ellps args are read first, in accordance with historical behaviour

* Add ellipsoid tests to CI targets

* Reduce some optimistic tolerances

OS/X appears to have a slightly off float handling, resulting in differences at the nanometer level. Switching to 10 nm.

10 files changed, 534 insertions, 41 deletions

diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index d8e1a921..5becfd74 100644
--- a/src/PJ_healpix.c
+++ b/src/PJ_healpix.c
@@ -30,6 +30,7 @@
  *****************************************************************************/
 # define PJ_LIB__
 # include <errno.h>
+# include "proj_internal.h"
 # include <proj.h>
 # include "projects.h"
 
@@ -624,7 +625,7 @@ PJ *PROJECTION(healpix) {
             return destructor(P, ENOMEM);
         Q->qp = pj_qsfn(1.0, P->e, P->one_es);  /* For auth_lat(). */
         P->a = P->a*sqrt(0.5*Q->qp);            /* Set P->a to authalic radius. */
-        P->ra = 1.0/P->a;
+        pj_calc_ellps_params (P, P->a, P->es);  /* Ensure we have a consistent parameter set */
         P->fwd = e_healpix_forward;
         P->inv = e_healpix_inverse;
     } else {
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index e1568423..19b8e18d 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -317,7 +317,7 @@ static void set_ellipsoid(PJ *P) {
 
     /* Check if there's any ellipsoid specification in the global params. */
     /* If not, use WGS84 as default                                       */
-    if (pj_ell_set(P->ctx, P->params, &P->a, &P->es)) {
+    if (0 != pj_ellipsoid (P)) {
         P->a  = 6378137.0;
         P->es = .00669438002290341575;
     }
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 8a9b39a1..8ba928b7 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,13 +1,424 @@
 /* set ellipsoid parameters a and es */
 #include <string.h>
+#include <errno.h>
+#include <proj.h>
 #include "proj_internal.h"
 #include "projects.h"
+
+/* series coefficients for calculating ellipsoid-equivalent spheres */
 #define SIXTH .1666666666666666667  /* 1/6 */
 #define RA4   .04722222222222222222 /* 17/360 */
 #define RA6   .02215608465608465608 /* 67/3024 */
 #define RV4   .06944444444444444444 /* 5/72 */
 #define RV6   .04243827160493827160 /* 55/1296 */
 
+/* Prototypes of the pj_ellipsoid helper functions */
+static int pj_ellps_handler (PJ *P);
+static int pj_size (PJ *P);
+static int pj_shape (PJ *P);
+static int pj_spherification (PJ *P);
+
+static paralist *pj_get_param (paralist *list, char *key);
+static char     *pj_param_value (paralist *list);
+static PJ_ELLPS *pj_find_ellps (char *name);
+
+
+/***************************************************************************************/
+int pj_ellipsoid (PJ *P) {
+/****************************************************************************************
+    This is a replacement for the clasic PROJ pj_ell_set function. The main difference
+    is that pj_ellipsoid augments the PJ object with a copy of the exact tags used to
+    define its related ellipsoid.
+
+    This makes it possible to let a new PJ object inherit the geometrical properties
+    of an existing one.
+
+    A complete ellipsoid definition comprises a size (primary) and a shape (secondary)
+    parameter.
+
+    Size parameters supported are:
+        R, defining the radius of a spherical planet
+        a, defining the semimajor axis of an ellipsoidal planet
+
+    Shape parameters supported are:
+        rf, the reverse flattening of the ellipsoid
+        f,  the flattening of the ellipsoid
+        es, the eccentricity squared
+        e,  the eccentricity
+        b,  the semiminor axis
+
+    The ellps=xxx parameter provides both size and shape for a number of built in
+    ellipsoid definitions.
+
+    The ellipsoid definition may be augmented with a spherification flag, turning
+    the ellipsoid into a sphere with features defined by the ellipsoid.
+
+    Spherification parameters supported are:
+        R_A, which gives a sphere with the same surface area as the ellipsoid
+        R_A, which gives a sphere with the same volume as the ellipsoid
+
+        R_a, which gives a sphere with R = (a + b)/2   (arithmetic mean)
+        R_g, which gives a sphere with R = sqrt(a*b)   (geometric mean)
+        R_h, which gives a sphere with R = 2*a*b/(a+b) (harmonic mean)
+
+        R_lat_a=phi, which gives a sphere with R being the arithmetic mean of
+             of the corresponding ellipsoid at latitude phi.
+        R_lat_g=phi, which gives a sphere with R being the geometric mean of
+             of the corresponding ellipsoid at latitude phi.
+
+    If R is given as size parameter, any shape and spherification parameters
+    given are ignored.
+
+    If size and shape is given as ellps=xxx, later shape and size parameters
+    are are taken into account as modifiers for the built in ellipsoid definition.
+
+    While this may seem strange, it is in accordance with historical PROJ
+    behaviour. It can e.g. be used to define coordinates on the ellipsoid
+    scaled to unit semimajor axis by specifying "+ellps=xxx +a=1"
+
+****************************************************************************************/
+    int   last_errno;
+
+    last_errno = proj_errno_reset (P);
+
+    P->def_size = P->def_shape = P->def_spherification = P->def_ellps = 0;
+
+    /* If an ellps argument is specified, start by using that */
+    if (0 != pj_ellps_handler (P))
+        return proj_errno (P);
+
+    /* We may overwrite the size */
+    if (0 != pj_size (P))
+        return proj_errno (P);
+
+    /* We may also overwrite the shape */
+    if (0 != pj_shape (P))
+        return proj_errno (P);
+
+    /* When we're done with it, we compute all related ellipsoid parameters */
+    pj_calc_ellps_params (P, P->a, P->es);
+
+    /* And finally, we may turn it into a sphere */
+    if (0 != pj_spherification (P))
+        return proj_errno (P);
+
+    if (proj_errno (P))
+        return proj_errno (P);
+
+    /* success - restore previous error status */
+    return proj_errno_restore (P, last_errno);
+}
+
+
+/***************************************************************************************/
+static int pj_ellps_handler (PJ *P) {
+/***************************************************************************************/
+    PJ B;
+    PJ_ELLPS *ellps;
+    paralist *par = 0;
+    char *def, *name;
+
+    /* ellps specified? */
+    par = pj_get_param (P->params, "ellps");
+    if (0==par)
+        return 0;
+
+    /* This amounts to zeroing all ellipsoidal parameters in P */
+    memset (&B, 0, sizeof (B));
+    pj_inherit_ellipsoid_defs(&B, P);
+
+    /* move B into P's context */
+    B.ctx = P->ctx;
+
+    /* Store definition */
+    P->def_ellps = def = par->param;
+    par->used = 1;
+
+    if (strlen (def) < 7)
+        return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+    name = def + 6;
+    ellps = pj_find_ellps (name);
+    if (0==ellps)
+        return proj_errno_set (P, PJD_ERR_UNKNOWN_ELLP_PARAM);
+
+    /* Now, we have the right size and shape parameters and can produce */
+    /* a fake PJ to make pj_shape do the hard work for us */
+    B = *P;
+    B.params = pj_mkparam (ellps->major);
+    B.params->next = pj_mkparam (ellps->ell);
+
+    pj_size (&B);
+    pj_shape (&B);
+
+    P->a = B.a;
+    P->b = B.b;
+    P->f = B.f;
+    P->e = B.e;
+    P->es = B.es;
+
+    pj_dealloc (B.params->next);
+    pj_dealloc (B.params);
+    if (proj_errno (&B))
+        return proj_errno (&B);
+     return 0;
+}
+
+
+/***************************************************************************************/
+static int pj_size (PJ *P) {
+/***************************************************************************************/
+    char *keys[]  =  {"R", "a"};
+    paralist *par = 0;
+    size_t i, len;
+    len = sizeof (keys) / sizeof (char *);
+    /* Check which size key is specified */
+    for (i = 0;  i < len;  i++) {
+        par = pj_get_param (P->params, keys[i]);
+        if (par)
+            break;
+    }
+
+    /* A size parameter *must* be given, but may have been given as ellps prior */
+    if (0==par)
+        return 0!=P->a? 0: proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+
+    P->def_size = par->param;
+    par->used = 1;
+    P->a = pj_atof (pj_param_value (par));
+    if (P->a <= 0)
+        return proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+    if (HUGE_VAL==P->a)
+        return proj_errno_set (P, PJD_ERR_MAJOR_AXIS_NOT_GIVEN);
+
+    if ('R'==par->param[0]) {
+        P->es = P->f = P->e = P->rf = 0;
+        P->b = P->a;
+    }
+    return 0;
+}
+
+
+/***************************************************************************************/
+static int pj_shape (PJ *P) {
+/***************************************************************************************/
+    char *keys[]  = {"rf", "f", "es", "e", "b"};
+    paralist *par = 0;
+    char *def = 0;
+    size_t i, len;
+
+    par = 0;
+    len = sizeof (keys) / sizeof (char *);
+
+    /* Check which shape key is specified */
+    for (i = 0;  i < len;  i++) {
+        par = pj_get_param (P->params, keys[i]);
+        if (par)
+            break;
+    }
+
+    /* Not giving a shape parameter means selecting a sphere, unless shape */
+    /* has been selected previously via ellps=xxx */
+    if (0==par && P->es != 0)
+        return 0;
+    if (0==par && P->es==0) {
+        P->es = P->f = 0;
+        P->b = P->a;
+        return 0;
+    }
+
+    P->def_shape = def = par->param;
+    par->used = 1;
+    P->es = P->f = P->b = P->e = P->rf = 0;
+
+    switch (i) {
+
+    /* reverse flattening, rf */
+    case 0:
+         P->rf = pj_atof (pj_param_value (par));
+         if (HUGE_VAL==P->rf)
+             return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+         P->f = 1 / P->rf;
+         P->es = 2*P->f - P->f*P->f;
+         break;
+
+    /* flattening, f */
+    case 1:
+        P->f = pj_atof (pj_param_value (par));
+        if (HUGE_VAL==P->f)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        if (0==P->f)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        P->rf = 1 / P->f;
+        P->es = 2*P->f - P->f*P->f;
+        break;
+
+    /* eccentricity squared, es */
+    case 2:
+        P->es = pj_atof (pj_param_value (par));
+        if (HUGE_VAL==P->es)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        if (1==P->es)
+            return proj_errno_set (P, PJD_ERR_ECCENTRICITY_IS_ONE);
+        break;
+
+    /* eccentricity, e */
+    case 3:
+        P->e = pj_atof (pj_param_value (par));
+        if (HUGE_VAL==P->e)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        if (0==P->e)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        if (1==P->e)
+             return proj_errno_set (P, PJD_ERR_ECCENTRICITY_IS_ONE);
+        P->es = P->e * P->e;
+        break;
+
+    /* semiminor axis, b */
+    case 4:
+        P->b = pj_atof (pj_param_value (par));
+        if (HUGE_VAL==P->b)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        if (0==P->b)
+            return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+        P->f = (P->a - P->b) / P->a;
+        P->es = 2*P->f - P->f*P->f;
+        break;
+    default:
+        return PJD_ERR_INVALID_ARG;
+
+    }
+    if (P->es < 0)
+        return proj_errno_set (P, PJD_ERR_ES_LESS_THAN_ZERO);
+    return 0;
+}
+
+
+/***************************************************************************************/
+static int pj_spherification (PJ *P) {
+/***************************************************************************************/
+    char *keys[] =  {"R_A", "R_V", "R_a", "R_g", "R_h", "R_lat_a", "R_lat_g"};
+    size_t len, i;
+    paralist *par = 0;
+    char *def = 0;
+
+    double t;
+    char *v, *endp;
+
+    len = sizeof (keys) /  sizeof (char *);
+    P->def_spherification = 0;
+
+    /* Check which spherification key is specified */
+    for (i = 0;  i < len;  i++) {
+        par = pj_get_param (P->params, keys[i]);
+        if (par)
+            break;
+    }
+
+    /* No spherification specified? Then we're done */
+    if (i==len)
+        return 0;
+
+    /* Store definition */
+    P->def_spherification = def = par->param;
+    par->used = 1;
+
+    switch (i) {
+
+    /* R_A - a sphere with same area as ellipsoid */
+    case 0:
+        P->a *= 1. - P->es * (SIXTH + P->es * (RA4 + P->es * RA6));
+        break;
+
+    /* R_V - a sphere with same volume as ellipsoid */
+    case 1:
+        P->a *= 1. - P->es * (SIXTH + P->es * (RV4 + P->es * RV6));
+        break;
+
+    /* R_a - a sphere with R = the arithmetic mean of the ellipsoid */
+    case 2:
+        P->a = (P->a + P->b) / 2;
+        break;
+
+    /* R_g - a sphere with R = the geometric mean of the ellipsoid */
+    case 3:
+        P->a = sqrt (P->a * P->b);
+        break;
+
+    /* R_h - a sphere with R = the harmonic mean of the ellipsoid */
+    case 4:
+        if (P->a + P->b == 0)
+            return proj_errno_set (P, PJD_ERR_TOLERANCE_CONDITION);
+        P->a = (2*P->a * P->b) / (P->a + P->b);
+        break;
+
+    /* R_lat_a - a sphere with R = the arithmetic mean of the ellipsoid at given latitude */
+    case 5:
+    /* R_lat_g - a sphere with R = the geometric  mean of the ellipsoid at given latitude */
+    case 6:
+        v = pj_param_value (par);
+        t = proj_dmstor (v, &endp);
+        if (fabs (t) > M_HALFPI)
+            return proj_errno_set (P, PJD_ERR_REF_RAD_LARGER_THAN_90);
+        t = sin (t);
+        t = 1 - P->es * t * t;
+        if (i==5)   /* arithmetic */
+            P->a *= (1. - P->es + t) / (2 * t * sqrt(t));
+        else        /* geometric */
+            P->a *= sqrt (1 - P->es) / t;
+        break;
+    }
+
+    /* Clean up the ellipsoidal parameters to reflect the sphere */
+    P->es = P->e = P->rf = P->f = 0;
+    P->b = P->a;
+    pj_calc_ellps_params (P, P->a, 0);
+
+    return 0;
+}
+
+
+
+
+
+/* locate parameter in list */
+static paralist *pj_get_param (paralist *list, char *key) {
+    size_t l = strlen(key);
+	while (list && !(0==strncmp(list->param, key, l) && (0==list->param[l] || list->param[l] == '=') ) )
+		list = list->next;
+	return list;
+}
+
+
+static char *pj_param_value (paralist *list) {
+    char *key, *value;
+    if (0==list)
+        return 0;
+
+    key = list->param;
+    value = strchr (key, '=');
+
+    /* a flag (i.e. a key without value) has its own name (key) as value */
+    return value? value + 1: key;
+}
+
+
+static PJ_ELLPS *pj_find_ellps (char *name) {
+    int i;
+    char *s;
+    if (0==name)
+        return 0;
+
+    /* Search through internal ellipsoid list for name */
+    for (i = 0; (s = pj_ellps[i].id) && strcmp(name, s) ; ++i);
+    if (0==s)
+        return 0;
+    return pj_ellps + i;
+}
+
+
+
+
+
 /* copy ellipsoidal parameters from src to dst */
 void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
 
@@ -44,13 +455,41 @@ void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
     dst->a_orig  = src->a_orig;
 }
 
+
+/***************************************************************************************/
 int pj_calc_ellps_params(PJ *P, double a, double es) {
+/****************************************************************************************
+    Calculate a large number of ancillary ellipsoidal parameters, in addition to
+    the two traditional PROJ defining parameters: Semimajor axis, a, and the
+    eccentricity squared, es.
+
+    Most of these parameters are fairly cheap to compute in comparison to the overall
+    effort involved in initializing a PJ object. They may, however, take a substantial
+    part of the time taken in computing an individual point transformation.
+
+    So by providing them up front, we can amortize the (already modest) cost over all
+    transformations carried out over the entire lifetime of a PJ object, rather than
+    incur that cost for every single transformation.
+
+    Most of the parameter calculations here are based on the "angular eccentricity",
+    i.e. the angle, measured from the semiminor axis, of a line going from the north
+    pole to one of the foci of the ellipsoid - or in other words: The arc sine of the
+    eccentricity.
+
+    The formulae used are mostly taken from:
+
+    Richard H. Rapp: Geometric Geodesy, Part I, (178 pp, 1991).
+    Columbus, Ohio:  Dept. of Geodetic Science
+    and Surveying, Ohio State University.
+
+****************************************************************************************/
 
     P->a = a;
     P->es = es;
 
     /* Compute some ancillary ellipsoidal parameters */
-    P->e = sqrt(P->es);      /* eccentricity */
+    if (P->e==0)
+        P->e = sqrt(P->es);  /* eccentricity */
     P->alpha = asin (P->e);  /* angular eccentricity */
 
     /* second eccentricity */
@@ -62,7 +501,8 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
     P->e3s = P->e3 * P->e3;
 
     /* flattening */
-    P->f  = 1 - cos (P->alpha);   /* = 1 - sqrt (1 - PIN->es); */
+    if (0==P->f)
+        P->f  = 1 - cos (P->alpha);   /* = 1 - sqrt (1 - PIN->es); */
     P->rf = P->f != 0.0 ? 1.0/P->f: HUGE_VAL;
 
     /* second flattening */
@@ -74,7 +514,8 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
     P->rn = P->n != 0.0 ? 1/P->n: HUGE_VAL;
 
     /* ...and a few more */
-    P->b  = (1 - P->f)*P->a;
+    if (0==P->b)
+        P->b  = (1 - P->f)*P->a;
     P->rb = 1. / P->b;
     P->ra = 1. / P->a;
 
@@ -89,8 +530,45 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
     return 0;
 }
 
-/* initialize geographic shape parameters */
+
+
+#ifndef KEEP_ORIGINAL_PJ_ELL_SET
 int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
+    PJ B;
+    int ret;
+
+    memset (&B, 0, sizeof (B));
+    B.ctx = ctx;
+    B.params = pl;
+
+    ret = pj_ellipsoid (&B);
+    if (ret)
+        return ret;
+
+    *a = B.a;
+    *es = B.es;
+    return 0;
+}
+#else
+
+
+
+/**************************************************************************************/
+int pj_ell_set (projCtx ctx, paralist *pl, double *a, double *es) {
+/***************************************************************************************
+    Initialize ellipsoidal parameters: This is the original ellipsoid setup
+    function by Gerald Evenden - significantly more compact than pj_ellipsoid and
+    its many helper functions, and still quite readable.
+
+    It is, however, also so tight that it is hard to modify and add functionality,
+    and equally hard to find the right place to add further commentary for improved
+    future maintainability.
+
+    Hence, when the need to store in the PJ object, the parameters actually used to
+    define the ellipsoid came up, rather than modifying this little gem of
+    "economy of expression", a much more verbose reimplementation, pj_ellipsoid,
+    was written.
+***************************************************************************************/
     int i;
     double b=0.0, e;
     char *name;
@@ -120,6 +598,7 @@ int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
         }
 
         *a = pj_param(ctx,pl, "da").f;
+
         if (pj_param(ctx,pl, "tes").i) /* eccentricity squared */
             *es = pj_param(ctx,pl, "des").f;
         else if (pj_param(ctx,pl, "te").i) { /* eccentricity */
@@ -142,6 +621,8 @@ int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
         }     /* else *es == 0. and sphere of radius *a */
         if (b == 0.0)
             b = *a * sqrt(1. - *es);
+
+
         /* following options turn ellipsoid into equivalent sphere */
         if (pj_param(ctx,pl, "bR_A").i) { /* sphere--area of ellipsoid */
             *a *= 1. - *es * (SIXTH + *es * (RA4 + *es * RA6));
@@ -192,3 +673,4 @@ bomb:
             { pj_ctx_set_errno( ctx, -13); return 1; }
     return 0;
 }
+#endif
diff --git a/src/pj_init.c b/src/pj_init.c
index 53fd1039..0cdfcc4d 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -536,7 +536,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
 
     if (PIN->need_ellps) {
-        if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) {
+        if (0 != pj_ellipsoid (PIN)) {
             pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
             return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
         }
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index c003c717..c9275074 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -103,7 +103,7 @@ void pj_dalloc(void *ptr) {
 /**********************************************************************/
     free(ptr);
 }
-    
+
 
 /**********************************************************************/
 void *pj_dealloc (void *ptr) {
@@ -150,7 +150,7 @@ void *pj_dealloc_params (projCtx ctx, paralist *start, int errlev) {
 
     Also called from pj_init_ctx when encountering errors before the PJ
     proper is allocated.
-******************************************************************************/    
+******************************************************************************/
     paralist *t, *n;
     for (t = start; t; t = n) {
         n = t->next;
@@ -177,12 +177,12 @@ void *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
 
     if (0==P)
         return 0;
-    
+
     /* free grid lists */
     pj_dealloc( P->gridlist );
     pj_dealloc( P->vgridlist_geoid );
     pj_dealloc( P->catalog_name );
-    
+
     /* We used to call pj_dalloc( P->catalog ), but this will leak */
     /* memory. The safe way to clear catalog and grid is to call */
     /* pj_gc_unloadall(pj_get_default_ctx()); and pj_deallocate_grids(); */
@@ -191,7 +191,7 @@ void *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
 
     /* free the interface to Charles Karney's geodesic library */
     pj_dealloc( P->geod );
-    
+
     /* free parameter list elements */
     pj_dealloc_params (pj_get_ctx(P), P->params, errlev);
 
diff --git a/src/pj_param.c b/src/pj_param.c
index 93f5cf53..ee952eca 100644
--- a/src/pj_param.c
+++ b/src/pj_param.c
@@ -2,8 +2,9 @@
 #include <projects.h>
 #include <stdio.h>
 #include <string.h>
-	paralist * /* create parameter list entry */
-pj_mkparam(char *str) {
+
+/* create parameter list entry */
+paralist *pj_mkparam(char *str) {
 	paralist *newitem;
 
 	if((newitem = (paralist *)pj_malloc(sizeof(paralist) + strlen(str))) != NULL) {
@@ -16,6 +17,7 @@ pj_mkparam(char *str) {
 	return newitem;
 }
 
+
 /************************************************************************/
 /*                              pj_param()                              */
 /*                                                                      */
diff --git a/src/proj.def b/src/proj.def
index c282525c..ca10489a 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -147,3 +147,6 @@ EXPORTS
 
     proj_angular_input             @133
     proj_angular_output            @134
+
+    pj_ellipsoid                   @135
+    pj_calc_ellps_params           @136
diff --git a/src/proj.h b/src/proj.h
index fbd59b45..f4640e15 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -382,9 +382,9 @@ double proj_xyz_dist (XYZ a, XYZ b);
 
 /* Set or read error level */
 int  proj_errno (PJ *P);
-void proj_errno_set (PJ *P, int err);
+int  proj_errno_set (PJ *P, int err);
 int  proj_errno_reset (PJ *P);
-void proj_errno_restore (PJ *P, int err);
+int  proj_errno_restore (PJ *P, int err);
 
 
 PJ_DERIVS  proj_derivatives(PJ *P, const LP lp);
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 7196b90c..4f044101 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -520,37 +520,33 @@ int proj_errno (PJ *P) {
 }
 
 /*****************************************************************************/
-void proj_errno_set (PJ *P, int err) {
+int proj_errno_set (PJ *P, int err) {
 /******************************************************************************
-    Sets errno at the context and bubble it up to the thread local errno
+    Set context-errno, bubble it up to the thread local errno, return err
 ******************************************************************************/
     /* Use proj_errno_reset to explicitly clear the error status */
     if (0==err)
-        return;
+        return 0;
 
     /* For P==0 err goes to the default context */
     proj_context_errno_set (pj_get_ctx (P), err);
     errno = err;
-    return;
+    return err;
 }
 
 /*****************************************************************************/
-void proj_errno_restore (PJ *P, int err) {
+int proj_errno_restore (PJ *P, int err) {
 /******************************************************************************
-    Reduce some mental impedance in the canonical reset/restore use case:
-    Basically, proj_errno_restore() is a synonym for proj_errno_set(),
-    but the use cases are very different (_set: indicate an error to higher
-    level user code, _restore: pass previously set error indicators in case
-    of no errors at this level).
-
-    Hence, although the inner working is identical, we provide both options,
-    to avoid some rather confusing real world code.
+    Use proj_errno_restore when the current function succeeds, but the
+    error flag was set on entry, and stored/reset using proj_errno_reset
+    in order to monitor for new errors.
 
     See usage example under proj_errno_reset ()
 ******************************************************************************/
     if (0==err)
-        return;
+        return 0;
     proj_errno_set (P, err);
+    return 0;
 }
 
 /*****************************************************************************/
@@ -562,17 +558,24 @@ int proj_errno_reset (PJ *P) {
     Returns the previous value of the errno, for convenient reset/restore
     operations:
 
-    void foo (PJ *P) {
+    int foo (PJ *P) {
+        // errno may be set on entry, but we need to reset it to be able to
+        // check for errors from "do_something_with_P(P)"
         int last_errno = proj_errno_reset (P);
 
+        // local failure
+        if (0==P)
+            return proj_errno_set (P, 42);
+
+        // call to function that may fail
         do_something_with_P (P);
 
-        // failure - keep latest error status
+        // failure in do_something_with_P? - keep latest error status
         if (proj_errno(P))
-            return;
-        // success - restore previous error status
-        proj_errno_restore (P, last_errno);
-        return;
+            return proj_errno (P);
+
+        // success - restore previous error status, return 0
+        return proj_errno_restore (P, last_errno);
     }
 ******************************************************************************/
     int last_errno;
diff --git a/src/projects.h b/src/projects.h
index 32b74b08..f8f1fa38 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -229,6 +229,11 @@ struct PJconsts {
     projCtx_t *ctx;
     const char *descr;             /* From pj_list.h or individual PJ_*.c file */
     paralist *params;              /* Parameter list */
+    char *def_size;                /* Shape and size parameters extracted from params */
+    char *def_shape;
+    char *def_spherification;
+    char *def_ellps;
+
     struct geod_geodesic *geod;    /* For geodesic computations */
     struct pj_opaque *opaque;      /* Projection specific parameters, Defined in PJ_*.c */
     int inverted;                  /* Tell high level API functions to swap inv/fwd */
@@ -413,11 +418,6 @@ struct ARG_list {
 typedef union { double  f; int  i; char *s; } PROJVALUE;
 
 
-struct PJ_SELFTEST_LIST {
-    char    *id;                 /* projection keyword */
-    int     (* testfunc)(void);  /* projection entry point */
-};
-
 struct PJ_ELLPS {
     char    *id;           /* ellipse keyword name */
     char    *major;        /* a= value */
@@ -534,6 +534,7 @@ enum deprecated_constants_for_now_dropped_analytical_factors {
 #define PJD_ERR_LAT_0_IS_ZERO           -55
 #define PJD_ERR_ELLIPSOIDAL_UNSUPPORTED -56
 #define PJD_ERR_TOO_MANY_INITS          -57
+#define PJD_ERR_INVALID_ARG             -58
 
 struct projFileAPI_t;
 
@@ -676,6 +677,7 @@ double aacos(projCtx,double), aasin(projCtx,double), asqrt(double), aatan2(doubl
 PROJVALUE pj_param(projCtx ctx, paralist *, const char *);
 paralist *pj_mkparam(char *);
 
+int pj_ellipsoid (PJ *);
 int pj_ell_set(projCtx ctx, paralist *, double *, double *);
 int pj_datum_set(projCtx,paralist *, PJ *);
 int pj_prime_meridian_set(paralist *, PJ *);