Merge pull request #786 from busstoptaktik/interdependent-wrap-ups

Interdependent wrap ups

5 files changed, 51 insertions, 46 deletions

diff --git a/src/gie.c b/src/gie.c
index 9931b2ef..95f50708 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -681,7 +681,7 @@ back/forward transformation pairs.
             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) F->lineno);
-        fprintf (T.fout, "     roundtrip deviation: %.3f mm, expected: %.3f mm\n", 1000*r, 1000*d);
+        fprintf (T.fout, "     roundtrip deviation: %.6f mm, expected: %.6f mm\n", 1000*r, 1000*d);
     }
     return another_failure ();
 }
@@ -700,13 +700,13 @@ static int expect_message (double d, const char *args) {
 
     fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) F->lineno);
     fprintf (T.fout, "     expected: %s\n", args);
-    fprintf (T.fout, "     got:      %.9f   %.9f", T.b.xy.x,  T.b.xy.y);
+    fprintf (T.fout, "     got:      %.12f   %.12f", T.b.xy.x,  T.b.xy.y);
     if (T.b.xyzt.t!=0 || T.b.xyzt.z!=0)
         fprintf (T.fout, "   %.9f", T.b.xyz.z);
     if (T.b.xyzt.t!=0)
         fprintf (T.fout, "   %.9f", T.b.xyzt.t);
     fprintf (T.fout, "\n");
-    fprintf (T.fout, "     deviation:  %.3f mm,  expected:  %.3f mm\n", 1000*d, 1000*T.tolerance);
+    fprintf (T.fout, "     deviation:  %.6f mm,  expected:  %.6f mm\n", 1000*d, 1000*T.tolerance);
     return 1;
 }
 
@@ -829,42 +829,35 @@ Tell GIE what to expect, when transforming the ACCEPTed input
     /* expected angular values, probably in degrees */
     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]);
+        printf ("EXPECTS  %.12f  %.12f  %.12f  %.12f\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.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]);
+        printf ("ACCEPTS  %.12f  %.12f  %.12f  %.12f\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 (T.P, T.dir, co): co;
     if (T.verbosity > 3)
-        printf ("GOT      %.4f  %.4f  %.4f  %.4f\n", co.v[0],co.v[1],co.v[2],co.v[3]);
+        printf ("GOT      %.12f  %.12f  %.12f  %.12f\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)) {
-        double e = HUGE_VAL;
-        d = proj_lpz_dist (T.P, ce.lpz, co.lpz);
-        /* check whether input was already in radians */
-        if (d > T.tolerance)
-            e = proj_lpz_dist (T.P, T.e.lpz, co.lpz);
-        if (e < d)
-            d = e;
-
-        /* or the tolerance may be based on euclidean distance */
-        if (d > T.tolerance)
-            e = proj_xyz_dist (T.b.xyz, T.e.xyz);
-        if (e < d)
-            d = e;
+#if 0
+    /* We need to handle unusual axis orders - that'll be an item for version 5.1 */
+    if (T.P->axisswap) {
+        ce = proj_trans (T.P->axisswap, T.dir, ce);
+        co = proj_trans (T.P->axisswap, T.dir, co);
     }
+#endif
+    if (proj_angular_output (T.P, T.dir))
+        d = proj_lpz_dist (T.P, ce.lpz, co.lpz);
     else
-        d = proj_xyz_dist (T.b.xyz, T.e.xyz);
+        d = proj_xyz_dist (co.xyz, ce.xyz);
+
     if (d > T.tolerance)
         return expect_message (d, args);
 
     another_success ();
-
     return 0;
 }
 
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 66b86aab..d238fa4d 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -73,11 +73,13 @@ static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD coo) {
             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->helmert,    PJ_INV, coo); /* Step into local frame */
             coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
         }
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
         if (P->vgridshift)
-            coo = proj_trans (P->vgridshift, PJ_FWD, coo);
+            coo = proj_trans (P->vgridshift, PJ_FWD, coo); /* Go orthometric from geometric */
 
         /* Distance from central meridian, taking system zero meridian into account */
         coo.lp.lam = (coo.lp.lam - P->from_greenwich) - P->lam0;
@@ -92,7 +94,7 @@ static PJ_COORD pj_fwd_prepare (PJ *P, PJ_COORD 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 proj_trans (P->helmert, PJ_INV, coo);
     return coo;
 }
 
