Merge pull request #731 from busstoptaktik/4D-API_cs2cs-style

[WIP] 4D API cs2cs style

20 files changed, 454 insertions, 123 deletions

diff --git a/src/PJ_axisswap.c b/src/PJ_axisswap.c
index d57834fc..44446d9c 100644
--- a/src/PJ_axisswap.c
+++ b/src/PJ_axisswap.c
@@ -274,8 +274,8 @@ PJ *CONVERSION(axisswap,0) {
     }
 
     if (pj_param(P->ctx, P->params, "tangularunits").i) {
-        P->left  = PJ_IO_UNITS_RADIANS;
-        P->right = PJ_IO_UNITS_RADIANS;
+        P->left  = PJ_IO_UNITS_ANGULAR;
+        P->right = PJ_IO_UNITS_ANGULAR;
     } else {
         P->left  = PJ_IO_UNITS_PROJECTED;
         P->right = PJ_IO_UNITS_PROJECTED;
diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index 8e11c030..0746ec08 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -213,7 +213,7 @@ PJ *CONVERSION(cart,1) {
     P->inv3d  =  geodetic;
     P->fwd    =  cart_forward;
     P->inv    =  cart_reverse;
-    P->left   =  PJ_IO_UNITS_RADIANS;
+    P->left   =  PJ_IO_UNITS_ANGULAR;
     P->right  =  PJ_IO_UNITS_CARTESIAN;
     return P;
 }
diff --git a/src/PJ_geoc.c b/src/PJ_geoc.c
index 865b1089..20064f65 100644
--- a/src/PJ_geoc.c
+++ b/src/PJ_geoc.c
@@ -48,8 +48,8 @@ static PJ *CONVERSION(geoc, 1) {
     P->inv4d = inverse;
     P->fwd4d = forward;
 
-    P->left   =  PJ_IO_UNITS_RADIANS;
-    P->right  =  PJ_IO_UNITS_RADIANS;
+    P->left   =  PJ_IO_UNITS_ANGULAR;
+    P->right  =  PJ_IO_UNITS_ANGULAR;
 
     P->is_latlong = 1;
     return P;
diff --git a/src/PJ_helmert.c b/src/PJ_helmert.c
index f726d946..229e30c2 100644
--- a/src/PJ_helmert.c
+++ b/src/PJ_helmert.c
@@ -486,6 +486,25 @@ PJ *TRANSFORMATION(helmert, 0) {
         P->right = PJ_IO_UNITS_PROJECTED;
     }
 
+    /* Support the classic PROJ towgs84 parameter, but allow later overrides.*/
+    /* Note that if towgs84 is specified, the datum_params array is set up   */
+    /* for us automagically by the pj_datum_set call in pj_init_ctx */
+    if (pj_param_exists (P->params, "towgs84")) {
+        Q->xyz_0.x = P->datum_params[0];
+        Q->xyz_0.y = P->datum_params[1];
+        Q->xyz_0.z = P->datum_params[2];
+
+        Q->opk_0.o = P->datum_params[3];
+        Q->opk_0.p = P->datum_params[4];
+        Q->opk_0.k = P->datum_params[5];
+
+        /* We must undo conversion to absolute scale from pj_datum_set */
+        if (0==P->datum_params[6])
+            Q->scale_0 = 0;
+        else
+            Q->scale_0 = (P->datum_params[6] - 1) * 1e6;
+    }
+
     /* Translations */
     if (pj_param (P->ctx, P->params, "tx").i)
         Q->xyz_0.x = pj_param (P->ctx, P->params, "dx").f;
@@ -570,14 +589,13 @@ PJ *TRANSFORMATION(helmert, 0) {
     /* Let's help with debugging */
     if (proj_log_level(P->ctx, PJ_LOG_TELL) >= PJ_LOG_DEBUG) {
         proj_log_debug(P, "Helmert parameters:");
-        proj_log_debug(P, "x=  % 3.5f  y=  % 3.5f  z=  % 3.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
-        proj_log_debug(P, "rx= % 3.5f  ry= % 3.5f  rz= % 3.5f",
+        proj_log_debug(P, "x=  %8.5f  y=  %8.5f  z=  %8.5f", Q->xyz.x, Q->xyz.y, Q->xyz.z);
+        proj_log_debug(P, "rx= %8.5f  ry= %8.5f  rz= %8.5f",
                 Q->opk.o / ARCSEC_TO_RAD, Q->opk.p / ARCSEC_TO_RAD, Q->opk.k / ARCSEC_TO_RAD);
-        proj_log_debug(P, "s=% 3.5f  exact=% d  transpose=% d",
-                Q->scale, Q->exact, Q->transpose);
-        proj_log_debug(P, "dx= % 3.5f  dy= % 3.5f  dz= % 3.5f", Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
-        proj_log_debug(P, "drx=% 3.5f  dry=% 3.5f  drz=% 3.5f", Q->dopk.o, Q->dopk.p, Q->dopk.k);
-        proj_log_debug(P, "ds=% 3.5f  t_epoch=% 5.5f  t_obs=% 5.5f", Q->dscale, Q->t_epoch, Q->t_obs);
+        proj_log_debug(P, "s=  %8.5f  exact=%d  transpose=%d", Q->scale, Q->exact, Q->transpose);
+        proj_log_debug(P, "dx= %8.5f  dy= %8.5f  dz= %8.5f",   Q->dxyz.x, Q->dxyz.y, Q->dxyz.z);
+        proj_log_debug(P, "drx=%8.5f  dry=%8.5f  drz=%8.5f",   Q->dopk.o, Q->dopk.p, Q->dopk.k);
+        proj_log_debug(P, "ds= %8.5f  t_epoch=%8.5f  t_obs=%8.5f", Q->dscale, Q->t_epoch, Q->t_obs);
     }
 
     if ((Q->opk.o==0)  && (Q->opk.p==0)  && (Q->opk.k==0) && (Q->scale==0) &&
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 5c2b944d..54440822 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -54,8 +54,8 @@ PJ *TRANSFORMATION(hgridshift,0) {
     P->fwd    = 0;
     P->inv    = 0;
 
-    P->left  = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
+    P->left  = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_ANGULAR;
 
     if (0==pj_param(P->ctx, P->params, "tgrids").i) {
         proj_log_error(P, "hgridshift: +grids parameter missing.");
diff --git a/src/PJ_latlong.c b/src/PJ_latlong.c
index b1909954..1331d59a 100644
--- a/src/PJ_latlong.c
+++ b/src/PJ_latlong.c
@@ -98,8 +98,8 @@ static PJ *latlong_setup (PJ *P) {
     P->fwd3d = latlong_forward_3d;
     P->inv4d = latlong_inverse_4d;
     P->fwd4d = latlong_forward_4d;
-    P->left = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
+    P->left = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_ANGULAR;
     return P;
 }
 
diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c
index 765e1d50..b537e802 100644
--- a/src/PJ_molodensky.c
+++ b/src/PJ_molodensky.c
@@ -281,8 +281,8 @@ PJ *TRANSFORMATION(molodensky,1) {
     P->fwd    = forward_2d;
     P->inv    = reverse_2d;
 
-    P->left   = PJ_IO_UNITS_RADIANS;
-    P->right  = PJ_IO_UNITS_RADIANS;
+    P->left   = PJ_IO_UNITS_ANGULAR;
+    P->right  = PJ_IO_UNITS_ANGULAR;
 
     /* read args */
     if (pj_param(P->ctx, P->params, "tdx").i)
diff --git a/src/PJ_ob_tran.c b/src/PJ_ob_tran.c
index debb211c..c447ac08 100644
--- a/src/PJ_ob_tran.c
+++ b/src/PJ_ob_tran.c
@@ -232,7 +232,7 @@ PJ *PROJECTION(ob_tran) {
 
     /* Support some rather speculative test cases, where the rotated projection */
     /* is actually latlong. We do not want scaling in that case... */
-    if (Q->link->right==PJ_IO_UNITS_RADIANS)
+    if (Q->link->right==PJ_IO_UNITS_ANGULAR)
         P->right = PJ_IO_UNITS_PROJECTED;
 
 
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index 850734fc..bb8b4d4d 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -69,8 +69,8 @@ PJ *TRANSFORMATION(vgridshift,0) {
     P->fwd    = 0;
     P->inv    = 0;
 
-    P->left  = PJ_IO_UNITS_RADIANS;
-    P->right = PJ_IO_UNITS_RADIANS;
+    P->left  = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_ANGULAR;
 
     return P;
 }
diff --git a/src/gie.c b/src/gie.c
index e0045ebd..b2673d13 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -387,6 +387,7 @@ Interpret <args> as a numeric followed by a linear decadal prefix.
 Return the properly scaled numeric
 ******************************************************************************/
     double s;
+    const double GRS80_DEG = 111319.4908; /* deg-to-m at equator of GRS80 */
     const char *endp = args;
     s = proj_strtod (args, (char **) &endp);
     if (args==endp)
@@ -408,6 +409,10 @@ Return the properly scaled numeric
         s /= 1e6;
     else if (0==strcmp(endp, "nm"))
         s /= 1e9;
+    else if (0==strcmp(endp, "rad"))
+        s = GRS80_DEG * proj_todeg (s);
+    else if (0==strcmp(endp, "deg"))
+        s = GRS80_DEG * s;
     else
         s *= default_scale;
     return s;
@@ -553,19 +558,29 @@ using the "builtins" command verb.
 }
 
 
-static PJ_COORD torad_coord (PJ_COORD a) {
-    PJ_COORD c = a;
-    c.lpz.lam = proj_torad (a.lpz.lam);
-    c.lpz.phi = proj_torad (a.lpz.phi);
-    return c;
+static PJ_COORD torad_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    size_t i, n;
+    char *axis = "enut";
+    paralist *l = pj_param_exists (P->params, "axis");
+    if (l && dir==PJ_INV)
+        axis = l->param + strlen ("axis=");
+    for (i = 0,  n = strlen (axis);  i < n;  i++)
+        if (strchr ("news", axis[i]))
+            a.v[i] = proj_torad (a.v[i]);
+    return a;
 }
 
 
-static PJ_COORD todeg_coord (PJ_COORD a) {
-    PJ_COORD c = a;
-    c.lpz.lam = proj_todeg (a.lpz.lam);
-    c.lpz.phi = proj_todeg (a.lpz.phi);
-    return c;
+static PJ_COORD todeg_coord (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
+    size_t i, n;
+    char *axis = "enut";
+    paralist *l = pj_param_exists (P->params, "axis");
+    if (l && dir==PJ_FWD)
+        axis = l->param + strlen ("axis=");
+    for (i = 0,  n = strlen (axis);  i < n;  i++)
+        if (strchr ("news", axis[i]))
+            a.v[i] = proj_todeg (a.v[i]);
+    return a;
 }
 
 
@@ -627,7 +642,7 @@ back/forward transformation pairs.
     d = d==HUGE_VAL?  T.tolerance:  d;
 
     /* input ("accepted") values - probably in degrees */
-    coo = proj_angular_input  (T.P, T.dir)? torad_coord (T.a):  T.a;
+    coo = proj_angular_input  (T.P, T.dir)? torad_coord (T.P, T.dir, T.a):  T.a;
 
     r = proj_roundtrip (T.P, T.dir, ntrips, &coo);
     if (r <= d)
@@ -747,7 +762,7 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         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;
+        ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
         co = expect_trans_n_dim (ci);
 
         /* Failed to fail? - that's a failure */
@@ -780,20 +795,20 @@ Tell GIE what to expect, when transforming the ACCEPTed input
         return expect_message_cannot_parse (args);
 
     /* expected angular values, probably in degrees */
-    ce = proj_angular_output (T.P, T.dir)? torad_coord (T.e): T.e;
+    ce = proj_angular_output (T.P, T.dir)? torad_coord (T.P, T.dir, T.e): T.e;
     if (T.verbosity > 3)
         printf ("EXPECTS  %.4f  %.4f  %.4f  %.4f\n", ce.v[0],ce.v[1],ce.v[2],ce.v[3]);
 
     /* input ("accepted") values, also probably in degrees */
-    ci = proj_angular_input (T.P, T.dir)? torad_coord (T.a): T.a;
+    ci = proj_angular_input (T.P, T.dir)? torad_coord (T.P, T.dir, T.a): T.a;
     if (T.verbosity > 3)
         printf ("ACCEPTS  %.4f  %.4f  %.4f  %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
 
     /* angular output from proj_trans comes in radians */
     co = expect_trans_n_dim (ci);
-    T.b = proj_angular_output (T.P, T.dir)? todeg_coord (co): co;
+    T.b = proj_angular_output (T.P, T.dir)? todeg_coord (T.P, T.dir, co): co;
     if (T.verbosity > 3)
-        printf ("GOT      %.4f  %.4f  %.4f  %.4f\n", ci.v[0],ci.v[1],ci.v[2],ci.v[3]);
+        printf ("GOT      %.4f  %.4f  %.4f  %.4f\n", co.v[0],co.v[1],co.v[2],co.v[3]);
 
     /* but there are a few more possible input conventions... */
     if (proj_angular_output (T.P, T.dir)) {
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index a1025aac..0509858a 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -32,58 +32,137 @@
 #include "proj_internal.h"
 #include "projects.h"
 
-PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
+#define INPUT_UNITS  P->left
+#define OUTPUT_UNITS P->right
+
+
+static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
     if (HUGE_VAL==coo.v[0])
         return proj_coord_error ();
 
+    /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
+    if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
+    if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
+
     /* Check validity of angular input coordinates */
-    if (P->left==PJ_IO_UNITS_RADIANS) {
+    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
         double t;
-        LP lp = coo.lp;
 
         /* check for latitude or longitude over-range */
-        t = (lp.phi < 0  ?  -lp.phi  :  lp.phi) - M_HALFPI;
-        if (t > PJ_EPS_LAT  ||  lp.lam > 10  ||  lp.lam < -10) {
+        t = (coo.lp.phi < 0  ?  -coo.lp.phi  :  coo.lp.phi) - M_HALFPI;
+        if (t > PJ_EPS_LAT  ||  coo.lp.lam > 10  ||  coo.lp.lam < -10) {
             proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
             return proj_coord_error ();
         }
 
         /* Clamp latitude to -90..90 degree range */
-        if (lp.phi > M_HALFPI)
-            lp.phi = M_HALFPI;
-        if (lp.phi < -M_HALFPI)
-            lp.phi = -M_HALFPI;
+        if (coo.lp.phi > M_HALFPI)
+            coo.lp.phi = M_HALFPI;
+        if (coo.lp.phi < -M_HALFPI)
+            coo.lp.phi = -M_HALFPI;
 
         /* If input latitude is geocentrical, convert to geographical */
         if (P->geoc)
             coo = proj_geoc_lat (P, PJ_INV, coo);
 
+        /* Ensure longitude is in the -pi:pi range */
+        if (0==P->over)
+            coo.lp.lam = adjlon(coo.lp.lam);
+
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_INV, coo);
+        else if (P->helmert) {
+            coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
+            coo = proj_trans (P->helmert,    PJ_FWD, coo); /* Step into local frame */
+            coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
+        }
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_FWD, coo);
+
         /* Distance from central meridian, taking system zero meridian into account */
-        lp.lam = (lp.lam - P->from_greenwich) - P->lam0;
+        coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
 
         /* Ensure longitude is in the -pi:pi range */
         if (0==P->over)
-            lp.lam = adjlon(lp.lam);
+            coo.lp.lam = adjlon(coo.lp.lam);
 
-        coo.lp = lp;
+        return coo;
     }
 
+
+    /* We do not support gridshifts on cartesian input */
+    if (INPUT_UNITS==PJ_IO_UNITS_CARTESIAN && P->helmert)
+            return proj_trans (P->helmert, PJ_FWD, coo);
     return coo;
 }
 
 
-PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
+
+static PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
+
+    switch (OUTPUT_UNITS) {
+
+    /* Handle false eastings/northings and non-metric linear units */
+    case PJ_IO_UNITS_CARTESIAN:
+
+        if (P->is_geocent) {
+            coo = proj_trans (P->cart, PJ_FWD, coo);
+        }
+
+        coo.xyz.x = P->fr_meter * (coo.xyz.x + P->x0);
+        coo.xyz.y = P->fr_meter * (coo.xyz.y + P->y0);
+        coo.xyz.z = P->fr_meter * (coo.xyz.z + P->z0);
+        break;
 
     /* Classic proj.4 functions return plane coordinates in units of the semimajor axis */
-    if (P->right==PJ_IO_UNITS_CLASSIC) {
+    case PJ_IO_UNITS_CLASSIC:
         coo.xy.x *= P->a;
         coo.xy.y *= P->a;
+
+    /* Falls through */ /* (<-- GCC warning silencer) */
+    /* to continue processing in common with PJ_IO_UNITS_PROJECTED */
+    case PJ_IO_UNITS_PROJECTED:
+        coo.xyz.x = P->fr_meter  * (coo.xyz.x + P->x0);
+        coo.xyz.y = P->fr_meter  * (coo.xyz.y + P->y0);
+        coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+        break;
+
+    case PJ_IO_UNITS_WHATEVER:
+        break;
+
+    case PJ_IO_UNITS_ANGULAR:
+        if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR)
+            break;
+
+        /* adjust longitude to central meridian */
+        if (0==P->over)
+            coo.lpz.lam = adjlon(coo.lpz.lam);
+
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+        else if (P->helmert) {
+            coo = proj_trans (P->cart_wgs84, PJ_FWD, coo); /* Go cartesian in WGS84 frame */
+            coo = proj_trans (P->helmert,    PJ_FWD, coo); /* Step into local frame */
+            coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
+        }
+
+        /* If input latitude was geocentrical, convert back to geocentrical */
+        if (P->geoc)
+            coo = proj_geoc_lat (P, PJ_FWD, coo);
+
+
+        /* Distance from central meridian, taking system zero meridian into account */
+        coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
+
+        /* adjust longitude to central meridian */
+        if (0==P->over)
+            coo.lpz.lam = adjlon(coo.lpz.lam);
     }
 
-    /* Handle false eastings/northings and non-metric linear units */
-    coo.xyz.x = P->fr_meter  * (coo.xyz.x + P->x0);
-    coo.xyz.y = P->fr_meter  * (coo.xyz.y + P->y0);
-    coo.xyz.z = P->vfr_meter * (coo.xyz.z + P->z0);
+    if (P->axisswap)
+        coo = proj_trans (P->axisswap, PJ_FWD, coo);
 
     return coo;
 }
diff --git a/src/pj_geocent.c b/src/pj_geocent.c
index eca62080..6530b103 100644
--- a/src/pj_geocent.c
+++ b/src/pj_geocent.c
@@ -54,6 +54,8 @@ PJ *PROJECTION(geocent) {
     P->y0 = 0.0;
     P->inv = inverse;
     P->fwd = forward;
+    P->left = PJ_IO_UNITS_ANGULAR;
+    P->right = PJ_IO_UNITS_CARTESIAN;
 
     return P;
 }
diff --git a/src/pj_init.c b/src/pj_init.c
index f25c0344..51a742c7 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -769,26 +769,3 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     proj_errno_restore (PIN, err);
     return PIN;
 }
-
-
-
-/************************************************************************/
-/*                              pj_free()                               */
-/*                                                                      */
-/*      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->destructor(), to do local work.               */
-/*      In most cases P->destructor()==pj_default_destructor.           */
-/************************************************************************/
-
-void pj_free(PJ *P) {
-    if (0==P)
-        return;
-    /* free projection parameters - all the hard work is done by */
-    /* pj_default_destructor (in pj_malloc.c), which is supposed */
-    /* to be called as the last step of the local destructor     */
-    /* pointed to by P->destructor. In most cases,               */
-    /* pj_default_destructor actually *is* what is pointed to    */
-    P->destructor (P, 0);
-}
diff --git a/src/pj_inv.c b/src/pj_inv.c
index e1fb05f6..45e43820 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -32,46 +32,131 @@
 #include "proj_internal.h"
 #include "projects.h"
 
+#define INPUT_UNITS  P->right
+#define OUTPUT_UNITS P->left
 
-
-PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
+static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
     if (coo.xyz.x == HUGE_VAL) {
         proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
         return proj_coord_error ();
     }
 
+    /* The helmert datum shift will choke unless it gets a sensible 4D coordinate */
+    if (HUGE_VAL==coo.v[2] && P->helmert) coo.v[2] = 0.0;
+    if (HUGE_VAL==coo.v[3] && P->helmert) coo.v[3] = 0.0;
+
+    if (P->axisswap)
+        coo = proj_trans (P->axisswap, PJ_INV, coo);
+
+    /* Check validity of angular input coordinates */
+    if (INPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+        double t;
+
+        /* check for latitude or longitude over-range */
+        t = (coo.lp.phi < 0  ?  -coo.lp.phi  :  coo.lp.phi) - M_HALFPI;
+        if (t > PJ_EPS_LAT  ||  coo.lp.lam > 10  ||  coo.lp.lam < -10) {
+            proj_errno_set (P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
+            return proj_coord_error ();
+        }
+
+        /* Clamp latitude to -90..90 degree range */
+        if (coo.lp.phi > M_HALFPI)
+            coo.lp.phi = M_HALFPI;
+        if (coo.lp.phi < -M_HALFPI)
+            coo.lp.phi = -M_HALFPI;
+
+        /* If input latitude is geocentrical, convert to geographical */
+        if (P->geoc)
+            coo = proj_geoc_lat (P, PJ_INV, coo);
+
+        /* Distance from central meridian, taking system zero meridian into account */
+        coo.lp.lam = (coo.lp.lam + P->from_greenwich) - P->lam0;
+
+        /* Ensure longitude is in the -pi:pi range */
+        if (0==P->over)
+            coo.lp.lam = adjlon(coo.lp.lam);
+
+        if (P->hgridshift)
+            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+        else if (P->helmert) {
+            coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
+            coo = proj_trans (P->helmert,    PJ_INV, coo); /* Step into WGS84 */
+            coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
+        }
+        if (P->vgridshift)
+            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+        return coo;
+    }
+
+    /* Handle remaining possible input types */
+    switch (INPUT_UNITS) {
+    case PJ_IO_UNITS_WHATEVER:
+        return coo;
+
     /* de-scale and de-offset */
-    coo.xyz.x = (coo.xyz.x * P->to_meter  - P->x0);
-    coo.xyz.y = (coo.xyz.y * P->to_meter  - P->y0);
-    coo.xyz.z = (coo.xyz.z * P->vto_meter - P->z0);
-
-    /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis  */
-    /* Multiplying by ra, rather than dividing by a because the CALCOFI projection       */
-    /* stomps on a and hence (apparently) depends on this to roundtrip correctly         */
-    /* (CALCOFI avoids further scaling by stomping - a better solution must be possible) */
-    if (P->right==PJ_IO_UNITS_CLASSIC) {
+    case PJ_IO_UNITS_CARTESIAN:
+        coo.xyz.x = P->to_meter * coo.xyz.x - P->x0;
+        coo.xyz.y = P->to_meter * coo.xyz.y - P->y0;
+        coo.xyz.z = P->to_meter * coo.xyz.z - P->z0;
+
+        if (P->is_geocent) {
+            coo = proj_trans (P->cart, PJ_INV, coo);
+        }
+
+        return coo;
+
+    case PJ_IO_UNITS_PROJECTED:
+    case PJ_IO_UNITS_CLASSIC:
+        coo.xyz.x = P->to_meter  * coo.xyz.x - P->x0;
+        coo.xyz.y = P->to_meter  * coo.xyz.y - P->y0;
+        coo.xyz.z = P->vto_meter * coo.xyz.z - P->z0;
+        if (INPUT_UNITS==PJ_IO_UNITS_PROJECTED)
+            return coo;
+
+        /* Classic proj.4 functions expect plane coordinates in units of the semimajor axis  */
+        /* Multiplying by ra, rather than dividing by a because the CALCOFI projection       */
+        /* stomps on a and hence (apparently) depends on this to roundtrip correctly         */
+        /* (CALCOFI avoids further scaling by stomping - but a better solution is possible)  */
         coo.xyz.x *= P->ra;
         coo.xyz.y *= P->ra;
+        return coo;
+    /* Silence some compiler warnings about PJ_IO_UNITS_ANGULAR not handled */
+    default:
+        break;
     }
 
+    /* Should not happen, so we could return pj_coord_err here */
     return coo;
 }
 
 
 
-PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
+static PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
     if (coo.xyz.x == HUGE_VAL) {
         proj_errno_set (P, PJD_ERR_INVALID_X_OR_Y);
         return proj_coord_error ();
     }
 
-    if (P->left==PJ_IO_UNITS_RADIANS) {
-        /* Distance from central meridian, taking system zero meridian into account */
-        coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
-
-         /* adjust longitude to central meridian */
-        if (0==P->over)
-            coo.lpz.lam = adjlon(coo.lpz.lam);
+    if (OUTPUT_UNITS==PJ_IO_UNITS_ANGULAR) {
+
+        if (INPUT_UNITS!=PJ_IO_UNITS_ANGULAR) {
+            /* Distance from central meridian, taking system zero meridian into account */
+            coo.lp.lam = coo.lp.lam + P->from_greenwich + P->lam0;
+
+            /* adjust longitude to central meridian */
+            if (0==P->over)
+                coo.lpz.lam = adjlon(coo.lpz.lam);
+
+            if (P->vgridshift)
+                coo = proj_trans (P->vgridshift, PJ_INV, coo);
+            if (P->hgridshift)
+                coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+            else if (P->helmert) {
+                coo = proj_trans (P->cart,       PJ_FWD, coo); /* Go cartesian in local frame */
+                coo = proj_trans (P->helmert,    PJ_INV, coo); /* Step into WGS84 */
+                coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
+            }
+        }
 
         /* If input latitude was geocentrical, convert back to geocentrical */
         if (P->geoc)
diff --git a/src/pj_malloc.c b/src/pj_malloc.c
index 63278a56..cfb7f658 100644
--- a/src/pj_malloc.c
+++ b/src/pj_malloc.c
@@ -162,6 +162,32 @@ void *pj_dealloc_params (PJ_CONTEXT *ctx, paralist *start, int errlev) {
 }
 
 
+
+
+/************************************************************************/
+/*                              pj_free()                               */
+/*                                                                      */
+/*      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->destructor(), to do local work.               */
+/*      In most cases P->destructor()==pj_default_destructor.           */
+/************************************************************************/
+
+void pj_free(PJ *P) {
+    if (0==P)
+        return;
+    /* free projection parameters - all the hard work is done by */
+    /* pj_default_destructor, which is supposed */
+    /* to be called as the last step of the local destructor     */
+    /* pointed to by P->destructor. In most cases,               */
+    /* pj_default_destructor actually *is* what is pointed to    */
+    P->destructor (P, 0);
+}
+
+
+
+
 /*****************************************************************************/
 void *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
 /*****************************************************************************
@@ -196,6 +222,14 @@ void *pj_default_destructor (PJ *P, int errlev) {   /* Destructor */
     /* free parameter list elements */
     pj_dealloc_params (pj_get_ctx(P), P->params, errlev);
 
+    /* free the cs2cs emulation elements */
+    pj_free (P->axisswap);
+    pj_free (P->helmert);
+    pj_free (P->cart);
+    pj_free (P->cart_wgs84);
+    pj_free (P->hgridshift);
+    pj_free (P->vgridshift);
+
     pj_dealloc (P->opaque);
     return pj_dealloc(P);
 }
diff --git a/src/proj.c b/src/proj.c
index 1d702ea4..1c86854e 100644
--- a/src/proj.c
+++ b/src/proj.c
@@ -152,12 +152,12 @@ static void process(FILE *fid) {
             }
         } else {    /* x-y or decimal degree ascii output, scale if warranted by output units */
             if (inverse) {
-                if (Proj->left == PJ_IO_UNITS_RADIANS) {
+                if (Proj->left == PJ_IO_UNITS_ANGULAR) {
                     data.uv.v *= RAD_TO_DEG;
                     data.uv.u *= RAD_TO_DEG;
                 }
             } else {
-                if (Proj->right == PJ_IO_UNITS_RADIANS) {
+                if (Proj->right == PJ_IO_UNITS_ANGULAR) {
                     data.uv.v *= RAD_TO_DEG;
                     data.uv.u *= RAD_TO_DEG;
                 }
@@ -262,7 +262,7 @@ static void vprocess(FILE *fid) {
         }
 
         /* apply rad->deg scaling in case the output from a pipeline has degrees as units */
-        if (!inverse && Proj->right == PJ_IO_UNITS_RADIANS) {
+        if (!inverse && Proj->right == PJ_IO_UNITS_ANGULAR) {
             dat_xy.x *= RAD_TO_DEG;
             dat_xy.y *= RAD_TO_DEG;
         }
diff --git a/src/proj.def b/src/proj.def
index 6ab0c007..099a78a5 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -155,3 +155,5 @@ EXPORTS
     pj_approx_2D_trans             @138
     pj_approx_3D_trans             @139
     pj_has_inverse                 @140
+    pj_param_exists                @141
+
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 2c6bac2d..2b250398 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -52,8 +52,8 @@ int proj_angular_input (PJ *P, enum PJ_DIRECTION dir) {
     dir: {PJ_FWD, PJ_INV}
 ******************************************************************************/
     if (PJ_FWD==dir)
-        return pj_left (P)==PJ_IO_UNITS_RADIANS;
-    return pj_right (P)==PJ_IO_UNITS_RADIANS;
+        return pj_left (P)==PJ_IO_UNITS_ANGULAR;
+    return pj_right (P)==PJ_IO_UNITS_ANGULAR;
 }
 
 /*****************************************************************************/
@@ -100,7 +100,7 @@ double proj_xyz_dist (XYZ a, XYZ b) {
 /* Measure numerical deviation after n roundtrips fwd-inv (or inv-fwd) */
 double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
     int i;
-    PJ_COORD o, u, org;
+    PJ_COORD t, org;
 
     if (0==P)
         return HUGE_VAL;
@@ -111,23 +111,23 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
     }
 
     /* in the first half-step, we generate the output value */
-    u = org = *coo;
-    o = *coo = proj_trans (P, direction, u);
+    org  = *coo;
+    *coo = proj_trans (P, direction, org);
+    t = *coo;
 
-    /* now we take n-1 full steps */
-    for (i = 0;  i < n - 1;  i++) {
-        u = proj_trans (P, -direction, o);
-        o = proj_trans (P,  direction, u);
-    }
+    /* now we take n-1 full steps in inverse direction: We are */
+    /* out of phase due to the half step already taken */
+    for (i = 0;  i < n - 1;  i++)
+        t = proj_trans (P,  direction,  proj_trans (P, -direction, t) );
 
     /* finally, we take the last half-step */
-    u = proj_trans (P, -direction, o);
+    t = proj_trans (P, -direction, t);
 
     /* checking for angular *input* since we do a roundtrip, and end where we begin */
     if (proj_angular_input (P, direction))
-        return proj_lpz_dist (P, org.lpz, u.lpz);
+        return proj_lpz_dist (P, org.lpz, t.lpz);
 
-    return proj_xyz_dist (org.xyz, u.xyz);
+    return proj_xyz_dist (org.xyz, t.xyz);
 }
 
 
@@ -380,7 +380,114 @@ char*  proj_rtodms(char *s, double r, int pos, int neg) {
     return rtodms(s, r, pos, neg);
 }
 
+/*************************************************************************************/
+static PJ* skip_prep_fin(PJ *P) {
+/**************************************************************************************
+Skip prepare and finalize function for the various "helper operations" added to P when
+in cs2cs compatibility mode.
+**************************************************************************************/
+    P->skip_fwd_prepare  = 1;
+    P->skip_fwd_finalize = 1;
+    P->skip_inv_prepare  = 1;
+    P->skip_inv_finalize = 1;
+    return P;
+}
 
+/*************************************************************************************/
+static int pj_cs2cs_emulation_setup (PJ *P) {
+/**************************************************************************************
+If any cs2cs style modifiers are given (axis=..., towgs84=..., ) create the 4D API
+equivalent operations, so the preparation and finalization steps in the pj_inv/pj_fwd
+invocators can emulate the behaviour of pj_transform and the cs2cs app.
+**************************************************************************************/
+    PJ *Q;
+    paralist *p;
+    char def[1000];
+    if (0==P)
+        return 0;
+
+    /* Don't recurse when calling proj_create (which calls us back) */
+    if (pj_param_exists (P->params, "break_cs2cs_recursion"))
+        return 1;
+
+    /* Swap axes? */
+    p = pj_param_exists (P->params, "axis");
+
+    /* Don't axisswap if data are already in "enu" order */
+    if (p && (0!=strcmp ("enu", p->param))) {
+        sprintf (def, "break_cs2cs_recursion     proj=axisswap  axis=%s", P->axis);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->axisswap = skip_prep_fin(Q);
+    }
+
+    /* Geoid grid(s) given? */
+    p = pj_param_exists (P->params, "geoidgrids");
+    if (p  &&  strlen (p->param) > strlen ("geoidgrids=")) {
+        char *gridnames = p->param + strlen ("geoidgrids=");
+        sprintf (def, "break_cs2cs_recursion     proj=vgridshift  grids=%s", gridnames);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->vgridshift = skip_prep_fin(Q);
+    }
+
+    /* Datum shift grid(s) given? */
+    p = pj_param_exists (P->params, "nadgrids");
+    if (p  &&  strlen (p->param) > strlen ("nadgrids=")) {
+        char *gridnames = p->param + strlen ("nadgrids=");
+        sprintf (def, "break_cs2cs_recursion     proj=hgridshift  grids=%s", gridnames);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->hgridshift = skip_prep_fin(Q);
+    }
+
+    /* We ignore helmert if we have grid shift */
+    p = P->hgridshift ? 0 : pj_param_exists (P->params, "towgs84");
+    while (p) {
+        char *s = p->param;
+        double *d = P->datum_params;
+        size_t n = strlen (s);
+
+        /* We ignore null helmert shifts (common in auto-translated resource files, e.g. epsg) */
+        if (0==d[0] && 0==d[1] && 0==d[2] && 0==d[3] && 0==d[4] && 0==d[5] && 0==d[6])
+            break;
+
+        if (n > 900)
+            return 0;
+        if (n <= 8) /* 8==strlen ("towgs84=") */
+            return 0;
+        sprintf (def, "break_cs2cs_recursion     proj=helmert %s", s);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->helmert = skip_prep_fin(Q);
+
+        break;
+    }
+
+    /* We also need cartesian/geographical transformations if we are working in */
+    /* geocentric/cartesian space or we need to do a Helmert transform.         */
+    if (P->is_geocent || P->helmert) {
+        char *wgs84 = "ellps=WGS84";
+        sprintf (def, "break_cs2cs_recursion     proj=cart");
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        pj_inherit_ellipsoid_def(P, Q);
+        P->cart = skip_prep_fin(Q);
+
+        sprintf (def, "break_cs2cs_recursion     proj=cart  %s", wgs84);
+        Q = proj_create (P->ctx, def);
+        if (0==Q)
+            return 0;
+        P->cart_wgs84 = skip_prep_fin(Q);
+    }
+
+    return 1;
+}
 
 
 /*************************************************************************************/
@@ -394,9 +501,10 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
     It may even use free formatting "proj  =  utm;  zone  =32  ellps= GRS80".
     Note that the semicolon separator is allowed, but not required.
 **************************************************************************************/
-    PJ   *P;
-    char *args, **argv;
+    PJ    *P;
+    char  *args, **argv;
     size_t argc, n;
+    int    ret;
 
     if (0==ctx)
         ctx = pj_get_default_ctx ();
@@ -421,6 +529,12 @@ PJ *proj_create (PJ_CONTEXT *ctx, const char *definition) {
 
     pj_dealloc (argv);
     pj_dealloc (args);
+
+    /* Support cs2cs-style modifiers */
+    ret = pj_cs2cs_emulation_setup  (P);
+    if (0==ret)
+        return proj_destroy (P);
+
     return P;
 }
 
@@ -449,6 +563,7 @@ indicator, as in {"+proj=utm", "+zone=32"}, or leave it out, as in {"proj=utm",
         return 0;
 
     P = proj_create (ctx, c);
+
     pj_dealloc ((char *) c);
     return P;
 }
diff --git a/src/proj_internal.h b/src/proj_internal.h
index 41f95c07..1b98c346 100644
--- a/src/proj_internal.h
+++ b/src/proj_internal.h
@@ -68,7 +68,7 @@ enum pj_io_units {
     PJ_IO_UNITS_CLASSIC   = 1,  /* Scaled meters (right), projected system */
     PJ_IO_UNITS_PROJECTED = 2,  /* Meters, projected system */
     PJ_IO_UNITS_CARTESIAN = 3,  /* Meters, 3D cartesian system */
-    PJ_IO_UNITS_RADIANS   = 4   /* Radians */
+    PJ_IO_UNITS_ANGULAR   = 4   /* Radians */
 };
 enum pj_io_units pj_left (PJ *P);
 enum pj_io_units pj_right (PJ *P);
@@ -84,10 +84,6 @@ void proj_context_inherit (PJ *parent, PJ *child);
 PJ_COORD pj_fwd4d (PJ_COORD coo, PJ *P);
 PJ_COORD pj_inv4d (PJ_COORD coo, PJ *P);
 
-PJ_COORD pj_fwd_prepare  (PJ *P, PJ_COORD coo);
-PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo);
-PJ_COORD pj_inv_prepare  (PJ *P, PJ_COORD coo);
-PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo);
 PJ_COORD pj_approx_2D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 PJ_COORD pj_approx_3D_trans (PJ *P, PJ_DIRECTION direction, PJ_COORD coo);
 
diff --git a/src/projects.h b/src/projects.h
index 86ca2524..e4bdca4f 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -196,7 +196,7 @@ enum pj_io_units {
     PJ_IO_UNITS_CLASSIC   = 1,  /* Scaled meters (right), projected system */
     PJ_IO_UNITS_PROJECTED = 2,  /* Meters, projected system */
     PJ_IO_UNITS_CARTESIAN = 3,  /* Meters, 3D cartesian system */
-    PJ_IO_UNITS_RADIANS   = 4   /* Radians */
+    PJ_IO_UNITS_ANGULAR   = 4   /* Radians */
 };
 #endif
 #ifndef PROJ_H
@@ -349,6 +349,14 @@ struct PJconsts {
     enum pj_io_units left;          /* Flags for input/output coordinate types */
     enum pj_io_units right;
 
+    /* These PJs are used for implementing cs2cs style coordinate handling in the 4D API */
+    PJ *axisswap;
+    PJ *cart;
+    PJ *cart_wgs84;
+    PJ *helmert;
+    PJ *hgridshift;
+    PJ *vgridshift;
+
 
     /*************************************************************************************
 
@@ -394,7 +402,7 @@ struct PJconsts {
     double  from_greenwich;            /* prime meridian offset (in radians) */
     double  long_wrap_center;          /* 0.0 for -180 to 180, actually in radians*/
     int     is_long_wrap_set;
-    char    axis[4];                   /* TODO: Description needed */
+    char    axis[4];                   /* Axis order, pj_transform/pj_adjust_axis */
 
     /* New Datum Shift Grid Catalogs */
     char   *catalog_name;
@@ -615,7 +623,7 @@ C_NAMESPACE PJ *pj_##name (PJ *P) {                          \
     P->destructor = pj_default_destructor;                   \
     P->descr = des_##name;                                   \
     P->need_ellps = NEED_ELLPS;                              \
-    P->left  = PJ_IO_UNITS_RADIANS;                          \
+    P->left  = PJ_IO_UNITS_ANGULAR;                          \
     P->right = PJ_IO_UNITS_CLASSIC;                          \
     return P;                                                \
 }                                                            \