Merge remote-tracking branch 'osgeo/master' into docs-release-4.10.0

45 files changed, 1640 insertions, 392 deletions

diff --git a/appveyor.yml b/appveyor.yml
index 9f66b22e..43669a40 100644
--- a/appveyor.yml
+++ b/appveyor.yml
@@ -62,6 +62,7 @@ test_script:
   - gie.exe ..\test\gie\more_builtins.gie
   - gie.exe ..\test\gie\deformation.gie
   - gie.exe ..\test\gie\axisswap.gie
+  - gie.exe ..\test\gie\ellipsoid.gie
 
 deploy: off
 
diff --git a/docs/plot/plotdefs.json b/docs/plot/plotdefs.json
index 7df9d390..b61466f0 100644
--- a/docs/plot/plotdefs.json
+++ b/docs/plot/plotdefs.json
@@ -143,6 +143,17 @@
         "type": "poly"
     },
     {
+        "filename": "ccon.png",
+        "latmax": 70,
+        "latmin": 34,
+        "lonmax": -30,
+        "lonmin": 68,
+        "name": "ccon",
+        "projstring": "+proj=ccon +lat_1=52 +lon_0=19",
+        "res": "med",
+        "type": "poly"
+    },
+    {
         "filename": "cea.png",
         "latmax": 90,
         "latmin": -90,
diff --git a/docs/source/references.rst b/docs/source/references.rst
index eaf5b17e..7d1d2216 100644
--- a/docs/source/references.rst
+++ b/docs/source/references.rst
@@ -29,3 +29,10 @@ References
 .. [Snyder1993] Snyder, 1993, Flattening the Earth, Chicago and London, The university of Chicago press
 
 .. [WeberMoore2013] Weber, E.D., and T.J. Moore. 2013. `Corrected Conversion Algorithms For The Calcofi Station Grid And Their Implementation In Several Computer Languages <http://calcofi.org/publications/calcofireports/v54/Vol_54_Weber.pdf>`__. California Cooperative Oceanic Fisheries Investigations Reports 54.
+
+
+.. [Zajac1978] A. Zajac, 1978, "Atlas of distribution of vascular plants in Poland (ATPOL)". Taxon 27(5/6), 481–484.
+
+.. [Komsta2016] L. Komsta, 2016, `ATPOL geobotanical grid revisited – a proposal of coordinate conversion algorithms <http://wydawnictwo.up.lublin.pl/annales/Agricultura/2016/1/03.pdf>`__. Annales UMCS Sectio E Agricultura 71(1), 31-37.
+
+.. [Verey2017] M. Verey, 2017, `Theoretical analysis and practical consequences of adopting an ATPOL grid model as a conical projection, defining the conversion of plane coordinates to the WGS - 84 ellipsoid <http://www.botany.pl/atpol/Siatka%20ATPOL%20w%20analitycznym%20ujeciu.pdf>`__. Fragmenta Floristica et Geobotanica Polonica (preprint submitted).
diff --git a/docs/source/usage/operations/projections/ccon.rst b/docs/source/usage/operations/projections/ccon.rst
new file mode 100644
index 00000000..6197f061
--- /dev/null
+++ b/docs/source/usage/operations/projections/ccon.rst
@@ -0,0 +1,118 @@
+.. _ccon:
+
+********************************************************************************
+Central Conic
+********************************************************************************
+
+.. image:: ./images/ccon.png
+   :scale: 50%
+   :alt:   Central Conic 
+
+This is central (centrographic) projection on cone tangent at ``lat_0`` latitude,
+identical with ``conic()`` projection from ``mapproj`` R package.
+
+Usage
+########
+
+This simple projection is rarely used, as it is not equidistant, equal-area, nor
+conformal. 
+
+An example of usage (and the main reason to implement this projection in proj4)
+is the ATPOL geobotanical grid of Poland, developed in Institute of Botany, 
+Jagiellonian University, Krakow, Poland in 1970s [Zajac1978]_. The grid was 
+originally handwritten on paper maps and further copied by hand. The projection
+(together with strange Earth radius) was chosen by its creators as the compromise
+fit to existing maps during first software development in DOS era. Many years later
+it is still de facto standard grid in Polish geobotanical research.
+
+The ATPOL coordinates can be achieved with with the following parameters:
+
+::
+
+   +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000 +x_0=330000 +y_0=-350000
+
+For more information see [Komsta2016]_ and [Verey2017]_.
+
+
+Forward projection
+==================
+
+.. math::
+
+   r = \cot \phi_0 - \tan (\phi - \phi_0)
+
+.. math::
+
+   x = r \sin (\lambda\sin\phi_0)
+
+.. math::
+
+   y = \cot \phi_0 - r \cos (\lambda\sin\phi_0) 
+
+
+Inverse projection
+==================
+
+.. math::
+
+   y = \cot \phi_0 - y
+
+.. math::
+
+   \phi = \phi_0 - \tan^{-1} ( \sqrt{x^2+y^2} - \cot \phi_0 )
+
+.. math::
+
+   \lambda = \frac{\tan^{-1} \sqrt{x^2+y^2}}{\sin \phi_0}
+
+Reference values
+==================
+
+For ATPOL to WGS84 test, run the following script:
+
+::
+
+   #!/bin/bash
+   cat << EOF | src/cs2cs -v -f "%E" +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000 +x_0=330000 +y_0=-350000 +to +proj=longlat +datum=WGS84 +no_defs
+   0 0
+   0 700000
+   700000 0
+   700000 700000
+   330000 350000
+   EOF
+
+It should result with
+
+::
+
+   1.384023E+01	5.503040E+01 0.000000E+00
+   1.451445E+01	4.877385E+01 0.000000E+00
+   2.478271E+01	5.500352E+01 0.000000E+00
+   2.402761E+01	4.875048E+01 0.000000E+00
+   1.900000E+01	5.200000E+01 0.000000E+00
+
+Analogous script can be run for reverse test:
+
+::
+
+   cat << EOF  | src/cs2cs -v -f "%E" +proj=longlat +datum=WGS84 +no_defs +to +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +axis=esu +a=6390000 +x_0=330000 +y_0=-350000
+   24 55
+   15 49
+   24 49
+   19 52
+   EOF
+
+and it should give the following results:
+
+::
+
+   6.500315E+05	4.106162E+03 0.000000E+00
+   3.707419E+04	6.768262E+05 0.000000E+00
+   6.960534E+05	6.722946E+05 0.000000E+00
+   3.300000E+05	3.500000E+05 0.000000E+00
+
+
+
+
+
+
diff --git a/docs/source/usage/operations/projections/images/ccon.png b/docs/source/usage/operations/projections/images/ccon.png
new file mode 100644
index 00000000..ccc0ec9e
--- /dev/null
+++ b/docs/source/usage/operations/projections/images/ccon.png
diff --git a/docs/source/usage/operations/projections/index.rst b/docs/source/usage/operations/projections/index.rst
index ae8ef1ee..5ccf0045 100644
--- a/docs/source/usage/operations/projections/index.rst
+++ b/docs/source/usage/operations/projections/index.rst
@@ -20,6 +20,7 @@ Projections
    calcofi
    cass
    cc
+   ccon
    cea
    chamb
    collg
diff --git a/src/Makefile.am b/src/Makefile.am
index 58514592..096cc672 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -46,11 +46,11 @@ libproj_la_SOURCES = \
 	pj_list.h proj_internal.h\
 	PJ_aeqd.c PJ_gnom.c PJ_laea.c PJ_mod_ster.c \
 	PJ_nsper.c PJ_nzmg.c PJ_ortho.c PJ_stere.c PJ_sterea.c \
-	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c \
+	PJ_aea.c PJ_bipc.c PJ_bonne.c PJ_eqdc.c PJ_isea.c PJ_ccon.c\
 	PJ_imw_p.c PJ_krovak.c PJ_lcc.c PJ_poly.c \
 	PJ_rpoly.c PJ_sconics.c proj_rouss.c \
-	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c \
-	PJ_gall.c PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
+	PJ_cass.c PJ_cc.c PJ_cea.c PJ_eqc.c PJ_gall.c PJ_geoc.c \
+	PJ_labrd.c PJ_lsat.c PJ_misrsom.c PJ_merc.c \
 	PJ_mill.c PJ_ocea.c PJ_omerc.c PJ_somerc.c \
 	PJ_tcc.c PJ_tcea.c PJ_times.c PJ_tmerc.c \
 	PJ_airy.c PJ_aitoff.c PJ_august.c PJ_bacon.c \
diff --git a/src/PJ_axisswap.c b/src/PJ_axisswap.c
index 53c93cc9..f8f17380 100644
--- a/src/PJ_axisswap.c
+++ b/src/PJ_axisswap.c
@@ -187,6 +187,10 @@ PJ *CONVERSION(axisswap,0) {
     /* read axes numbers and signs */
     for ( s = order, n = 0; *s != '\0' && n < 4; ) {
         Q->axis[n] = abs(atoi(s))-1;
+        if (Q->axis[n] >= 4) {
+            proj_log_error(P, "swapaxis: invalid axis '%s'", s);
+            return pj_default_destructor(P, PJD_ERR_AXIS);
+        }
         Q->sign[n++] = sign(atoi(s));
         while ( *s != '\0' && *s != ',' )
             s++;
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 2bcf3b4b..cd1995c1 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -187,7 +187,7 @@ static LPZ geodetic (XYZ cartesian,  PJ *P) {
 
 /* In effect, 2 cartesian coordinates of a point on the ellipsoid. Rather pointless, but... */
 static XY cart_forward (LP lp, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point;
     point.lp = lp;
     point.lpz.z = 0;
 
@@ -197,7 +197,7 @@ static XY cart_forward (LP lp, PJ *P) {
 
 /* And the other way round. Still rather pointless, but... */
 static LP cart_reverse (XY xy, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point;
     point.xy = xy;
     point.xyz.z = 0;
 
diff --git a/src/PJ_ccon.c b/src/PJ_ccon.c
new file mode 100644
index 00000000..cf8a9f68
--- /dev/null
+++ b/src/PJ_ccon.c
@@ -0,0 +1,102 @@
+/******************************************************************************
+ * Copyright (c) 2017, Lukasz Komsta
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#define PJ_LIB__
+#include <errno.h>
+#include <proj.h>
+#include "projects.h"
+
+struct pj_opaque {
+    double phi1;
+    double ctgphi1;
+    double sinphi1;
+    double cosphi1;
+    double *en;	
+};
+
+PROJ_HEAD(ccon, "Central Conic")
+    "\n\tCentral Conic, Sph.\n\tlat_1=";
+
+
+
+static XY forward (LP lp, PJ *P) {          
+    XY xy = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+    double r;
+
+    r = Q->ctgphi1 - tan(lp.phi - Q->phi1);
+    xy.x = r * sin(lp.lam * Q->sinphi1);
+    xy.y = Q->ctgphi1 - r * cos(lp.lam * Q->sinphi1);	
+
+    return xy;
+}
+
+
+static LP inverse (XY xy, PJ *P) {         
+    LP lp = {0.0,0.0};
+    struct pj_opaque *Q = P->opaque;
+
+    xy.y = Q->ctgphi1 - xy.y;
+    lp.phi = Q->phi1 - atan(hypot(xy.x,xy.y) - Q->ctgphi1);
+    lp.lam = atan2(xy.x,xy.y)/Q->sinphi1;
+
+    return lp;
+}
+
+
+static void *destructor (PJ *P, int errlev) {                       
+    if (0==P)
+        return 0;
+
+    if (0==P->opaque)
+        return pj_default_destructor (P, errlev);
+
+    pj_dealloc (P->opaque->en);
+    return pj_default_destructor (P, errlev);
+}
+
+
+PJ *PROJECTION(ccon) {
+
+    struct pj_opaque *Q = pj_calloc (1, sizeof (struct pj_opaque));
+    if (0==Q)
+        return pj_default_destructor (P, ENOMEM);
+    P->opaque = Q;
+    P->destructor = destructor;
+
+    Q->phi1 = pj_param(P->ctx, P->params, "rlat_1").f;
+
+    if (!(Q->en = pj_enfn(P->es)))
+        return pj_default_destructor(P, ENOMEM);
+
+    Q->sinphi1 = sin(Q->phi1);
+    Q->cosphi1 = cos(Q->phi1);
+    Q->ctgphi1 = Q->cosphi1/Q->sinphi1;
+
+
+    P->inv = inverse;
+    P->fwd = forward;
+    
+    return P;
+}
+
+
diff --git a/src/PJ_deformation.c b/src/PJ_deformation.c
index 81dd602d..09692ccb 100644
--- a/src/PJ_deformation.c
+++ b/src/PJ_deformation.c
@@ -232,9 +232,8 @@ PJ *TRANSFORMATION(deformation,1) {
     if (Q->cart == 0)
         return destructor(P, ENOMEM);
 
-    /* inherit ellipsoid definition from P to Q->cart (simpler than guessing */
-    /* how the ellipsoid was specified in original definition)               */
-    pj_inherit_ellipsoid_defs(P, Q->cart);
+    /* inherit ellipsoid definition from P to Q->cart */
+    pj_inherit_ellipsoid_def (P, Q->cart);
 
     Q->has_xy_grids = pj_param(P->ctx, P->params, "txy_grids").i;
     Q->has_z_grids  = pj_param(P->ctx, P->params, "tz_grids").i;
diff --git a/src/PJ_geoc.c b/src/PJ_geoc.c
new file mode 100644
index 00000000..865b1089
--- /dev/null
+++ b/src/PJ_geoc.c
@@ -0,0 +1,56 @@
+/******************************************************************************
+ * Project:  PROJ.4
+ * Purpose:  Conversion from geographic to geocentric latitude and back.
+ * Author:   Thomas Knudsen (2017)
+ *
+ ******************************************************************************
+ * Copyright (c) 2017, SDFE, http://www.sdfe.dk
+ * Copyright (c) 2017, Thomas Knudsen
+ *
+ * Permission is hereby granted, free of charge, to any person obtaining a
+ * copy of this software and associated documentation files (the "Software"),
+ * to deal in the Software without restriction, including without limitation
+ * the rights to use, copy, modify, merge, publish, distribute, sublicense,
+ * and/or sell copies of the Software, and to permit persons to whom the
+ * Software is furnished to do so, subject to the following conditions:
+ *
+ * The above copyright notice and this permission notice shall be included
+ * in all copies or substantial portions of the Software.
+ *
+ * THE SOFTWARE IS PROVIDED "AS IS", WITHOUT WARRANTY OF ANY KIND, EXPRESS
+ * OR IMPLIED, INCLUDING BUT NOT LIMITED TO THE WARRANTIES OF MERCHANTABILITY,
+ * FITNESS FOR A PARTICULAR PURPOSE AND NONINFRINGEMENT. IN NO EVENT SHALL
+ * THE AUTHORS OR COPYRIGHT HOLDERS BE LIABLE FOR ANY CLAIM, DAMAGES OR OTHER
+ * LIABILITY, WHETHER IN AN ACTION OF CONTRACT, TORT OR OTHERWISE, ARISING
+ * FROM, OUT OF OR IN CONNECTION WITH THE SOFTWARE OR THE USE OR OTHER
+ * DEALINGS IN THE SOFTWARE.
+ *****************************************************************************/
+
+#define PJ_LIB__
+#include <proj.h>
+#include <errno.h>
+#include "projects.h"
+
+PROJ_HEAD(geoc, "Geocentric Latitude");
+
+/* Geographical to geocentric */
+static PJ_COORD forward(PJ_COORD coo, PJ *P) {
+    return proj_geoc_lat (P, PJ_FWD, coo);
+}
+
+/* Geocentric to geographical */
+static PJ_COORD inverse(PJ_COORD coo, PJ *P) {
+    return proj_geoc_lat (P, PJ_INV, coo);
+}
+
+
+static PJ *CONVERSION(geoc, 1) {
+    P->inv4d = inverse;
+    P->fwd4d = forward;
+
+    P->left   =  PJ_IO_UNITS_RADIANS;
+    P->right  =  PJ_IO_UNITS_RADIANS;
+
+    P->is_latlong = 1;
+    return P;
+}
diff --git a/src/PJ_healpix.c b/src/PJ_healpix.c
index d8e1a921..41f62ee4 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_ellipsoid_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_helmert.c b/src/PJ_helmert.c
index 478a2c82..d0fd8222 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -310,7 +310,7 @@ static void build_rot_matrix(PJ *P) {
 static XY helmert_forward (LP lp, PJ *P) {
 /***********************************************************************/
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     double x, y, cr, sr;
     point.lp = lp;
 
@@ -330,7 +330,7 @@ static XY helmert_forward (LP lp, PJ *P) {
 static LP helmert_reverse (XY xy, PJ *P) {
 /***********************************************************************/
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     double x, y, sr, cr;
     point.xy = xy;
 
@@ -350,7 +350,7 @@ static LP helmert_reverse (XY xy, PJ *P) {
 static XYZ helmert_forward_3d (LPZ lpz, PJ *P) {
 /***********************************************************************/
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     double X, Y, Z, scale;
 
     point.lpz = lpz;
@@ -390,7 +390,7 @@ static XYZ helmert_forward_3d (LPZ lpz, PJ *P) {
 static LPZ helmert_reverse_3d (XYZ xyz, PJ *P) {
 /***********************************************************************/
     struct pj_opaque_helmert *Q = (struct pj_opaque_helmert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     double X, Y, Z, scale;
 
     point.xyz = xyz;
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 92275e51..5c2b944d 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -5,7 +5,7 @@
 PROJ_HEAD(hgridshift, "Horizontal grid shift");
 
 static XYZ forward_3d(LPZ lpz, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
     if (P->gridlist != NULL) {
@@ -19,7 +19,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
 
 
 static LPZ reverse_3d(XYZ xyz, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
     if (P->gridlist != NULL) {
diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c
index 1b0eb3a1..765e1d50 100644
--- a/src/PJ_molodensky.c
+++ b/src/PJ_molodensky.c
@@ -193,7 +193,7 @@ static LPZ calc_abridged_params(LPZ lpz, PJ *P) {
 
 
 static XY forward_2d(LP lp, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
 
     point.lp = lp;
     point.xyz = forward_3d(point.lpz, P);
@@ -203,7 +203,7 @@ static XY forward_2d(LP lp, PJ *P) {
 
 
 static LP reverse_2d(XY xy, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
 
     point.xy = xy;
     point.xyz.z = 0;
@@ -215,7 +215,7 @@ static LP reverse_2d(XY xy, PJ *P) {
 
 static XYZ forward_3d(LPZ lpz, PJ *P) {
     struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
 
     point.lpz = lpz;
 
@@ -243,7 +243,7 @@ static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
 
 static LPZ reverse_3d(XYZ xyz, PJ *P) {
     struct pj_opaque_molodensky *Q = (struct pj_opaque_molodensky *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     LPZ lpz;
 
     /* calculate parameters depending on the mode we are in */
diff --git a/src/PJ_pipeline.c b/src/PJ_pipeline.c
index e1568423..35f79213 100644
--- a/src/PJ_pipeline.c
+++ b/src/PJ_pipeline.c
@@ -305,6 +305,7 @@ static char **argv_params (paralist *params, size_t argc) {
 /* re-initialize P->geod.                                                    */
 static void set_ellipsoid(PJ *P) {
     paralist *cur, *attachment;
+    int err = proj_errno_reset (P);
 
     /* Break the linked list after the global args */
     attachment = 0;
@@ -317,16 +318,24 @@ 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;
+
+        /* reset an "unerror": In this special use case, the errno is    */
+        /* not an error signal, but just a reply from pj_ellipsoid,      */
+        /* telling us that "No - there was no ellipsoid definition in    */
+        /* the PJ you provided".                                         */
+        proj_errno_reset (P);
     }
 
-    pj_calc_ellps_params(P, P->a, P->es);
+    pj_calc_ellipsoid_params (P, P->a, P->es);
+
     geod_init(P->geod, P->a,  (1 - sqrt (1 - P->es)));
 
     /* Re-attach the dangling list */
     cur->next = attachment;
+    proj_errno_restore (P, err);
 }
 
 
@@ -345,7 +354,7 @@ PJ *OPERATION(pipeline,0) {
 
     P->opaque = pj_calloc (1, sizeof(struct pj_opaque));
     if (0==P->opaque)
-        return pj_default_destructor(P, ENOMEM);
+        return destructor(P, ENOMEM);
 
     argc = (int)argc_params (P->params);
     P->opaque->argv = argv = argv_params (P->params, argc);
@@ -396,6 +405,7 @@ PJ *OPERATION(pipeline,0) {
     for (i_current_step = i_first_step, i = 0;  i < nsteps;  i++) {
         int j;
         int  current_argc = 0;
+        int  err;
         PJ     *next_step = 0;
 
         /* Build a set of setup args for the current step */
@@ -415,14 +425,22 @@ PJ *OPERATION(pipeline,0) {
         for (j = 1;  j < current_argc; j++)
             proj_log_trace (P, "    %s", current_argv[j]);
 
-        next_step = pj_init_ctx (P->ctx, current_argc, current_argv);
+        err = proj_errno_reset (P);
 
+        next_step = proj_create_argv (P->ctx, current_argc, current_argv);
         proj_log_trace (P, "Pipeline: Step %d at %p", i, next_step);
+
         if (0==next_step) {
-            proj_log_error (P, "Pipeline: Bad step definition: %s", current_argv[0]);
-            return destructor (P, PJD_ERR_MALFORMED_PIPELINE); /* ERROR: bad pipeline def */
+            /* The step init failed, but possibly without setting errno. If so, we say "malformed" */
+            int err_to_report = proj_errno(P);
+            if (0==err_to_report)
+                err_to_report = PJD_ERR_MALFORMED_PIPELINE;
+            proj_log_error (P, "Pipeline: Bad step definition: %s (%s)", current_argv[0], pj_strerrno (err_to_report));
+            return destructor (P, err_to_report); /* ERROR: bad pipeline def */
         }
 
+        proj_errno_restore (P, err);
+
         /* Is this step inverted? */
         for (j = 0;  j < current_argc; j++)
             if (0==strcmp("inv", current_argv[j]))
@@ -430,7 +448,7 @@ PJ *OPERATION(pipeline,0) {
 
         P->opaque->pipeline[i+1] = next_step;
 
-        proj_log_trace (P, "Pipeline:    step done");
+        proj_log_trace (P, "Pipeline at [%p]:    step at [%p] done", P, next_step);
     }
 
     proj_log_trace (P, "Pipeline: %d steps built. Determining i/o characteristics", nsteps);
@@ -440,7 +458,5 @@ PJ *OPERATION(pipeline,0) {
 
     /* Now, correspondingly determine forward output (= reverse input) data type */
     P->right = pj_right (P->opaque->pipeline[nsteps]);
-
     return P;
 }
-
diff --git a/src/PJ_unitconvert.c b/src/PJ_unitconvert.c
index 523d8897..f951ebb6 100644
--- a/src/PJ_unitconvert.c
+++ b/src/PJ_unitconvert.c
@@ -204,7 +204,7 @@ static XY forward_2d(LP lp, PJ *P) {
     Forward unit conversions in the plane
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.lp = lp;
 
     point.xy.x *= pj_units[Q->xy_in_id].factor / pj_units[Q->xy_out_id].factor;
@@ -220,7 +220,7 @@ static LP reverse_2d(XY xy, PJ *P) {
     Reverse unit conversions in the plane
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.xy = xy;
 
     point.xy.x *= pj_units[Q->xy_out_id].factor / pj_units[Q->xy_in_id].factor;
@@ -236,7 +236,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
     Forward unit conversions the vertical component
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
     /* take care of the horizontal components in the 2D function */
@@ -253,7 +253,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
     Reverse unit conversions the vertical component
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
     /* take care of the horizontal components in the 2D function */
@@ -271,7 +271,7 @@ static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
     Forward conversion of time units
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_COORD out;
+    PJ_COORD out = {{0,0,0,0}};
 
     /* delegate unit conversion of physical dimensions to the 3D function */
     out.xyz = forward_3d(obs.lpz, P);
@@ -291,7 +291,7 @@ static PJ_COORD reverse_4d(PJ_COORD obs, PJ *P) {
     Reverse conversion of time units
 ************************************************************************/
     struct pj_opaque_unitconvert *Q = (struct pj_opaque_unitconvert *) P->opaque;
-    PJ_COORD out;
+    PJ_COORD out = {{0,0,0,0}};
 
     /* delegate unit conversion of physical dimensions to the 3D function */
     out.lpz = reverse_3d(obs.xyz, P);
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index 97faa240..850734fc 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -6,7 +6,7 @@ PROJ_HEAD(vgridshift, "Vertical grid shift");
 
 
 static XYZ forward_3d(LPZ lpz, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.lpz = lpz;
 
     if (P->vgridlist_geoid != NULL) {
@@ -20,7 +20,7 @@ static XYZ forward_3d(LPZ lpz, PJ *P) {
 
 
 static LPZ reverse_3d(XYZ xyz, PJ *P) {
-    PJ_TRIPLET point;
+    PJ_COORD point = {{0,0,0,0}};
     point.xyz = xyz;
 
     if (P->vgridlist_geoid != NULL) {
@@ -34,7 +34,7 @@ static LPZ reverse_3d(XYZ xyz, PJ *P) {
 
 
 static PJ_COORD forward_4d(PJ_COORD obs, PJ *P) {
-    PJ_COORD point;
+    PJ_COORD point = {{0,0,0,0}};
     point.xyz = forward_3d (obs.lpz, P);
     return point;
 }
diff --git a/src/cct.c b/src/cct.c
index 1bd25f1d..a1c275b5 100644
--- a/src/cct.c
+++ b/src/cct.c
@@ -171,7 +171,7 @@ int main(int argc, char **argv) {
     verbose   = opt_given (o, "v");
 
     if (opt_given (o, "o"))
-        fout = fopen (opt_arg (o, "output"), "rt");
+        fout = fopen (opt_arg (o, "output"), "wt");
     if (0==fout) {
         fprintf (stderr, "%s: Cannot open '%s' for output\n", o->progname, opt_arg (o, "output"));
         free (o);
@@ -206,7 +206,8 @@ int main(int argc, char **argv) {
     /* Setup transformation */
     P = proj_create_argv (0, o->pargc, o->pargv);
     if ((0==P) || (0==o->pargc)) {
-        fprintf (stderr, "%s: Bad transformation arguments. '%s -h' for help\n", o->progname, o->progname);
+        fprintf (stderr, "%s: Bad transformation arguments - (%s)\n    '%s -h' for help\n",
+                 o->progname, pj_strerrno (proj_errno(P)), o->progname);
         free (o);
         if (stdout != fout)
             fclose (fout);
@@ -232,6 +233,7 @@ int main(int argc, char **argv) {
 
     /* Loop over all records of all input files */
     while (opt_input_loop (o, optargs_file_format_text)) {
+        int err;
         void *ret = fgets (buf, 10000, o->input);
         opt_eof_handler (o);
         if (0==ret) {
@@ -259,20 +261,26 @@ int main(int argc, char **argv) {
             point.lpzt.lam = proj_torad (point.lpzt.lam);
             point.lpzt.phi = proj_torad (point.lpzt.phi);
         }
+        err = proj_errno_reset (P);
         point = proj_trans (P, direction, point);
-        if (proj_angular_output (P, direction)) {
-            point.lpzt.lam = proj_todeg (point.lpzt.lam);
-            point.lpzt.phi = proj_todeg (point.lpzt.phi);
-        }
 
         if (HUGE_VAL==point.xyzt.x) {
-            /* transformation error (TODO provide existing internal errmsg here) */
-            fprintf (fout, "# Record %d TRANSFORMATION ERROR: %s", (int) o->record_index, buf);
+            /* transformation error */
+            fprintf (fout, "# Record %d TRANSFORMATION ERROR: %s (%s)",
+                            (int) o->record_index, buf, pj_strerrno (proj_errno(P)));
+            proj_errno_restore (P, err);
             continue;
         }
+        proj_errno_restore (P, err);
 
         /* Time to print the result */
-        fprintf (fout, "%20.15f  %20.15f  %20.15f  %20.15f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
+        if (proj_angular_output (P, direction)) {
+            point.lpzt.lam = proj_todeg (point.lpzt.lam);
+            point.lpzt.phi = proj_todeg (point.lpzt.phi);
+            fprintf (fout, "%14.10f  %14.10f  %12.4f  %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
+        }
+        else
+            fprintf (fout, "%13.4f  %13.4f  %12.4f  %12.4f\n", point.xyzt.x, point.xyzt.y, point.xyzt.z, point.xyzt.t);
     }
 
     if (stdout != fout)
diff --git a/src/gie.c b/src/gie.c
index 4b89fe8f..9e377648 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -127,13 +127,16 @@ double proj_atof(const char *str);
 
 int   main(int argc, char **argv);
 
-static int   process_file (char *fname);
-static int   errmsg (int errlev, char *msg, ...);
+static int   process_file (const char *fname);
+static int   errmsg (int errlev, const char *msg, ...);
 static int   get_inp (FILE *f, char *inp, int size);
-static int   get_cmnd (char *inp, char *cmnd, int len);
-static char *get_args (char *inp);
-static int   dispatch (char *cmnd, char *args);
-static char *column (char *buf, int n);
+static int   get_cmnd (const char *inp, char *cmnd, int len);
+static const char *get_args (const char *inp);
+static int   dispatch (const char *cmnd, const char *args);
+static const char *column (const char *buf, int n);
+static int errno_from_err_const (const char *err_const);
+static const char *err_const_from_errno (int err);
+static void list_err_codes (void);
 
 
 
@@ -152,7 +155,7 @@ typedef struct {
     int grand_ok, grand_ko;
     size_t operation_lineno;
     double tolerance;
-    char *curr_file;
+    const char *curr_file;
     FILE *fout;
 } gie_ctx;
 
@@ -187,12 +190,14 @@ static const char usage[] = {
     "    -q                Quiet: Opposite of verbose. In quiet mode not even errors\n"
     "                      are reported. Only interaction is through the return code\n"
     "                      (0 on success, non-zero indicates number of FAILED tests)\n"
+    "    -l                List the PROJ internal system error codes\n"
     "--------------------------------------------------------------------------------\n"
     "Long Options:\n"
     "--------------------------------------------------------------------------------\n"
     "    --output          Alias for -o\n"
     "    --verbose         Alias for -v\n"
     "    --help            Alias for -h\n"
+    "    --list            Alias for -l\n"
     "--------------------------------------------------------------------------------\n"
     "Examples:\n"
     "--------------------------------------------------------------------------------\n"
@@ -205,10 +210,10 @@ static const char usage[] = {
 
 int main (int argc, char **argv) {
     int  i;
-    const char *longflags[]  = {"v=verbose", "q=quiet", "h=help", 0};
+    const char *longflags[]  = {"v=verbose", "q=quiet", "h=help", "l=list", 0};
     const char *longkeys[]   = {"o=output", 0};
 
-    o = opt_parse (argc, argv, "hvq", "o", longflags, longkeys);
+    o = opt_parse (argc, argv, "hlvq", "o", longflags, longkeys);
     if (0==o)
         return 0;
 
@@ -217,6 +222,8 @@ int main (int argc, char **argv) {
         return 0;
     }
 
+    if (opt_given (o, "l"))
+        list_err_codes ();
 
 
     T.verbosity = opt_given (o, "q");
@@ -273,11 +280,11 @@ static int another_success (void) {
 }
 
 
-static int process_file (char *fname) {
+static int process_file (const char *fname) {
     FILE *f;
     char inp[CMDLEN];
     char cmnd[1000];
-    char *args;
+    const char *args;
 
     lineno = level = 0;
     T.op_ok = T.total_ok = 0;
@@ -326,8 +333,11 @@ static int process_file (char *fname) {
 }
 
 
-/* return a pointer to the n'th column of buf */
-char *column (char *buf, int n) {
+/*****************************************************************************/
+const char *column (const char *buf, int n) {
+/*****************************************************************************
+    Return a pointer to the n'th column of buf. Coulmn numbers start at 0.
+******************************************************************************/
     int i;
     if (n <= 0)
         return buf;
@@ -343,11 +353,15 @@ char *column (char *buf, int n) {
 }
 
 
-/* interpret <args> as a numeric followed by a linear decadal prefix - return the properly scaled numeric */
-static double strtod_scaled (char *args, double default_scale) {
+/*****************************************************************************/
+static double strtod_scaled (const char *args, double default_scale) {
+/*****************************************************************************
+    Interpret <args> as a numeric followed by a linear decadal prefix.
+    Return the properly scaled numeric
+******************************************************************************/
     double s;
-    char *endp = args;
-    s = proj_strtod (args, &endp);
+    const char *endp = args;
+    s = proj_strtod (args, (char **) &endp);
     if (args==endp)
         return HUGE_VAL;
 
@@ -373,10 +387,7 @@ static double strtod_scaled (char *args, double default_scale) {
 }
 
 
-
-
-
-static int banner (char *args) {
+static int banner (const char *args) {
     char dots[] = {"..."}, nodots[] = {""}, *thedots = nodots;
     if (strlen(args) > 70)
         thedots = dots;
@@ -385,11 +396,7 @@ static int banner (char *args) {
 }
 
 
-
-
-
-
-static int tolerance (char *args) {
+static int tolerance (const char *args) {
     T.tolerance = strtod_scaled (args, 1);
     if (HUGE_VAL==T.tolerance) {
         T.tolerance = 0.0005;
@@ -399,8 +406,8 @@ static int tolerance (char *args) {
 }
 
 
-static int direction (char *args) {
-    char *endp = args;
+static int direction (const char *args) {
+    const char *endp = args;
     while (isspace (*endp))
         endp++;
     switch (*endp) {
@@ -421,14 +428,14 @@ static int direction (char *args) {
 }
 
 
-
-
-static void finish_previous_operation (char *args) {
+static void finish_previous_operation (const char *args) {
     if (T.verbosity > 1 && T.op_id > 1 && T.op_ok+T.op_ko)
         fprintf (T.fout, "%s     %d tests succeeded,  %d tests %s\n", delim, T.op_ok, T.op_ko, T.op_ko? "FAILED!": "failed.");
     (void) args;
 }
 
+
+
 /*****************************************************************************/
 static int operation (char *args) {
 /*****************************************************************************
@@ -482,8 +489,9 @@ static int operation (char *args) {
 static int pj_unitconvert_selftest (void);
 static int pj_cart_selftest (void);
 static int pj_horner_selftest (void);
+
 /*****************************************************************************/
-static int builtins (char *args) {
+static int builtins (const char *args) {
 /*****************************************************************************
     There are still a few tests that cannot be described using gie
     primitives. Instead, they are implemented as builtins, and invoked
@@ -526,9 +534,6 @@ static int builtins (char *args) {
 }
 
 
-
-
-
 static PJ_COORD torad_coord (PJ_COORD a) {
     PJ_COORD c = a;
     c.lpz.lam = proj_torad (a.lpz.lam);
@@ -536,6 +541,7 @@ static PJ_COORD torad_coord (PJ_COORD a) {
     return c;
 }
 
+
 static PJ_COORD todeg_coord (PJ_COORD a) {
     PJ_COORD c = a;
     c.lpz.lam = proj_todeg (a.lpz.lam);
@@ -543,14 +549,19 @@ static PJ_COORD todeg_coord (PJ_COORD a) {
     return c;
 }
 
-/* try to parse args as a PJ_COORD */
-static PJ_COORD parse_coord (char *args) {
+
+
+/*****************************************************************************/
+static PJ_COORD parse_coord (const char *args) {
+/*****************************************************************************
+    Attempt to interpret args as a PJ_COORD.
+******************************************************************************/
     int i;
-    char *endp, *prev = args;
+    const char *endp, *prev = args;
     PJ_COORD a = proj_coord (0,0,0,0);
 
     for (i = 0; i < 4; i++) {
-        double d = proj_strtod (prev, &endp);
+        double d = proj_strtod (prev,  (char **) &endp);
         if (prev==endp)
             return i > 1? a: proj_coord_error ();
         a.v[i] = d;
@@ -562,7 +573,7 @@ static PJ_COORD parse_coord (char *args) {
 
 
 /*****************************************************************************/
-static int accept (char *args) {
+static int accept (const char *args) {
 /*****************************************************************************
     Read ("ACCEPT") a 2, 3, or 4 dimensional input coordinate.
 ******************************************************************************/
@@ -573,9 +584,8 @@ static int accept (char *args) {
 }
 
 
-
 /*****************************************************************************/
-static int roundtrip (char *args) {
+static int roundtrip (const char *args) {
 /*****************************************************************************
     Check how far we go from the ACCEPTed point when doing successive
     back/forward transformation pairs.
@@ -612,8 +622,7 @@ static int roundtrip (char *args) {
 }
 
 
-
-static int expect_message (double d, char *args) {
+static int expect_message (double d, const char *args) {
     another_failure ();
 
     if (T.verbosity < 0)
@@ -636,7 +645,8 @@ static int expect_message (double d, char *args) {
     return 1;
 }
 
-static int expect_message_cannot_parse (char *args) {
+
+static int expect_message_cannot_parse (const char *args) {
     another_failure ();
     if (T.verbosity > -1) {
         if (0==T.op_ko && T.verbosity < 2)
@@ -647,37 +657,78 @@ static int expect_message_cannot_parse (char *args) {
     return 1;
 }
 
+static int expect_failure_with_errno_message (int expected, int got) {
+    another_failure ();
+
+    if (T.verbosity < 0)
+        return 1;
+    if (0==T.op_ko && T.verbosity < 2)
+        banner (T.operation);
+    fprintf (T.fout, "%s", T.op_ko? "     -----\n": delim);
+    fprintf (T.fout, "     FAILURE in %s(%d):\n",    opt_strip_path (T.curr_file), (int) lineno);
+    fprintf (T.fout, "     got errno %s (%d): %s\n", err_const_from_errno(got), got,  pj_strerrno (got));
+    fprintf (T.fout, "     expected %s (%d):  %s",   err_const_from_errno(expected), expected, pj_strerrno (expected));
+    fprintf (T.fout, "\n");
+    return 1;
+}
+
 
 /*****************************************************************************/
-static int expect (char *args) {
+static int expect (const char *args) {
 /*****************************************************************************
     Tell GIE what to expect, when transforming the ACCEPTed input
 ******************************************************************************/
     PJ_COORD ci, co, ce;
     double d;
     int expect_failure = 0;
+    int expect_failure_with_errno = 0;
 
-    if (0==strcmp (args, "failure"))
+    if (0==strncmp (args, "failure", 7)) {
         expect_failure = 1;
 
-    if (0==T.P && !expect_failure) {
+        /* Option: Fail with an expected errno (syntax: expect failure errno -33) */
+        if (0==strncmp (column (args, 2), "errno", 5))
+            expect_failure_with_errno = errno_from_err_const (column (args, 3));
+    }
+
+    if (0==T.P) {
+        /* If we expect failure, and fail, then it's a success... */
+        if (expect_failure) {
+            /* Failed to fail correctly? */
+            if (expect_failure_with_errno && proj_errno (T.P)!=expect_failure_with_errno)
+                return expect_failure_with_errno_message (expect_failure_with_errno, proj_errno(T.P));
+            return another_success ();
+        }
+
+        /* Otherwise, it's a true failure */
         banner (T.operation);
-        errmsg(3, "%sInvalid operation definition in line no. %d\n", delim, (int) T.operation_lineno);
+        errmsg(3, "%sInvalid operation definition in line no. %d: %s\n",
+            delim, (int) T.operation_lineno, pj_strerrno(proj_errno(T.P)));
         return another_failure ();
     }
 
-
+    /* We may still successfully fail even if the proj_create succeeded */
     if (expect_failure) {
-        /* If we expect failure, and fail, then it's a success... */
-        if (0==T.P)
-            return another_success ();
-        /* We may still successfully fail even if the proj_create succeeded */
+        proj_errno_reset (T.P);
+
+        /* Try to carry out the operation - and expect failure */
         ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
         co = proj_trans (T.P, T.dir, ci);
-        if (co.xyz.x==HUGE_VAL)
-            return another_success ();
-        /* no - we didn't manage to purportedly fail */
-        return another_failure ();
+
+        /* Failed to fail? - that's a failure */
+        if (co.xyz.x!=HUGE_VAL)
+            return another_failure ();
+
+        if (expect_failure_with_errno) {
+            printf ("errno=%d, expected=%d\n", proj_errno (T.P), expect_failure_with_errno);
+            if (proj_errno (T.P)==expect_failure_with_errno)
+                return another_success ();
+
+            return another_failure ();
+        }
+
+        /* Yes, we failed successfully */
+        return another_success ();
     }
 
 
@@ -724,7 +775,6 @@ static int expect (char *args) {
             e = proj_xyz_dist (T.b.xyz, T.e.xyz);
         if (e < d)
             d = e;
-
     }
     else
         d = proj_xyz_dist (T.b.xyz, T.e.xyz);
@@ -739,7 +789,7 @@ static int expect (char *args) {
 
 
 /*****************************************************************************/
-static int verbose (char *args) {
+static int verbose (const char *args) {
 /*****************************************************************************
     Tell the system how noisy it should be
 ******************************************************************************/
@@ -757,7 +807,7 @@ static int verbose (char *args) {
 }
 
 /*****************************************************************************/
-static int comment (char *args) {
+static int comment (const char *args) {
 /*****************************************************************************
     in line comment. Equivalent to #
 ******************************************************************************/
@@ -767,7 +817,7 @@ static int comment (char *args) {
 
 
 /*****************************************************************************/
-static int echo (char *args) {
+static int echo (const char *args) {
 /*****************************************************************************
     Add user defined noise to the output stream
 ******************************************************************************/
@@ -777,17 +827,19 @@ fprintf (T.fout, "%s\n", args);
 
 
 
-static int dispatch (char *cmnd, char *args) {
+static int dispatch (const char *cmnd, const char *args) {
+#if 0
     int last_errno = proj_errno_reset (T.P);
+#endif
 
     if  (0==level%2) {
-        if (0==strcmp (cmnd, "BEGIN"))
+        if (0==strcmp (cmnd, "BEGIN") || 0==strcmp (cmnd, "<begin>") || 0==strcmp (cmnd, "<gie>"))
            level++;
         return 0;
     }
 
-    if  (0==strcmp (cmnd, "OPERATION")) return  operation (args);
-    if  (0==strcmp (cmnd, "operation")) return  operation (args);
+    if  (0==strcmp (cmnd, "OPERATION")) return  operation ((char *) args);
+    if  (0==strcmp (cmnd, "operation")) return  operation ((char *) args);
     if  (0==strcmp (cmnd, "ACCEPT"))    return  accept    (args);
     if  (0==strcmp (cmnd, "accept"))    return  accept    (args);
     if  (0==strcmp (cmnd, "EXPECT"))    return  expect    (args);
@@ -807,22 +859,160 @@ static int dispatch (char *cmnd, char *args) {
     if  (0==strcmp (cmnd, "ECHO"))      return  echo      (args);
     if  (0==strcmp (cmnd, "echo"))      return  echo      (args);
     if  (0==strcmp  (cmnd, "END"))      return          finish_previous_operation (args), level++, 0;
+    if  (0==strcmp  (cmnd, "<end>"))    return          finish_previous_operation (args), level++, 0;
+    if  (0==strcmp  (cmnd, "</gie>"))   return          finish_previous_operation (args), level++, 0;
     if  ('#'==cmnd[0])                  return  comment   (args);
 
+#if 0
     if (proj_errno(T.P))
         printf ("#####***** ERRNO=%d\n", proj_errno(T.P));
     proj_errno_restore (T.P, last_errno);
+#endif
     return 0;
 }
 
 
 
 
+struct errno_vs_err_const {const char *the_err_const; int the_errno;};
+static const struct errno_vs_err_const lookup[] = {
+    {"pjd_err_no_args"                  ,  -1},
+    {"pjd_err_no_option_in_init_file"   ,  -2},
+    {"pjd_err_no_colon_in_init_string"  ,  -3},
+    {"pjd_err_proj_not_named"           ,  -4},
+    {"pjd_err_unknown_projection_id"    ,  -5},
+    {"pjd_err_eccentricity_is_one"      ,  -6},
+    {"pjd_err_unknow_unit_id"           ,  -7},
+    {"pjd_err_invalid_boolean_param"    ,  -8},
+    {"pjd_err_unknown_ellp_param"       ,  -9},
+    {"pjd_err_rev_flattening_is_zero"   ,  -10},
+    {"pjd_err_ref_rad_larger_than_90"   ,  -11},
+    {"pjd_err_es_less_than_zero"        ,  -12},
+    {"pjd_err_major_axis_not_given"     ,  -13},
+    {"pjd_err_lat_or_lon_exceed_limit"  ,  -14},
+    {"pjd_err_invalid_x_or_y"           ,  -15},
+    {"pjd_err_wrong_format_dms_value"   ,  -16},
+    {"pjd_err_non_conv_inv_meri_dist"   ,  -17},
+    {"pjd_err_non_con_inv_phi2"         ,  -18},
+    {"pjd_err_acos_asin_arg_too_large"  ,  -19},
+    {"pjd_err_tolerance_condition"      ,  -20},
+    {"pjd_err_conic_lat_equal"          ,  -21},
+    {"pjd_err_lat_larger_than_90"       ,  -22},
+    {"pjd_err_lat1_is_zero"             ,  -23},
+    {"pjd_err_lat_ts_larger_than_90"    ,  -24},
+    {"pjd_err_control_point_no_dist"    ,  -25},
+    {"pjd_err_no_rotation_proj"         ,  -26},
+    {"pjd_err_w_or_m_zero_or_less"      ,  -27},
+    {"pjd_err_lsat_not_in_range"        ,  -28},
+    {"pjd_err_path_not_in_range"        ,  -29},
+    {"pjd_err_h_less_than_zero"         ,  -30},
+    {"pjd_err_k_less_than_zero"         ,  -31},
+    {"pjd_err_lat_1_or_2_zero_or_90"    ,  -32},
+    {"pjd_err_lat_0_or_alpha_eq_90"     ,  -33},
+    {"pjd_err_ellipsoid_use_required"   ,  -34},
+    {"pjd_err_invalid_utm_zone"         ,  -35},
+    {"pjd_err_tcheby_val_out_of_range"  ,  -36},
+    {"pjd_err_failed_to_find_proj"      ,  -37},
+    {"pjd_err_failed_to_load_grid"      ,  -38},
+    {"pjd_err_invalid_m_or_n"           ,  -39},
+    {"pjd_err_n_out_of_range"           ,  -40},
+    {"pjd_err_lat_1_2_unspecified"      ,  -41},
+    {"pjd_err_abs_lat1_eq_abs_lat2"     ,  -42},
+    {"pjd_err_lat_0_half_pi_from_mean"  ,  -43},
+    {"pjd_err_unparseable_cs_def"       ,  -44},
+    {"pjd_err_geocentric"               ,  -45},
+    {"pjd_err_unknown_prime_meridian"   ,  -46},
+    {"pjd_err_axis"                     ,  -47},
+    {"pjd_err_grid_area"                ,  -48},
+    {"pjd_err_invalid_sweep_axis"       ,  -49},
+    {"pjd_err_malformed_pipeline"       ,  -50},
+    {"pjd_err_unit_factor_less_than_0"  ,  -51},
+    {"pjd_err_invalid_scale"            ,  -52},
+    {"pjd_err_non_convergent"           ,  -53},
+    {"pjd_err_missing_args"             ,  -54},
+    {"pjd_err_lat_0_is_zero"            ,  -55},
+    {"pjd_err_ellipsoidal_unsupported"  ,  -56},
+    {"pjd_err_too_many_inits"           ,  -57},
+    {"pjd_err_invalid_arg"              ,  -58},
+    {"pjd_err_unknown"                  ,  9999},
+    {"pjd_err_enomem"                   ,  ENOMEM},
+};
+
+static const struct errno_vs_err_const unknown = {"PJD_ERR_UNKNOWN", 9999};
+
+
+static void list_err_codes (void) {
+    int i;
+    const int n = sizeof lookup / sizeof lookup[0];
+
+    for (i = 0;  i < n;  i++) {
+        if (9999==lookup[i].the_errno)
+            break;
+        printf ("%25s  (%2.2d):  %s\n", lookup[i].the_err_const + 8, lookup[i].the_errno, pj_strerrno(lookup[i].the_errno));
+    }
+    exit (0);
+}
 
 
+static const char *err_const_from_errno (int err) {
+    size_t i;
+    const size_t n = sizeof lookup / sizeof lookup[0];
+
+    for (i = 0;  i < n;  i++) {
+        if (err==lookup[i].the_errno)
+            return lookup[i].the_err_const + 8;
+    }
+    return unknown.the_err_const;
+}
 
 
-static int errmsg (int errlev, char *msg, ...) {
+
+static int errno_from_err_const (const char *err_const) {
+    const size_t n = sizeof lookup / sizeof lookup[0];
+    size_t i, len;
+    int ret;
+    char tolower_err_const[100];
+
+    /* Make a lower case copy for matching */
+    for (i = 0;  i < 99; i++) {
+        if (0==err_const[i] || isspace (err_const[i]))
+             break;
+        tolower_err_const[i] = (char) tolower (err_const[i]);
+    }
+    tolower_err_const[i] = 0;
+
+    /* If it looks numeric, return that numeric */
+    ret = (int) pj_atof (err_const);
+    if (0!=ret)
+        return ret;
+
+    /* Else try to find a matching identifier */
+    len = strlen (tolower_err_const);
+
+    /* First try to find a match excluding the PJD_ERR_ prefix */
+    for (i = 0;  i < n;  i++) {
+        if (0==strncmp (lookup[i].the_err_const + 8, err_const, len))
+            return lookup[i].the_errno;
+    }
+
+    /* If that did not work, try with the full name */
+    for (i = 0;  i < n;  i++) {
+        if (0==strncmp (lookup[i].the_err_const, err_const, len))
+            return lookup[i].the_errno;
+    }
+
+    /* On failure, return something unlikely */
+    return 9999;
+}
+
+
+
+
+
+
+
+
+static int errmsg (int errlev, const char *msg, ...) {
     va_list args;
     va_start(args, msg);
     vfprintf(stdout, msg, args);
@@ -894,7 +1084,7 @@ static int get_inp (FILE *f, char *inp, int size) {
     return (int) strlen(inp);
 }
 
-static int get_cmnd (char *inp, char *cmnd, int len) {
+static int get_cmnd (const char *inp, char *cmnd, int len) {
     cmnd[0] = 0;
     while (isspace(*inp++));
     inp--;
@@ -904,8 +1094,8 @@ static int get_cmnd (char *inp, char *cmnd, int len) {
     return len;
 }
 
-static char *get_args (char *inp) {
-    char *args = inp;
+static const char *get_args (const char *inp) {
+    const char *args = inp;
     while (isspace(*args++))
         if (0==*args)
             return args;
@@ -1038,7 +1228,6 @@ static int pj_cart_selftest (void) {
     PJ_GRID_INFO grid_info;
     PJ_INIT_INFO init_info;
 
-    PJ_DERIVS derivs;
     PJ_FACTORS factors;
 
     const PJ_OPERATIONS *oper_list;
@@ -1333,25 +1522,14 @@ static int pj_cart_selftest (void) {
     a.lp.lam = PJ_TORAD(12);
     a.lp.phi = PJ_TORAD(55);
 
-    derivs = proj_derivatives(P, a.lp);
-    if (proj_errno(P))
-        return 80; /* derivs not created correctly */
-
-    if ( fabs(derivs.x_l - 1.0)     > 1e-5 )   return 81;
-    if ( fabs(derivs.x_p - 0.0)     > 1e-5 )   return 82;
-    if ( fabs(derivs.y_l - 0.0)     > 1e-5 )   return 83;
-    if ( fabs(derivs.y_p - 1.73959) > 1e-5 )   return 84;
-
-
     factors = proj_factors(P, a.lp);
     if (proj_errno(P))
         return 85; /* factors not created correctly */
 
     /* check a few key characteristics of the Mercator projection */
-    if (factors.omega != 0.0)       return 86; /* angular distortion should be 0 */
-    if (factors.thetap != M_PI_2)   return 87; /* Meridian/parallel angle should be 90 deg */
-    if (factors.conv != 0.0)        return 88; /* meridian convergence should be 0 */
-
+    if (factors.angular_distortion != 0.0)  return 86; /* angular distortion should be 0 */
+    if (factors.meridian_parallel_angle != M_PI_2)  return 87; /* Meridian/parallel angle should be 90 deg */
+    if (factors.meridian_convergence != 0.0)  return 88; /* meridian convergence should be 0 */
 
     proj_destroy(P);
 
@@ -1446,7 +1624,7 @@ static int pj_cart_selftest (void) {
 
 
 
-static int test_time(char* args, double tol, double t_in, double t_exp) {
+static int test_time(const char* args, double tol, double t_in, double t_exp) {
     PJ_COORD in, out;
     PJ *P = proj_create(PJ_DEFAULT_CTX, args);
     int ret = 0;
@@ -1472,8 +1650,8 @@ static int test_time(char* args, double tol, double t_in, double t_exp) {
     return ret;
 }
 
-static int test_xyz(char* args, double tol, PJ_TRIPLET in, PJ_TRIPLET exp) {
-    PJ_COORD out, obs_in;
+static int test_xyz(const char* args, double tol, PJ_COORD in, PJ_COORD exp) {
+    PJ_COORD out = {{0,0,0,0}}, obs_in = {{0,0,0,0}};
     PJ *P = proj_create(PJ_DEFAULT_CTX, args);
     int ret = 0;
 
@@ -1515,10 +1693,10 @@ static int pj_unitconvert_selftest (void) {
     double in4 = 1877.71428, exp4 = 2016.0;
 
     char args5[] = "+proj=unitconvert +xy_in=m +xy_out=dm +z_in=cm +z_out=mm";
-    PJ_TRIPLET in5 = {{55.25, 23.23, 45.5}}, exp5 = {{552.5, 232.3, 455.0}};
+    PJ_COORD in5 = {{55.25, 23.23, 45.5, 0}}, exp5 = {{552.5, 232.3, 455.0, 0}};
 
     char args6[] = "+proj=unitconvert +xy_in=m +xy_out=m +z_in=m +z_out=m";
-    PJ_TRIPLET in6 = {{12.3, 45.6, 7.89}};
+    PJ_COORD in6 = {{12.3, 45.6, 7.89, 0}};
 
     ret = test_time(args1, 1e-6, in1, in1);   if (ret) return ret + 10;
     ret = test_time(args2, 1e-6, in2, in2);   if (ret) return ret + 20;
@@ -1531,5 +1709,3 @@ static int pj_unitconvert_selftest (void) {
 
 }
 
-
-
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 81de451e..88d88a97 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -54,6 +54,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_cart.c
         PJ_cass.c
         PJ_cc.c
+        PJ_ccon.c
         PJ_cea.c
         PJ_chamb.c
         PJ_collg.c
@@ -71,6 +72,7 @@ SET(SRC_LIBPROJ_PJ
         PJ_fahey.c
         PJ_fouc_s.c
         PJ_gall.c
+        PJ_geoc.c
         PJ_geos.c
         PJ_gins8.c
         PJ_gnom.c
diff --git a/src/makefile.vc b/src/makefile.vc
index 6bc06bd1..ef460719 100644
--- a/src/makefile.vc
+++ b/src/makefile.vc
@@ -11,7 +11,7 @@ azimuthal = \
 conic = \
 	PJ_aea.obj PJ_bipc.obj PJ_bonne.obj PJ_eqdc.obj \
 	PJ_imw_p.obj PJ_lcc.obj PJ_poly.obj \
-	PJ_rpoly.obj PJ_sconics.obj PJ_lcca.obj
+	PJ_rpoly.obj PJ_sconics.obj PJ_lcca.obj PJ_ccon.obj
 
 cylinder = \
 	PJ_cass.obj PJ_cc.obj PJ_cea.obj PJ_eqc.obj \
@@ -26,7 +26,7 @@ misc = \
 	PJ_lask.obj PJ_nocol.obj PJ_ob_tran.obj PJ_oea.obj \
 	PJ_sch.obj PJ_tpeqd.obj PJ_vandg.obj PJ_vandg2.obj \
     PJ_vandg4.obj PJ_wag7.obj PJ_latlong.obj PJ_krovak.obj \
-    pj_geocent.obj PJ_healpix.obj PJ_qsc.obj
+    PJ_geoc.obj pj_geocent.obj PJ_healpix.obj PJ_qsc.obj
 
 pseudo = \
 	PJ_boggs.obj PJ_collg.obj PJ_crast.obj PJ_denoy.obj \
diff --git a/src/optargpm.h b/src/optargpm.h
index 23e375f8..acb96583 100644
--- a/src/optargpm.h
+++ b/src/optargpm.h
@@ -207,7 +207,7 @@ static int opt_raise_flag (OPTARGS *opt, int ordinal);
 static int opt_ordinal (OPTARGS *opt, char *option);
 int opt_given (OPTARGS *opt, char *option);
 char *opt_arg (OPTARGS *opt, char *option);
-char *opt_strip_path (char *full_name);
+const char *opt_strip_path (const char *full_name);
 OPTARGS *opt_parse (int argc, char **argv, const char *flags, const char *keys, const char **longflags, const char **longkeys);
 
 #define opt_eof_handler(opt) if (opt_eof (opt)) {continue;} else {;}
@@ -219,7 +219,7 @@ struct OPTARGS {
     FILE *input;
     int   input_index;
     int   record_index;
-    char  *progname;         /* argv[0], stripped from /path/to, if present */
+    const char  *progname;   /* argv[0], stripped from /path/to, if present */
     char   flaglevel[21];    /* if flag -f is specified n times, its optarg pointer is set to flaglevel + n */
     char  *optarg[256];      /* optarg[(int) 'f'] holds a pointer to the argument of option "-f" */
     char  *flags;            /* a list of flag style options supported, e.g. "hv" (help and verbose) */
@@ -285,6 +285,8 @@ int opt_input_loop (OPTARGS *opt, int binary) {
 
     /* otherwise, open next input file */
     opt->input = fopen (opt->fargv[opt->input_index++], binary? "rb": "rt");
+    if (0 != opt->input)
+        return 1;
 
     /* ignore non-existing files - go on! */
     if (0==opt->input)
@@ -399,8 +401,8 @@ char *opt_arg (OPTARGS *opt, char *option) {
     return opt->optarg[ordinal];
 }
 
-char *opt_strip_path (char *full_name) {
-    char *last_path_delim, *stripped_name = full_name;
+const char *opt_strip_path (const char *full_name) {
+    const char *last_path_delim, *stripped_name = full_name;
 
     last_path_delim = strrchr (stripped_name, '\\');
     if (last_path_delim > stripped_name)
diff --git a/src/pj_ctx.c b/src/pj_ctx.c
index a8edaf43..326249ac 100644
--- a/src/pj_ctx.c
+++ b/src/pj_ctx.c
@@ -159,7 +159,7 @@ void pj_ctx_set_debug( projCtx ctx, int new_debug )
 
 {
     if (0==ctx)
-        pj_get_default_ctx ()->debug_level = new_debug;
+        return;
     ctx->debug_level = new_debug;
 }
 
diff --git a/src/pj_deriv.c b/src/pj_deriv.c
index 1c969a9d..7c1fcded 100644
--- a/src/pj_deriv.c
+++ b/src/pj_deriv.c
@@ -2,8 +2,12 @@
 #define PJ_LIB__
 #include "projects.h"
 
-int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
+int pj_deriv(LP lp, double h, const PJ *P, struct DERIVS *der) {
     XY t;
+    /* get rid of constness until we can do it for real */
+    PJ *Q = (PJ *) P;
+    if (0==Q->fwd)
+        return 1;
 
     lp.lam += h;
     lp.phi += h;
@@ -11,7 +15,7 @@ int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
         return 1;
 
     h += h;
-    t = (*P->fwd)(lp, P);
+    t = (*Q->fwd)(lp, Q);
     if (t.x == HUGE_VAL)
         return 1;
 
@@ -19,11 +23,12 @@ int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
     der->y_p = t.y;
     der->x_p = t.x;
     der->y_l = t.y;
+
     lp.phi -= h;
     if (fabs(lp.phi) > M_HALFPI)
         return 1;
 
-    t = (*P->fwd)(lp, P);
+    t = (*Q->fwd)(lp, Q);
     if (t.x == HUGE_VAL)
         return 1;
 
@@ -31,8 +36,9 @@ int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
     der->y_p -= t.y;
     der->x_p -= t.x;
     der->y_l += t.y;
+
     lp.lam -= h;
-    t = (*P->fwd)(lp, P);
+    t = (*Q->fwd)(lp, Q);
     if (t.x == HUGE_VAL)
         return 1;
 
@@ -40,8 +46,9 @@ int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
     der->y_p -= t.y;
     der->x_p -= t.x;
     der->y_l -= t.y;
+
     lp.phi += h;
-    t = (*P->fwd)(lp, P);
+    t = (*Q->fwd)(lp, Q);
     if (t.x == HUGE_VAL)
         return 1;
 
@@ -49,7 +56,9 @@ int pj_deriv(LP lp, double h, PJ *P, struct DERIVS *der) {
     der->y_p += t.y;
     der->x_p += t.x;
     der->y_l -= t.y;
-    der->x_l /= (h += h);
+
+    h += h;
+    der->x_l /= h;
     der->y_p /= h;
     der->x_p /= h;
     der->y_l /= h;
diff --git a/src/pj_ell_set.c b/src/pj_ell_set.c
index 8a9b39a1..d21122e3 100644
--- a/src/pj_ell_set.c
+++ b/src/pj_ell_set.c
@@ -1,15 +1,460 @@
 /* set ellipsoid parameters a and es */
 #include <string.h>
+#include <errno.h>
+#include <proj.h>
 #include "proj_internal.h"
 #include "projects.h"
-#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 */
 
-/* copy ellipsoidal parameters from src to dst */
-void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst) {
+
+/* Prototypes of the pj_ellipsoid helper functions */
+static int ellps_ellps (PJ *P);
+static int ellps_size (PJ *P);
+static int ellps_shape (PJ *P);
+static int ellps_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 are 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 err = proj_errno_reset (P);
+    char *empty = {""};
+
+    P->def_size = P->def_shape = P->def_spherification = P->def_ellps = 0;
+
+    /* Specifying R overrules everything */
+    if (pj_get_param (P->params, "R")) {
+        ellps_size (P);
+        pj_calc_ellipsoid_params (P, P->a, 0);
+        if (proj_errno (P))
+            return 1;
+        return proj_errno_restore (P, err);
+    }
+
+
+    /* If an ellps argument is specified, start by using that */
+    if (0 != ellps_ellps (P))
+        return 1;
+
+    /* We may overwrite the size */
+    if (0 != ellps_size (P))
+        return 2;
+
+    /* We may also overwrite the shape */
+    if (0 != ellps_shape (P))
+        return 3;
+
+    /* When we're done with it, we compute all related ellipsoid parameters */
+    pj_calc_ellipsoid_params (P, P->a, P->es);
+
+    /* And finally, we may turn it into a sphere */
+    if (0 != ellps_spherification (P))
+        return 4;
+
+    proj_log_debug (P, "pj_ellipsoid - final: a=%.3f f=1/%7.3f, errno=%d",
+                        P->a,  P->f!=0? 1/P->f: 0,  proj_errno (P));
+    proj_log_debug (P, "pj_ellipsoid - final: %s %s %s %s",
+                        P->def_size?           P->def_size: empty,
+                        P->def_shape?          P->def_shape: empty,
+                        P->def_spherification? P->def_spherification: empty,
+                        P->def_ellps?          P->def_ellps: empty            );
+
+    if (proj_errno (P))
+        return 5;
+
+    /* success */
+    return proj_errno_restore (P, err);
+}
+
+
+/***************************************************************************************/
+static int ellps_ellps (PJ *P) {
+/***************************************************************************************/
+    PJ B;
+    PJ_ELLPS *ellps;
+    paralist *par = 0;
+    char *name;
+    int err;
+
+    /* Sail home if ellps=xxx is not specified */
+    par = pj_get_param (P->params, "ellps");
+    if (0==par)
+        return 0;
+
+    /* Otherwise produce a fake PJ to make ellps_size/ellps_shape do the hard work for us */
+
+    /* First move B into P's context to get error messages onto the right channel */
+    B.ctx = P->ctx;
+
+    /* Then look up the right size and shape parameters from the builtin list */
+    if (strlen (par->param) < 7)
+        return proj_errno_set (P, PJD_ERR_INVALID_ARG);
+    name = par->param + 6;
+    ellps = pj_find_ellps (name);
+    if (0==ellps)
+        return proj_errno_set (P, PJD_ERR_UNKNOWN_ELLP_PARAM);
+
+    /* Now, get things ready for ellps_size/ellps_shape, make them do their thing, and clean up */
+    err = proj_errno_reset (P);
+    B = *P;
+    pj_erase_ellipsoid_def (&B);
+    B.params = pj_mkparam (ellps->major);
+    B.params->next = pj_mkparam (ellps->ell);
+
+    ellps_size (&B);
+    ellps_shape (&B);
+
+    pj_dealloc (B.params->next);
+    pj_dealloc (B.params);
+    if (proj_errno (&B))
+        return proj_errno (&B);
+
+    /* Finally update P and sail home */
+    pj_inherit_ellipsoid_def (&B, P);
+    P->def_ellps = par->param;
+    par->used = 1;
+
+    return proj_errno_restore (P, err);
+}
+
+
+/***************************************************************************************/
+static int ellps_size (PJ *P) {
+/***************************************************************************************/
+    paralist *par = 0;
+    int a_was_set = 0;
+
+    /* A size parameter *must* be given, but may have been given as ellps prior */
+    if (P->a != 0)
+        a_was_set = 1;
+
+    /* Check which size key is specified */
+    par = pj_get_param (P->params, "R");
+    if (0==par)
+        par = pj_get_param (P->params, "a");
+    if (0==par)
+        return a_was_set? 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 ellps_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);
+         if (0==P->rf)
+             return proj_errno_set (P, PJD_ERR_REV_FLATTENING_IS_ZERO);
+         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_ECCENTRICITY_IS_ONE);
+        if (P->b==P->a)
+            break;
+        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;
+}
+
+
+/* series coefficients for calculating ellipsoid-equivalent spheres */
+static const double SIXTH = 1/6.;
+static const double RA4   = 17/360.;
+static const double RA6   = 67/3024.;
+static const double RV4   = 5/72.;
+static const double RV6   = 55/1296.;
+
+/***************************************************************************************/
+static int ellps_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->f = 0;
+    P->rf = HUGE_VAL;
+    P->b = P->a;
+    pj_calc_ellipsoid_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;
+}
+
+
+/**************************************************************************************/
+void pj_erase_ellipsoid_def (PJ *P) {
+/***************************************************************************************
+    Erase all ellipsoidal parameters in P
+***************************************************************************************/
+    PJ B;
+
+    /* Make a blank PJ to copy from */
+    memset (&B, 0, sizeof (B));
+
+    /* And use it to overwrite all existing ellipsoid defs */
+    pj_inherit_ellipsoid_def (&B, P);
+}
+
+
+/**************************************************************************************/
+void pj_inherit_ellipsoid_def (const PJ *src, PJ *dst) {
+/***************************************************************************************
+    Brute force copy the ellipsoidal parameters from src to dst.  This code was
+    written before the actual ellipsoid setup parameters were kept available in
+    the PJ->def_xxx elements.
+***************************************************************************************/
 
     /* The linear parameters */
     dst->a  = src->a;
@@ -44,13 +489,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) {
+
+/***************************************************************************************/
+int pj_calc_ellipsoid_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 +535,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 +548,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 +564,49 @@ int pj_calc_ellps_params(PJ *P, double a, double es) {
     return 0;
 }
 
-/* initialize geographic shape parameters */
-int pj_ell_set(projCtx ctx, paralist *pl, double *a, double *es) {
+
+
+#ifndef KEEP_ORIGINAL_PJ_ELL_SET
+/**************************************************************************************/
+int pj_ell_set (PJ_CONTEXT *ctx, paralist *pl, double *a, double *es) {
+/***************************************************************************************
+    Initialize ellipsoidal parameters by emulating the original ellipsoid setup
+    function by Gerald Evenden, through a call to pj_ellipsoid
+***************************************************************************************/
+    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 +636,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 +659,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 +711,4 @@ bomb:
             { pj_ctx_set_errno( ctx, -13); return 1; }
     return 0;
 }
+#endif
diff --git a/src/pj_factors.c b/src/pj_factors.c
index 69d59e49..31c0e539 100644
--- a/src/pj_factors.c
+++ b/src/pj_factors.c
@@ -11,9 +11,11 @@
 
 #define EPS 1.0e-12
 
-int pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
+int pj_factors(LP lp, const PJ *P, double h, struct FACTORS *fac) {
     double cosphi, t, n, r;
     int err;
+    PJ_COORD coo = {{0, 0, 0, 0}};
+    coo.lp = lp;
 
     err = proj_errno_reset (P);
 
@@ -38,7 +40,7 @@ int pj_factors(LP lp, PJ *P, double h, struct FACTORS *fac) {
 
     /* If input latitudes are geocentric, convert to geographic */
     if (P->geoc)
-        lp.phi = atan(P->rone_es * tan(lp.phi));
+        lp = proj_geoc_lat (P, PJ_INV, coo).lp;
 
     /* If latitude + one step overshoots the pole, move it slightly inside, */
     /* so the numerical derivative still exists */
diff --git a/src/pj_init.c b/src/pj_init.c
index 534f0827..2d588ee4 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -219,23 +219,30 @@ get_opt(projCtx ctx, paralist **start, PAFile fid, char *name, paralist *next,
             strncpy(sword+1, start_of_word, word_len);
             sword[word_len+1] = '\0';
 
-            /* do not override existing parameter value of same name - unless in pipeline definition */
+            /* do not override existing parameter value of same name */
             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
-                        && !pj_param(ctx, *start, "tf").i) )
-                {
-                    next = next->next = pj_mkparam(sword+1);
+
+                /* don't default ellipse if datum, ellps or any earth model information is set */
+                if (0==strncmp(sword,"tellps=", 7)) {
+                    int n = 0;
+
+                    n += pj_param(ctx, *start, "tdatum").i;
+                    n += pj_param(ctx, *start, "tellps").i;
+                    n += pj_param(ctx, *start, "ta").i;
+                    n += pj_param(ctx, *start, "tb").i;
+                    n += pj_param(ctx, *start, "trf").i;
+                    n += pj_param(ctx, *start, "tf").i;
+                    n += pj_param(ctx, *start, "te").i;
+                    n += pj_param(ctx, *start, "tes").i;
+
+                    if (0==n)
+                        next = next->next = pj_mkparam(sword+1);
                 }
                 else
                     next = next->next = pj_mkparam(sword+1);
             }
+
+
         }
         else
         {
@@ -434,6 +441,7 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     PJ *(*proj)(PJ *);
     paralist *curr;
     int i;
+    int err;
     int found_def = 0;
     PJ *PIN = 0;
     int n_pipelines = 0;
@@ -533,17 +541,17 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 
     /* set datum parameters */
     if (pj_datum_set(ctx, start, PIN))
-        return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
+        return pj_default_destructor (PIN, proj_errno(PIN));
 
     if (PIN->need_ellps) {
-        if (pj_ell_set(ctx, start, &PIN->a, &PIN->es)) {
+        int ret = pj_ellipsoid (PIN);
+        if (0 != ret) {
             pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
-            return pj_default_destructor (PIN, PJD_ERR_MISSING_ARGS);
+            return pj_default_destructor (PIN, proj_errno(PIN));
         }
-
         PIN->a_orig = PIN->a;
         PIN->es_orig = PIN->es;
-        if (pj_calc_ellps_params(PIN, PIN->a, PIN->es))
+        if (pj_calc_ellipsoid_params (PIN, PIN->a, PIN->es))
             return pj_default_destructor (PIN, PJD_ERR_ECCENTRICITY_IS_ONE);
     }
 
@@ -628,12 +636,22 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
         s = pj_units[i].to_meter;
     }
     if (s || (s = pj_param(ctx, start, "sto_meter").s)) {
-        PIN->to_meter = pj_strtod(s, &s);
-        if (*s == '/') /* ratio number */
-            PIN->to_meter /= pj_strtod(++s, 0);
-        if (PIN->to_meter <= 0.0)
+        double factor;
+        int ratio = 0;
+
+        /* ratio number? */
+        if (*s == '/') {
+            ratio = 1;
+            s++;
+        }
+
+        factor = pj_strtod(s, &s);
+        if ((factor <= 0.0) || (1/factor==0))
             return pj_default_destructor (PIN, PJD_ERR_UNIT_FACTOR_LESS_THAN_0);
-        PIN->fr_meter = 1. / PIN->to_meter;
+
+        PIN->to_meter = ratio?  1 / factor: factor;
+        PIN->fr_meter = 1 / PIN->to_meter;
+
     } else
         PIN->to_meter = PIN->fr_meter = 1.;
 
@@ -691,12 +709,13 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     geod_init(PIN->geod, PIN->a,  (1 - sqrt (1 - PIN->es)));
 
     /* projection specific initialization */
+    err = proj_errno_reset (PIN);
     PIN = proj(PIN);
-    if ((0==PIN) || ctx->last_errno)
-    {
+    if (proj_errno (PIN)) {
         pj_free(PIN);
         return 0;
     }
+    proj_errno_restore (PIN, err);
     return PIN;
 }
 
@@ -706,8 +725,8 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
 /*      This is the application callable entry point for destroying     */
 /*      a projection definition.  It does work generic to all           */
 /*      projection types, and then calls the projection specific        */
-/*      free function (P->pfree()) to do local work.  In most cases     */
-/*      P->pfree()==pj_default_destructor.                              */
+/*      free function, P->destructor(), to do local work.               */
+/*      In most cases P->destructor()==pj_default_destructor.           */
 /************************************************************************/
 
 void pj_free(PJ *P) {
diff --git a/src/pj_list.h b/src/pj_list.h
index 0720b456..09a66034 100644
--- a/src/pj_list.h
+++ b/src/pj_list.h
@@ -21,6 +21,7 @@ PROJ_HEAD(calcofi, "Cal Coop Ocean Fish Invest Lines/Stations")
 PROJ_HEAD(cart,    "Geodetic/cartesian conversions")
 PROJ_HEAD(cass, "Cassini")
 PROJ_HEAD(cc, "Central Cylindrical")
+PROJ_HEAD(ccon, "Central Conic")
 PROJ_HEAD(cea, "Equal Area Cylindrical")
 PROJ_HEAD(chamb, "Chamberlin Trimetric")
 PROJ_HEAD(collg, "Collignon")
@@ -42,6 +43,7 @@ PROJ_HEAD(fahey, "Fahey")
 PROJ_HEAD(fouc, "Foucaut")
 PROJ_HEAD(fouc_s, "Foucaut Sinusoidal")
 PROJ_HEAD(gall, "Gall (Gall Stereographic)")
+PROJ_HEAD(geoc, "Geocentric Latitude")
 PROJ_HEAD(geocent, "Geocentric")
 PROJ_HEAD(geos, "Geostationary Satellite View")
 PROJ_HEAD(gins8, "Ginsburg VIII (TsNIIGAiK)")
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/pj_release.c b/src/pj_release.c
index 1758a536..41d02a8b 100644
--- a/src/pj_release.c
+++ b/src/pj_release.c
@@ -2,7 +2,7 @@
 
 #include <projects.h>
 
-char const pj_release[]="Rel. 4.9.3, 15 August 2016";
+char const pj_release[]="Rel. 5.0.0, development version";
 
 const char *pj_get_release()
 
diff --git a/src/pj_strerrno.c b/src/pj_strerrno.c
index fe137b87..3e726895 100644
--- a/src/pj_strerrno.c
+++ b/src/pj_strerrno.c
@@ -63,6 +63,7 @@ pj_err_list[] = {
     "lat_0 = 0",                                                       /* -55 */
     "ellipsoidal usage unsupported",                                   /* -56 */
     "only one +init allowed for non-pipeline operations",              /* -57 */
+    "argument not numerical or out of range",                          /* -58 */
 
     /* When adding error messages, remember to update ID defines in
        projects.h, and transient_error array in pj_transform                  */
@@ -91,7 +92,7 @@ char *pj_strerrno(int err) {
     adjusted_err = - err - 1;
     if (adjusted_err < (sizeof(pj_err_list) / sizeof(char *)))
         return(pj_err_list[adjusted_err]);
-    
+
     sprintf( note, "invalid projection system error (%d)", (err > -9999)? err: -9999);
     return note;
 }
diff --git a/src/pj_transform.c b/src/pj_transform.c
index 9f532188..21861331 100644
--- a/src/pj_transform.c
+++ b/src/pj_transform.c
@@ -79,7 +79,7 @@ static const int transient_error[60] = {
     /* 20 to 29 */ 1, 0, 0, 0, 0, 0, 0, 1, 0, 0,
     /* 30 to 39 */ 0, 0, 0, 0, 0, 0, 0, 0, 0, 0,
     /* 40 to 49 */ 0, 0, 0, 0, 0, 0, 0, 0, 1, 0,
-    /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 0, 0, 0 };
+    /* 50 to 59 */ 1, 0, 1, 0, 1, 1, 1, 1, 0, 0 };
 
 /************************************************************************/
 /*                            pj_transform()                            */
diff --git a/src/proj.def b/src/proj.def
index dd4a456c..887f4719 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -133,10 +133,10 @@ EXPORTS
 
     proj_torad                     @122
     proj_todeg                     @123
-    proj_rtodms                    @124
-    proj_dmstor                    @125
+    proj_geoc_lat                  @124
+    proj_rtodms                    @125
+    proj_dmstor                    @126
 
-    proj_derivatives               @126
     proj_factors                   @127
 
     proj_list_operations           @128
@@ -146,3 +146,6 @@ EXPORTS
 
     proj_angular_input             @132
     proj_angular_output            @133
+
+    pj_ellipsoid                   @134
+    pj_calc_ellipsoid_params       @135
diff --git a/src/proj.h b/src/proj.h
index cb689621..aa87ae7d 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -135,40 +135,15 @@
 extern "C" {
 #endif
 
-/************************************************************************
-    These version numbers should be updated with every release!
-    The format of PJ_VERSION is
-
-    * Before version 4.10.0: PJ_VERSION=MNP
-      where M, N, and P are the major, minor, and patch numbers;
-      e.g., PJ_VERSION=493 for version 4.9.3.
-
-    * Version 4.10.0 and later: PJ_VERSION=MMMNNNPP later
-      where MMM, NNN, PP are the major, minor, and patch numbers.
-      The minor and patch numbers are padded with leading zeros if
-      necessary);
-      e.g., PJ_VERSION=401000 for version 4.10.0.
-************************************************************************/
-#define PROJ_VERSION_MAJOR 10
+/* The version numbers should be updated with every release! **/
+#define PROJ_VERSION_MAJOR 5
 #define PROJ_VERSION_MINOR 0
 #define PROJ_VERSION_PATCH 0
 
-#ifndef PJ_VERSION
-#define PJ_VERSION 10##000##00
-#endif
-#ifndef PROJ_VERSION
-#define PROJ_VERSION PJ_VERSION
-#endif
-
-
 extern char const pj_release[]; /* global release id string */
 
 /* first forward declare everything needed */
 
-/* Data type for generic geodetic 3D data */
-union PJ_TRIPLET;
-typedef union PJ_TRIPLET PJ_TRIPLET;
-
 /* Data type for generic geodetic 3D data plus epoch information */
 union PJ_COORD;
 typedef union PJ_COORD PJ_COORD;
@@ -176,11 +151,24 @@ typedef union PJ_COORD PJ_COORD;
 struct PJ_AREA;
 typedef struct PJ_AREA PJ_AREA;
 
-struct DERIVS;
-typedef struct DERIVS PJ_DERIVS;
-
-struct FACTORS;
-typedef struct FACTORS PJ_FACTORS;
+/* The slimmed down PROJ 5.0.0 version of struct FACTORS  */
+/* Will take over the world and the name when we can rid  */
+/* the library for deprecated stuff, but it's the typedef */
+/* which is userspace useful, so it does not do much of a */
+/* difference */
+struct P5_FACTORS {                  /* Common designation */
+    double meridional_scale;               /* h */
+    double parallel_scale;                 /* k */
+    double areal_scale;                    /* s */
+
+    double angular_distortion;             /* omega */
+    double meridian_parallel_angle;        /* theta-prime */
+    double meridian_convergence;           /* alpha */
+
+    double tissot_semimajor;               /* a */
+    double tissot_semiminor;               /* b */
+};
+typedef struct P5_FACTORS PJ_FACTORS;
 
 /* Data type for projection/transformation information */
 struct PJconsts;
@@ -213,17 +201,11 @@ struct PJ_PRIME_MERIDIANS;
 typedef struct PJ_PRIME_MERIDIANS PJ_PRIME_MERIDIANS;
 
 
-/* Omega, Phi, Kappa: Rotations */
-typedef struct {double o, p, k;}  PJ_OPK;
-
-/* Easting, Northing, and some kind of height (orthometric or ellipsoidal) */
-typedef struct {double e, n, h;}  PJ_ENH;
-
-/* Geodetic spatiotemporal coordinate types */
+/* Geodetic, mostly spatiotemporal coordinate types */
 typedef struct { double   x,   y,  z, t; }  PJ_XYZT;
-typedef struct { double   e,   n,  h, t; }  PJ_ENHT;
 typedef struct { double   u,   v,  w, t; }  PJ_UVWT;
 typedef struct { double lam, phi,  z, t; }  PJ_LPZT;
+typedef struct { double o, p, k; }          PJ_OPK;  /* Rotations: omega, phi, kappa */
 
 /* Classic proj.4 pair/triplet types */
 typedef struct { double   u,   v; }  UV;
@@ -234,60 +216,22 @@ typedef struct { double   x,   y,  z; }  XYZ;
 typedef struct { double   u,   v,  w; }  UVW;
 typedef struct { double lam, phi,  z; }  LPZ;
 
-/* Ancillary pairs and triplets for geodetic computations */
-
-/* easting and northing */
-typedef struct { double   e, n; }  PJ_EN;
-
-/* Degrees, minutes, and seconds */
-typedef struct { double d, m,  s; }  PJ_DMS;
-
-/* Geoid undulation (N) and deflections of the vertical (eta, zeta) */
-typedef struct { double   e,  z, N; }  PJ_EZN;
-
-/* Ellipsoidal parameters */
-typedef struct { double   a,   f; }  PJ_AF;
 
-/* Avoid preprocessor renaming and implicit type-punning: Use unions to make it explicit */
+/* Avoid preprocessor renaming and implicit type-punning: Use a union to make it explicit */
 union PJ_COORD {
+     double v[4];   /* First and foremost, it really is "just 4 numbers in a vector" */
     PJ_XYZT xyzt;
     PJ_UVWT uvwt;
-    PJ_ENHT enht;
     PJ_LPZT lpzt;
-    PJ_ENH  enh;
-    PJ_EN   en;
-    double v[4]; /* It's just a vector */
-    XYZ  xyz;
-    UVW  uvw;
-    LPZ  lpz;
-    XY   xy;
-    UV   uv;
-    LP   lp;
+     PJ_OPK opk;
+        XYZ xyz;
+        UVW uvw;
+        LPZ lpz;
+         XY xy;
+         UV uv;
+         LP lp;
 };
 
-union PJ_TRIPLET {
-    PJ_OPK  opk;
-    PJ_ENH  enh;
-    PJ_EZN  ezn;
-    PJ_DMS  dms;
-    double v[3]; /* It's just a vector */
-    XYZ    xyz;
-    LPZ    lpz;
-    UVW    uvw;
-    XY     xy;
-    LP     lp;
-    UV     uv;
-    PJ_AF  af;
-};
-
-union PJ_PAIR {
-    XY     xy;
-    LP     lp;
-    UV     uv;
-    PJ_AF  af;
-    PJ_EN  en;
-    double v[2]; /* Yes - It's really just a vector! */
-};
 
 struct PJ_INFO {
     char        release[64];        /* Release info. Version + date         */
@@ -341,7 +285,6 @@ PJ_CONTEXT *proj_context_create (void);
 PJ_CONTEXT *proj_context_destroy (PJ_CONTEXT *ctx);
 
 
-
 /* Manage the transformation definition object PJ */
 PJ  *proj_create (PJ_CONTEXT *ctx, const char *definition);
 PJ  *proj_create_argv (PJ_CONTEXT *ctx, int argc, char **argv);
@@ -363,15 +306,12 @@ enum PJ_DIRECTION {
 typedef enum PJ_DIRECTION PJ_DIRECTION;
 
 
-
-
 int proj_angular_input (PJ *P, enum PJ_DIRECTION dir);
 int proj_angular_output (PJ *P, enum PJ_DIRECTION dir);
 
 
 PJ_COORD proj_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coord);
-
-
+int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord);
 size_t proj_trans_generic (
     PJ *P,
     PJ_DIRECTION direction,
@@ -381,7 +321,6 @@ size_t proj_trans_generic (
     double *t, size_t st, size_t nt
 );
 
-int proj_trans_array (PJ *P, PJ_DIRECTION direction, size_t n, PJ_COORD *coord);
 
 /* Initializers */
 PJ_COORD proj_coord (double x, double y, double z, double t);
@@ -403,14 +342,13 @@ 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_reset (PJ *P);
-void proj_errno_restore (PJ *P, int err);
-
+int  proj_errno (const PJ *P);
+int  proj_errno_set (const PJ *P, int err);
+int  proj_errno_reset (const PJ *P);
+int  proj_errno_restore (const PJ *P, int err);
 
-PJ_DERIVS  proj_derivatives(PJ *P, const LP lp);
-PJ_FACTORS proj_factors(PJ *P, const LP lp);
+/* Scaling and angular distortion factors */
+PJ_FACTORS proj_factors(PJ *P, LP lp);
 
 /* Info functions - get information about various PROJ.4 entities */
 PJ_INFO      proj_info(void);
@@ -431,6 +369,9 @@ const PJ_PRIME_MERIDIANS  *proj_list_prime_meridians(void);
 double proj_torad (double angle_in_degrees);
 double proj_todeg (double angle_in_radians);
 
+/* Geographical to geocentric latitude - another of the "simple, but useful" */
+PJ_COORD proj_geoc_lat (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
+
 double proj_dmstor(const char *is, char **rs);
 char*  proj_rtodms(char *s, double r, int pos, int neg);
 
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 6e89fb60..9a576997 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -334,6 +334,48 @@ size_t proj_trans_generic (
 
 
 /*************************************************************************************/
+PJ_COORD proj_geoc_lat (const PJ *P, PJ_DIRECTION direction, PJ_COORD coo) {
+/**************************************************************************************
+    Convert geographical latitude to geocentric (or the other way round if
+    direction = PJ_INV)
+
+    The conversion involves a call to the tangent function, which goes through the
+    roof at the poles, so very close (the last few micrometers) to the poles no
+    conversion takes place and the input latitude is copied directly to the output.
+
+    Fortunately, the geocentric latitude converges to the geographical at the
+    poles, so the difference is negligible.
+
+    For the spherical case, the geographical latitude equals the geocentric, and
+    consequently, the input is copied directly to the output.
+**************************************************************************************/
+    const double limit = M_HALFPI - 1e-12;
+    PJ_COORD res = coo;
+    if ((coo.lp.phi > limit) || (coo.lp.phi < -limit) || (P->es==0))
+        return res;
+    if (direction==PJ_FWD)
+        res.lp.phi = atan (P->one_es * tan (coo.lp.phi) );
+    else
+        res.lp.phi = atan (P->rone_es * tan (coo.lp.phi) );
+
+    return res;
+}
+
+double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);}
+double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);}
+
+double proj_dmstor(const char *is, char **rs) {
+    return dmstor(is, rs);
+}
+
+char*  proj_rtodms(char *s, double r, int pos, int neg) {
+    return rtodms(s, r, pos, neg);
+}
+
+
+
+
+/*************************************************************************************/
 PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
 /**************************************************************************************
     Create a new PJ object in the context ctx, using the given definition. If ctx==0,
@@ -390,7 +432,6 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
     argv = (char **) calloc (argc, sizeof (char *));
     if (0==argv)
         return pj_dealloc (args);
-
     argv[0] = args;
     for (i = 0, j = 1;  i < n;  i++) {
         if (0==args[i])
@@ -474,46 +515,42 @@ PJ *proj_destroy (PJ *P) {
     return 0;
 }
 
-int proj_errno (PJ *P) {
-    return pj_ctx_get_errno (pj_get_ctx (P));
+int proj_errno (const PJ *P) {
+    return pj_ctx_get_errno (pj_get_ctx ((PJ *) P));
 }
 
 /*****************************************************************************/
-void proj_errno_set (PJ *P, int err) {
+int proj_errno_set (const 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);
+    proj_context_errno_set (pj_get_ctx ((PJ *) P), err);
     errno = err;
-    return;
+    return err;
 }
 
 /*****************************************************************************/
-void proj_errno_restore (PJ *P, int err) {
+int proj_errno_restore (const 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;
 }
 
 /*****************************************************************************/
-int proj_errno_reset (PJ *P) {
+int proj_errno_reset (const PJ *P) {
 /******************************************************************************
     Clears errno in the context and thread local levels
     through the low level pj_ctx interface.
@@ -521,23 +558,30 @@ 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;
     last_errno = proj_errno (P);
 
-    pj_ctx_set_errno (pj_get_ctx (P), 0);
+    pj_ctx_set_errno (pj_get_ctx ((PJ *) P), 0);
     errno = 0;
     return last_errno;
 }
@@ -782,28 +826,9 @@ PJ_INIT_INFO proj_init_info(const char *initname){
 }
 
 
-/*****************************************************************************/
-PJ_DERIVS proj_derivatives(PJ *P, const LP lp) {
-/******************************************************************************
-    Derivatives of coordinates.
-
-    returns PJ_DERIVS. If unsuccessfull error number is set and the returned
-    struct contains NULL data.
-
-******************************************************************************/
-    PJ_DERIVS derivs;
-
-    if (pj_deriv(lp, 1e-5, P, &derivs)) {
-        /* errno set in pj_derivs */
-        memset(&derivs, 0, sizeof(PJ_DERIVS));
-    }
-
-    return derivs;
-}
-
 
 /*****************************************************************************/
-PJ_FACTORS proj_factors(PJ *P, const LP lp) {
+PJ_FACTORS proj_factors(PJ *P, LP lp) {
 /******************************************************************************
     Cartographic characteristics at point lp.
 
@@ -814,12 +839,22 @@ PJ_FACTORS proj_factors(PJ *P, const LP lp) {
     struct contains NULL data.
 
 ******************************************************************************/
-    PJ_FACTORS factors;
+    PJ_FACTORS factors = {0,0,0, 0,0,0, 0,0};
+    struct FACTORS f;
 
-    if (pj_factors(lp, P, 0.0, &factors)) {
-        /* errno set in pj_factors */
-        memset(&factors, 0, sizeof(PJ_FACTORS));
-    }
+    if (pj_factors(lp, P, 0.0, &f))
+        return factors;
+
+    factors.meridional_scale  =  f.h;
+    factors.parallel_scale    =  f.k;
+    factors.areal_scale       =  f.s;
+
+    factors.angular_distortion        =  f.omega;
+    factors.meridian_parallel_angle   =  f.thetap;
+    factors.meridian_convergence      =  f.conv;
+
+    factors.tissot_semimajor  =  f.a;
+    factors.tissot_semiminor  =  f.b;
 
     return factors;
 }
@@ -841,14 +876,3 @@ const PJ_PRIME_MERIDIANS *proj_list_prime_meridians(void) {
     return pj_get_prime_meridians_ref();
 }
 
-double proj_torad (double angle_in_degrees) { return PJ_TORAD (angle_in_degrees);}
-double proj_todeg (double angle_in_radians) { return PJ_TODEG (angle_in_radians);}
-
-
-double proj_dmstor(const char *is, char **rs) {
-    return dmstor(is, rs);
-}
-
-char*  proj_rtodms(char *s, double r, int pos, int neg) {
-    return rtodms(s, r, pos, neg);
-}
diff --git a/src/proj_api.h b/src/proj_api.h
index b0c4b71f..b57d9f38 100644
--- a/src/proj_api.h
+++ b/src/proj_api.h
@@ -28,19 +28,14 @@
 
 
  /*
-  * This version number should be updated with every release!  The format of
-  * PJ_VERSION is
+  * This version number should be updated with every release!
   *
-  * * Before version 4.10.0: PJ_VERSION=MNP where M, N, and P are the major,
-  *   minor, and patch numbers; e.g., PJ_VERSION=493 for version 4.9.3.
+  * This file is expected to be removed from the PROJ distribution
+  * when a few minor-version releases has been made.
   *
-  * * Version 4.10.0 and later: PJ_VERSION=MMMNNNPP later where MMM, NNN, PP
-  *   are the major, minor, and patch numbers (the minor and patch numbers
-  *   are padded with leading zeros if necessary); e.g., PJ_VERSION=401000
-  *   for version 4.10.0.
   */
  #ifndef PJ_VERSION
- #define PJ_VERSION 493
+ #define PJ_VERSION 500
  #endif
 
 
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 95f34994..fd6dc75d 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -101,11 +101,11 @@ typedef void (*PJ_LOG_FUNCTION)(void *, int, const char *);
 void proj_log_error (PJ *P, const char *fmt, ...);
 void proj_log_debug (PJ *P, const char *fmt, ...);
 void proj_log_trace (PJ *P, const char *fmt, ...);
-/*void proj_log_func (PJ *P, void *app_data, void (*log)(void *, int, const char *));*/
 void proj_log_func (PJ_CONTEXT *ctx, void *app_data, PJ_LOG_FUNCTION log);
 
-void pj_inherit_ellipsoid_defs(const PJ *src, PJ *dst);
-int pj_calc_ellps_params(PJ *P, double a, double es);
+void pj_inherit_ellipsoid_def (const PJ *src, PJ *dst);
+void pj_erase_ellipsoid_def (PJ *P);
+int pj_calc_ellipsoid_params (PJ *P, double a, double es);
 
 /* Lowest level: Minimum support for fileapi */
 void proj_fileapi_set (PJ *P, void *fileapi);
diff --git a/src/projects.h b/src/projects.h
index 32b74b08..99df0e3e 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -178,7 +178,6 @@ union  PJ_COORD;
 struct geod_geodesic;
 struct pj_opaque;
 struct ARG_list;
-struct FACTORS;
 struct PJ_REGION_S;
 typedef struct PJ_REGION_S  PJ_Region;
 typedef struct ARG_list paralist;   /* parameter list */
@@ -229,6 +228,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 +417,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 */
@@ -479,14 +478,14 @@ enum deprecated_constants_for_now_dropped_analytical_factors {
 /* library errors */
 #define PJD_ERR_NO_ARGS                  -1
 #define PJD_ERR_NO_OPTION_IN_INIT_FILE   -2
-#define PJD_ERR_NO_COLOR_IN_INIT_STRING  -3
+#define PJD_ERR_NO_COLON_IN_INIT_STRING  -3
 #define PJD_ERR_PROJ_NOT_NAMED           -4
 #define PJD_ERR_UNKNOWN_PROJECTION_ID    -5
 #define PJD_ERR_ECCENTRICITY_IS_ONE      -6
 #define PJD_ERR_UNKNOW_UNIT_ID           -7
 #define PJD_ERR_INVALID_BOOLEAN_PARAM    -8
 #define PJD_ERR_UNKNOWN_ELLP_PARAM       -9
-#define PJD_ERR_REC_FLATTENING_IS_ZERO  -10
+#define PJD_ERR_REV_FLATTENING_IS_ZERO  -10
 #define PJD_ERR_REF_RAD_LARGER_THAN_90  -11
 #define PJD_ERR_ES_LESS_THAN_ZERO       -12
 #define PJD_ERR_MAJOR_AXIS_NOT_GIVEN    -13
@@ -534,6 +533,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 +676,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 *);
@@ -702,8 +703,8 @@ double  pj_authlat(double, double *);
 COMPLEX pj_zpoly1(COMPLEX, COMPLEX *, int);
 COMPLEX pj_zpolyd1(COMPLEX, COMPLEX *, int, COMPLEX *);
 
-int pj_deriv(LP, double, PJ *, struct DERIVS *);
-int pj_factors(LP, PJ *, double, struct FACTORS *);
+int pj_deriv(LP, double, const PJ *, struct DERIVS *);
+int pj_factors(LP, const PJ *, double, struct FACTORS *);
 
 struct PW_COEF {    /* row coefficient structure */
     int m;          /* number of c coefficients (=0 for none) */
diff --git a/test/gie/builtins.gie b/test/gie/builtins.gie
index c6a801f8..6b820c01 100644
--- a/test/gie/builtins.gie
+++ b/test/gie/builtins.gie
@@ -502,6 +502,37 @@ expect  -0.001790493 0.000895247
 accept  -200 -100
 expect  -0.001790493 -0.000895247
 
+===============================================================================
+Central Conic
+	Sph
+	lat_1 
+===============================================================================
+
+-------------------------------------------------------------------------------
+operation +proj=pipeline +step +proj=ccon +lat_1=52 +lat_0=52 +lon_0=19 +a=6390000 +x_0=330000 +y_0=-350000 +step +proj=axisswap +order=1,-2,3,4
+-------------------------------------------------------------------------------
+tolerance 0.00010 mm
+accept 24 55
+expect 650031.54109413219363 4106.1617770643609028
+accept 15 49
+expect 37074.189007307473069 676826.23559270039774
+accept 24 49
+expect 696053.36061617843913 672294.56795827199940
+accept 19 52
+expect 330000.00000000000000 350000.00000000000000
+
+direction inverse
+accept 0 0
+expect 13.840227318521004431 55.030403993648806391
+accept 0 700000
+expect 14.514453594615022781 48.773847834747808675
+accept 700000 0
+expect 24.782707184271129766 55.003515505218481835
+accept 700000 700000
+expect 24.027610763560529927 48.750476070495021286
+accept 330000 350000
+expect 19.000000000000000000 52.000000000000000000
+
 
 ===============================================================================
 Central Cylindrical
diff --git a/test/gie/ellipsoid.gie b/test/gie/ellipsoid.gie
new file mode 100644
index 00000000..65101f6c
--- /dev/null
+++ b/test/gie/ellipsoid.gie
@@ -0,0 +1,161 @@
+===============================================================================
+
+Test pj_ellipsoid, the reimplementation of pj_ell_set
+
+===============================================================================
+
+
+BEGIN
+
+-------------------------------------------------------------------------------
+First a spherical example
+-------------------------------------------------------------------------------
+operation proj=merc R=6400000
+-------------------------------------------------------------------------------
+tolerance 10 nm
+accept    1 2
+expect    111701.0721276371  223447.5262032605
+
+accept    12 55
+expect    1340412.8655316452  7387101.1430967357
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Then an explicitly defined ellipsoidal example
+-------------------------------------------------------------------------------
+operation proj=merc a=6400000 rf=297
+-------------------------------------------------------------------------------
+tolerance 10 nm
+accept    1 2
+expect    111701.0721276371  221945.9681832088
+
+accept    12 55
+expect    1340412.8655316452  7351803.9151705895
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Then try using a built in ellipsoid
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80
+-------------------------------------------------------------------------------
+tolerance 10 nm
+accept    1 2
+expect    111319.4907932736  221194.0771604237
+
+accept    12 55
+expect    1335833.8895192828  7326837.7148738774
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Then try to fail deliberately
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80000000000
+expect    failure errno unknown_ellp_param
+operation proj=merc +a=-1
+expect    failure errno major_axis_not_given
+operation proj=merc +no_defs
+expect    failure errno major_axis_not_given
+
+# This one should succeed due to ellps=WGS84 in proj_def.dat
+operation proj=merc
+accept    0 0
+expect    0 0
+
+operation proj=merc +es=-1
+expect    failure errno major_axis_not_given
+
+operation
+expect    failure
+operation cobra
+expect    failure
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+Finally test the spherification functionality
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_A
+tolerance 10 nm
+accept    12 55
+expect    1334340.6237297705  7353636.6296552019
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_V
+tolerance 10 nm
+accept    12 55
+expect    1334339.2852675652  7353629.2533042720
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_a
+tolerance 10 nm
+accept    12 55
+expect    1333594.4904527504  7349524.6413825499
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_g
+tolerance 10 nm
+accept    12 55
+expect    1333592.6102291327  7349514.2793497816
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_h
+tolerance 10 nm
+accept    12 55
+expect    1333590.7300081658  7349503.9173316229
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_lat_a=60
+tolerance 10 nm
+accept    12 55
+expect    1338073.7436268919  7374210.0924803326
+-------------------------------------------------------------------------------
+operation proj=merc ellps=GRS80 R_lat_g=60
+tolerance 10 nm
+accept    12 55
+expect    1338073.2696101593  7374207.4801437631
+-------------------------------------------------------------------------------
+
+
+-------------------------------------------------------------------------------
+This one from testvarious failed at first version of the pull request
+-------------------------------------------------------------------------------
+operation  proj=healpix a=1 lon_0=0 ellps=WGS84
+-------------------------------------------------------------------------------
+accept     0 41.937853904844985
+expect     0 0.78452
+accept    -90 0
+expect    -1.56904	0
+-------------------------------------------------------------------------------
+
+-------------------------------------------------------------------------------
+Shape parameters
+-------------------------------------------------------------------------------
+operation  proj=utm zone=32   ellps=GRS80 rf=0
+expect     failure errno rev_flattening_is_zero
+
+operation  proj=utm zone=32   ellps=GRS80 es=1
+expect     failure errno eccentricity_is_one
+
+operation  proj=utm zone=32   ellps=GRS80 b=0
+expect     failure errno eccentricity_is_one
+
+operation  proj=utm zone=32   ellps=GRS80 b=6000000
+accept     12 55
+expect     699293.0880   5674591.5295
+
+operation  proj=utm zone=32   ellps=GRS80 rf=300
+accept     12 55
+expect     691873.1212   6099054.9661
+
+operation  proj=utm zone=32   ellps=GRS80 f=0.00333333333333
+accept     12 55
+expect     691873.1212   6099054.9661
+
+operation  proj=utm zone=32   ellps=GRS80 b=6000000
+accept     12 55
+expect     699293.0880   5674591.5295
+
+operation  proj=utm zone=32   a=6400000 b=6000000
+accept     12 55
+expect     700416.5900   5669475.8884
+-------------------------------------------------------------------------------
+
+END
diff --git a/test/gie/more_builtins.gie b/test/gie/more_builtins.gie
index 1e3ce185..8ac8a667 100644
--- a/test/gie/more_builtins.gie
+++ b/test/gie/more_builtins.gie
@@ -1,4 +1,3 @@
-
 ===============================================================================
 
 Various test material, mostly converted from selftest entries in PJ_xxx.c
@@ -251,14 +250,47 @@ accept     3370658.378    711877.314    5349787.086    2017.0
 expect     3370658.18890  711877.42370  5349787.12430  2017.0
 accept     3370658.378    711877.314    5349787.086    2018.0
 expect     3370658.18087  711877.42750  5349787.12648  2018.0
+-------------------------------------------------------------------------------
 
 
 
 -------------------------------------------------------------------------------
-builtins
+geocentric latitude
 -------------------------------------------------------------------------------
+operation  proj=geoc ellps=GRS80
+accept     12   55                  0 0
+expect     12   54.818973308324573  0 0
+roundtrip  1000
 
+accept     12   90                  0 0
+expect     12   90                  0 0
 
+accept     12  -90                  0 0
+expect     12  -90                  0 0
 
-END
+accept     12   89.99999999999      0 0
+expect     12   89.999999999989996  0 0
+-------------------------------------------------------------------------------
 
+
+-------------------------------------------------------------------------------
+some less used options
+-------------------------------------------------------------------------------
+operation  proj=utm  ellps=GRS80  zone=32  to_meter=0
+expect     failure
+
+operation  proj=utm  ellps=GRS80  zone=32  to_meter=10
+accept     12   55
+expect        69187.5632    609890.7825
+-------------------------------------------------------------------------------
+
+
+
+-------------------------------------------------------------------------------
+run the few gie-builtin tests, which are currently either awkward or impossible
+to express in the gie command set
+-------------------------------------------------------------------------------
+builtins
+-------------------------------------------------------------------------------
+
+END
diff --git a/travis/install.sh b/travis/install.sh
index 049669d4..f71108c6 100755
--- a/travis/install.sh
+++ b/travis/install.sh
@@ -84,6 +84,7 @@ PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/builtins.gie
 PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/more_builtins.gie
 PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/deformation.gie
 PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/axisswap.gie
+PROJ_LIB=$GRIDDIR ./src/gie ./test/gie/ellipsoid.gie
 
 # install & run the working GIGS test
   # create locations that pyproj understands