Improved IO predicates (#648)

* enter proj_angular_input and proj_angular_output, exit proj_angular_left and proj_angular_right

* remove unused variable 'unit'

* In gie: remove unused func 'torad_if_needed', and add static keyword where needed

* In gie: add some comments

6 files changed, 215 insertions, 107 deletions

diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index d8f4a9b9..52237f24 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -567,6 +567,63 @@ int pj_cart_selftest (void) {
     for (pm_list = proj_list_prime_meridians(); pm_list->id; ++pm_list) n++;
     if (n == 0) return 93;
 
+
+    /* check io-predicates */
+
+    /* angular in on fwd, linear out */
+    P = proj_create (PJ_DEFAULT_CTX, "+proj=cart +ellps=GRS80");
+    if (0==P) return 0;
+    if (!proj_angular_input (P, PJ_FWD))  return 100;
+    if ( proj_angular_input (P, PJ_INV))  return 101;
+    if ( proj_angular_output (P, PJ_FWD)) return 102;
+    if (!proj_angular_output (P, PJ_INV)) return 103;
+    P->inverted = 1;
+    if ( proj_angular_input (P, PJ_FWD))  return 104;
+    if (!proj_angular_input (P, PJ_INV))  return 105;
+    if (!proj_angular_output (P, PJ_FWD)) return 106;
+    if ( proj_angular_output (P, PJ_INV)) return 107;
+    proj_destroy(P);
+
+    /* angular in and out */
+    P = proj_create(PJ_DEFAULT_CTX,
+        "+proj=molodensky +a=6378160 +rf=298.25 "
+        "+da=-23 +df=-8.120449e-8 +dx=-134 +dy=-48 +dz=149 "
+        "+abridged "
+    );
+    if (0==P) return 0;
+    if (!proj_angular_input (P, PJ_FWD))  return 108;
+    if (!proj_angular_input (P, PJ_INV))  return 109;
+    if (!proj_angular_output (P, PJ_FWD)) return 110;
+    if (!proj_angular_output (P, PJ_INV)) return 111;
+    P->inverted = 1;
+    if (!proj_angular_input (P, PJ_FWD))  return 112;
+    if (!proj_angular_input (P, PJ_INV))  return 113;
+    if (!proj_angular_output (P, PJ_FWD)) return 114;
+    if (!proj_angular_output (P, PJ_INV)) return 115;
+    proj_destroy(P);
+
+    /* linear in and out */
+    P = proj_create(PJ_DEFAULT_CTX,
+        " +proj=helmert +ellps=GRS80"
+        " +x=0.0127 +y=0.0065 +z=-0.0209 +s=0.00195"
+        " +rx=-0.00039 +ry=0.00080 +rz=-0.00114"
+        " +dx=-0.0029 +dy=-0.0002 +dz=-0.0006 +ds=0.00001"
+        " +drx=-0.00011 +dry=-0.00019 +drz=0.00007"
+        " +epoch=1988.0 +transpose"
+    );
+    if (0==P) return 0;
+    if (proj_angular_input (P, PJ_FWD))  return 116;
+    if (proj_angular_input (P, PJ_INV))  return 117;
+    if (proj_angular_output (P, PJ_FWD)) return 118;
+    if (proj_angular_output (P, PJ_INV)) return 119;
+    P->inverted = 1;
+    if (proj_angular_input (P, PJ_FWD))  return 120;
+    if (proj_angular_input (P, PJ_INV))  return 121;
+    if (proj_angular_output (P, PJ_FWD)) return 122;
+    if (proj_angular_output (P, PJ_INV)) return 123;
+    proj_destroy(P);
+
+
     return 0;
 }
 
diff --git a/src/cct.c b/src/cct.c
index de97f728..1bd25f1d 100644
--- a/src/cct.c
+++ b/src/cct.c
@@ -151,7 +151,7 @@ int main(int argc, char **argv) {
     OPTARGS *o;
     FILE *fout = stdout;
     char *buf;
-    int input_unit, output_unit, nfields = 4, direction = 1, verbose;
+    int nfields = 4, direction = 1, verbose;
     double fixed_z = HUGE_VAL, fixed_time = HUGE_VAL;
     int columns_xyzt[] = {1, 2, 3, 4};
     const char *longflags[]  = {"v=verbose", "h=help", "I=inverse", 0};
@@ -218,9 +218,6 @@ int main(int argc, char **argv) {
         P->inverted = !(P->inverted);
     direction = 1;
 
-    input_unit   =  proj_angular_left (P)?  PJ_IO_UNITS_RADIANS: PJ_IO_UNITS_METERS;
-    output_unit  =  proj_angular_right (P)? PJ_IO_UNITS_RADIANS: PJ_IO_UNITS_METERS;
-
     /* Allocate input buffer */
     buf = calloc (1, 10000);
     if (0==buf) {
@@ -258,12 +255,12 @@ int main(int argc, char **argv) {
             continue;
         }
 
-        if (PJ_IO_UNITS_RADIANS==input_unit) {
+        if (proj_angular_input (P, direction)) {
             point.lpzt.lam = proj_torad (point.lpzt.lam);
             point.lpzt.phi = proj_torad (point.lpzt.phi);
         }
         point = proj_trans (P, direction, point);
-        if (PJ_IO_UNITS_RADIANS==output_unit) {
+        if (proj_angular_output (P, direction)) {
             point.lpzt.lam = proj_todeg (point.lpzt.lam);
             point.lpzt.phi = proj_todeg (point.lpzt.phi);
         }
diff --git a/src/gie.c b/src/gie.c
index 2bca2f66..0a344043 100644
--- a/src/gie.c
+++ b/src/gie.c
@@ -421,7 +421,13 @@ static void finish_previous_operation (char *args) {
     (void) args;
 }
 
+/*****************************************************************************/
 static int operation (char *args) {
+/*****************************************************************************
+    Define the operation to apply to the input data (in ISO 19100 lingo,
+    an operation is the general term describing something that can be
+    either a conversion or a transformation)
+******************************************************************************/
     T.op_id++;
     strcpy (&(T.operation[0]), args);
     if (T.verbosity > 1) {
@@ -445,38 +451,44 @@ static int operation (char *args) {
 }
 
 
-static PJ_COORD torad_if_needed (PJ *P, PJ_DIRECTION dir, PJ_COORD a) {
-    enum pj_io_units u = P->left;
-    PJ_COORD c;
-    if (dir==PJ_INV)
-        u = P->right;
-    if (u==PJ_IO_UNITS_CLASSIC || u==PJ_IO_UNITS_METERS)
-        return a;
-
-    c.lpz.lam = proj_torad (T.a.lpz.lam);
-    c.lpz.phi = proj_torad (T.a.lpz.phi);
-    c.v[2] = T.a.v[2];
-    c.v[3] = T.a.v[3];
-
+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 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 int accept (char *args) {
-    int n, i;
+/* try to parse args as a PJ_COORD */
+static PJ_COORD parse_coord (char *args) {
+    int i;
     char *endp, *prev = args;
-    T.a = proj_coord (0,0,0,0);
-    n = 4;
+    PJ_COORD a = proj_coord (0,0,0,0);
+
     for (i = 0; i < 4; i++) {
-        T.a.v[i] = proj_strtod (prev, &endp);
-        if (prev==endp) {
-            n--;
-            T.a.v[i] = 0;
-            break;
-        }
+        double d = proj_strtod (prev, &endp);
+        if (prev==endp)
+            return i > 1? a: proj_coord_error ();
+        a.v[i] = d;
         prev = endp;
     }
-    T.a = torad_if_needed (T.P, T.dir, T.a);
+
+    return a;
+}
+
+
+/*****************************************************************************/
+static int accept (char *args) {
+/*****************************************************************************
+    Read ("ACCEPT") a 2, 3, or 4 dimensional input coordinate.
+******************************************************************************/
+    T.a = parse_coord (args);
     if (T.verbosity > 3)
         printf ("#  %s", args);
     return 0;
@@ -484,7 +496,12 @@ static int accept (char *args) {
 
 
 
+/*****************************************************************************/
 static int roundtrip (char *args) {
+/*****************************************************************************
+    Check how far we go from the ACCEPTed point when doing successive
+    back/forward transformation pairs.
+******************************************************************************/
     int ntrips;
     double d, r, ans;
     char *endp;
@@ -512,84 +529,99 @@ static int roundtrip (char *args) {
     return 0;
 }
 
-static int expect (char *args) {
-    double d;
-    enum pj_io_units unit;
-    char *endp, *prev = args;
-    int i;
-
-    T.e     =  proj_coord (0,0,0,0);
-    T.b     =  proj_coord (0,0,0,0);
-    T.nargs = 4;
-    for (i = 0; i < 4; i++) {
-        T.e.v[i] = proj_strtod (prev, &endp);
-        if (prev==endp) {
-            T.nargs--;
-            T.e.v[i] = 0;
-            break;
-        }
-        prev = endp;
-    }
-    T.e = torad_if_needed (T.P, T.dir==PJ_FWD? PJ_INV:PJ_FWD, T.e);
+static int expect_message (double d, char *args) {
+    T.op_ko++;
+    T.total_ko++;
 
-    T.b = proj_trans (T.P, T.dir, T.a);
-    T.b = torad_if_needed (T.P, T.dir==PJ_FWD? PJ_INV:PJ_FWD, T.b);
+    if (T.verbosity < 0)
+        return 1;
+    if (d > 1e6)
+        d = 999999.999999;
+    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, "     expected: %s\n", args);
+    fprintf (T.fout, "     got:      %.9f   %.9f", 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);
+    return 1;
+}
 
-    if (T.nargs < 2) {
-        T.op_ko++;
-        T.total_ko++;
-        if (T.verbosity > -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     Too few args: %s\n", opt_strip_path (T.curr_file), (int) lineno, args);
-        }
-     return 1;
+static int expect_message_cannot_parse (char *args) {
+    T.op_ko++;
+    T.total_ko++;
+    if (T.verbosity > -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     Too few args: %s\n", opt_strip_path (T.curr_file), (int) lineno, args);
     }
+    return 1;
+}
 
-    unit = T.dir==PJ_FWD? T.P->right: T.P->left;
-    if (PJ_IO_UNITS_CLASSIC==unit)
-        unit = PJ_IO_UNITS_METERS;
 
-    if (unit==PJ_IO_UNITS_RADIANS)
-        d = proj_lp_dist (T.P, T.b.lp, T.e.lp);
+/*****************************************************************************/
+static int expect (char *args) {
+/*****************************************************************************
+    Tell GIE what to expect, when transforming the ACCEPTed input
+******************************************************************************/
+    PJ_COORD ci, co, ce;
+    double d;
+
+    T.e  =  parse_coord (args);
+    if (HUGE_VAL==T.e.v[0])
+        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;
+
+    /* input ("accepted") values also probably in degrees */
+    ci = proj_angular_input  (T.P, T.dir)? torad_coord (T.a): T.a;
+
+    /* angular output from proj_trans comes in radians */
+    co = proj_trans (T.P, T.dir, ci);
+    T.b = proj_angular_output (T.P, T.dir)? todeg_coord (co): co;
+
+    /* but there are a few more possible input conventions... */
+    if (proj_angular_output (T.P, T.dir)) {
+        double e = HUGE_VAL;
+        d = hypot (proj_lp_dist (T.P, ce.lp, co.lp), ce.lpz.z - co.lpz.z);
+        /* check whether input was already in radians */
+        if (d > T.tolerance)
+            e = hypot (proj_lp_dist (T.P, T.e.lp, co.lp), T.e.lpz.z - co.lpz.z);
+        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;
+    }
     else
         d = proj_xyz_dist (T.b.xyz, T.e.xyz);
 
-    if (d > T.tolerance) {
-        if (d > 1e6)
-            d = 999999.999999;
-        if (T.verbosity > -1) {
-            if (0==T.op_ko && T.verbosity < 2)
-                banner (T.operation);
-            fprintf (T.fout, "%s", T.op_ko? "     -----\n": delim);
+    if (d > T.tolerance)
+        return expect_message (d, args);
+
+    T.op_ok++;
+    T.total_ok++;
 
-            fprintf (T.fout, "     FAILURE in %s(%d):\n", opt_strip_path (T.curr_file), (int) lineno);
-            fprintf (T.fout, "     expected: %s\n", args);
-            fprintf (T.fout, "     got:      %.9f   %.9f", T.b.xy.x,  T.b.xy.y);
-            if (T.nargs > 2)
-                fprintf (T.fout, "   %.9f", T.b.xyz.z);
-            if (T.nargs > 3)
-                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);
-        }
-        T.op_ko++;
-        T.total_ko++;
-    }
-    else {
-        T.op_ok++;
-        T.total_ok++;
-    }
     return 0;
 }
 
 
 
-
-
-
+/*****************************************************************************/
 static int verbose (char *args) {
+/*****************************************************************************
+    Tell the system how noisy it should be
+******************************************************************************/
     int i = (int) proj_atof (args);
 
     /* if -q/--quiet flag has been given, we do nothing */
@@ -603,14 +635,22 @@ static int verbose (char *args) {
     return 0;
 }
 
+/*****************************************************************************/
 static int comment (char *args) {
+/*****************************************************************************
+    in line comment. Equivalent to #
+******************************************************************************/
     (void) args;
     return 0;
 }
 
 
+/*****************************************************************************/
 static int echo (char *args) {
-    fprintf (T.fout, "%s\n", args);
+/*****************************************************************************
+    Add user defined noise to the output stream
+******************************************************************************/
+fprintf (T.fout, "%s\n", args);
     return 0;
 }
 
diff --git a/src/proj.def b/src/proj.def
index 998a260e..8b1882ee 100644
--- a/src/proj.def
+++ b/src/proj.def
@@ -144,5 +144,5 @@ EXPORTS
     proj_list_units                @130
     proj_list_prime_meridians      @131
 
-    proj_angular_left              @132
-    proj_angular_right             @133
+    proj_angular_input             @132
+    proj_angular_output            @133
diff --git a/src/proj.h b/src/proj.h
index 7607c3d8..bc0e1e06 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -371,8 +371,9 @@ typedef enum PJ_DIRECTION PJ_DIRECTION;
 
 
 
-int proj_angular_left (PJ *P);
-int proj_angular_right (PJ *P);
+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);
 
diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 9ae4147f..07208698 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -44,14 +44,29 @@ PJ_COORD proj_coord (double x, double y, double z, double t) {
     return res;
 }
 
-/* We do not want to bubble enum pj_io_units up to the API level, so we provide these predicates instead */
-int proj_angular_left (PJ *P) {
-    return pj_left (P)==PJ_IO_UNITS_RADIANS;
-}
-int proj_angular_right (PJ *P) {
+/*****************************************************************************/
+int proj_angular_input (PJ *P, enum PJ_DIRECTION dir) {
+/******************************************************************************
+    Returns 1 if the operator P expects angular input coordinates when
+    operating in direction dir, 0 otherwise.
+    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;
 }
 
+/*****************************************************************************/
+int proj_angular_output (PJ *P, enum PJ_DIRECTION dir) {
+/******************************************************************************
+    Returns 1 if the operator P provides angular output coordinates when
+    operating in direction dir, 0 otherwise.
+    dir: {PJ_FWD, PJ_INV}
+******************************************************************************/
+    return proj_angular_input (P, -dir);
+}
+
+
 /* Geodesic distance (in meter) between two points with angular 2D coordinates */
 double proj_lp_dist (const PJ *P, LP a, LP b) {
     double s12, azi1, azi2;
@@ -76,7 +91,6 @@ double proj_xyz_dist (XYZ a, XYZ b) {
 double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) {
     int i;
     PJ_COORD o, u;
-    enum pj_io_units unit;
 
     if (0==P)
         return HUGE_VAL;
@@ -106,9 +120,8 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) {
             return HUGE_VAL;
     }
 
-    /* left when forward, because we do a roundtrip, and end where we begin */
-    unit = direction==PJ_FWD?  P->left:  P->right;
-    if (unit==PJ_IO_UNITS_RADIANS)
+    /* checking for angular input since we do a roundtrip, and end where we begin */
+    if (proj_angular_input (P, direction))
         return hypot (proj_lp_dist (P, coo.lp, o.lp), coo.lpz.z - o.lpz.z);
 
     return proj_xyz_dist (coo.xyz, coo.xyz);