@@ -139,14 +141,18 @@ static PJ_COORD pj_fwd_finalize (PJ *P, PJ_COORD coo) {
             coo.lpz.lam = adjlon(coo.lpz.lam);
 
         if (P->vgridshift)
-            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+            coo = proj_trans (P->vgridshift, PJ_FWD, coo); /* Go orthometric from geometric */
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
         if (P->hgridshift)
-            coo = proj_trans (P->hgridshift, PJ_FWD, coo);
+            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->helmert,    PJ_INV, coo); /* Step into local frame */
             coo = proj_trans (P->cart,       PJ_INV, coo); /* Go back to angular using local ellps */
         }
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
 
         /* If input latitude was geocentrical, convert back to geocentrical */
         if (P->geoc)
diff --git a/src/pj_init.c b/src/pj_init.c
index 0acc6c71..f54a2a92 100644
--- a/src/pj_init.c
+++ b/src/pj_init.c
@@ -593,12 +593,12 @@ pj_init_ctx(projCtx ctx, int argc, char **argv) {
     if (pj_datum_set(ctx, start, PIN))
         return pj_default_destructor (PIN, proj_errno(PIN));
 
-    if (PIN->need_ellps) {
-        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, proj_errno(PIN));
-        }
+    err = pj_ellipsoid (PIN);
+    if (PIN->need_ellps && 0 != err) {
+        pj_log (ctx, PJ_LOG_DEBUG_MINOR, "pj_init_ctx: Must specify ellipsoid or sphere");
+        return pj_default_destructor (PIN, proj_errno(PIN));
+    }
+    if (0==err) {
         PIN->a_orig = PIN->a;
         PIN->es_orig = PIN->es;
         if (pj_calc_ellipsoid_params (PIN, PIN->a, PIN->es))
diff --git a/src/pj_inv.c b/src/pj_inv.c
index 5e44eb34..4ea88b69 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -80,11 +80,13 @@ static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
             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->helmert,    PJ_FWD, coo); /* Step into WGS84 */
             coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
         }
+        if (coo.lp.lam==HUGE_VAL)
+            return coo;
         if (P->vgridshift)
-            coo = proj_trans (P->vgridshift, PJ_INV, coo);
+            coo = proj_trans (P->vgridshift, PJ_INV, coo); /* Go geometric from orthometric */
         return coo;
     }
 
@@ -99,9 +101,8 @@ static PJ_COORD pj_inv_prepare (PJ *P, PJ_COORD coo) {
         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) {
+        if (P->is_geocent)
             coo = proj_trans (P->cart, PJ_INV, coo);
-        }
 
         return coo;
 
@@ -148,14 +149,18 @@ static PJ_COORD pj_inv_finalize (PJ *P, PJ_COORD coo) {
                 coo.lpz.lam = adjlon(coo.lpz.lam);
 
             if (P->vgridshift)
-                coo = proj_trans (P->vgridshift, PJ_INV, coo);
+                coo = proj_trans (P->vgridshift, PJ_INV, coo); /* Go geometric from orthometric */
+            if (coo.lp.lam==HUGE_VAL)
+                return 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->helmert,    PJ_FWD, coo); /* Step into WGS84 */
                 coo = proj_trans (P->cart_wgs84, PJ_INV, coo); /* Go back to angular using WGS84 ellps */
             }
+            if (coo.lp.lam==HUGE_VAL)
+                return coo;
         }
 
         /* If input latitude was geocentrical, convert back to geocentrical */
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 42e3cbf9..8f6305a4 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -473,12 +473,13 @@ invocators can emulate the behaviour of pj_transform and the cs2cs app.
         def = malloc (100+n);
         if (0==def)
             return 0;
-        sprintf (def, "break_cs2cs_recursion     proj=helmert %s", s);
+        sprintf (def, "break_cs2cs_recursion     proj=helmert %s transpose", s);
         Q = proj_create (P->ctx, def);
+        pj_inherit_ellipsoid_def (P, Q);
         free (def);
         if (0==Q)
             return 0;
-        P->helmert = skip_prep_fin(Q);
+        P->helmert = skip_prep_fin (Q);
 
         break;
     }
@@ -491,14 +492,14 @@ invocators can emulate the behaviour of pj_transform and the cs2cs app.
         Q = proj_create (P->ctx, def);
         if (0==Q)
             return 0;
-        pj_inherit_ellipsoid_def(P, Q);
-        P->cart = skip_prep_fin(Q);
+        pj_inherit_ellipsoid_def (P, Q);
+        P->cart = skip_prep_fin (Q);
 
         sprintf (def, "break_cs2cs_recursion     proj=cart  ellps=WGS84");
         Q = proj_create (P->ctx, def);
         if (0==Q)
             return 0;
-        P->cart_wgs84 = skip_prep_fin(Q);
+        P->cart_wgs84 = skip_prep_fin (Q);
     }
 
     return 1;