Merge pull request #592 from busstoptaktik/obs_api-improvements

Cleaned up the use of PJ_IO_UNITS_CLASSIC, so it means only one thing: That output is linear and measured in units of the semimajor axis. This made it possible to remove a number of checks that had now become superfluous

Corrected pj_roundtrip, so it takes PJ_COORD args, rather than PJ_OBS, and made it measure geodesic distances, rather than cartesian, where it makes sense

13 files changed, 59 insertions, 38 deletions

diff --git a/src/PJ_cart.c b/src/PJ_cart.c
index eb2a66bc..d34996eb 100644
--- a/src/PJ_cart.c
+++ b/src/PJ_cart.c
@@ -315,8 +315,8 @@ int pj_cart_selftest (void) {
     b = proj_trans_obs (P, PJ_FWD, a);
 
     /* Check roundtrip precision for 10000 iterations each way */
-    dist = proj_roundtrip (P, PJ_FWD, 10000, a);
-    dist = proj_roundtrip (P, PJ_INV, 10000, b);
+    dist = proj_roundtrip (P, PJ_FWD, 10000, a.coo);
+    dist = proj_roundtrip (P, PJ_INV, 10000, b.coo);
     if (dist > 2e-9)
         return 7;
 
@@ -328,7 +328,7 @@ int pj_cart_selftest (void) {
     a.coo.lpz.z   = 100;
 
     /* Forward projection: Ellipsoidal-to-3D-Cartesian */
-    dist = proj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
     if (dist > 1e-12)
         return 8;
 
@@ -339,7 +339,7 @@ int pj_cart_selftest (void) {
     a.coo.lpz.z   = 100;
 
     /* Forward projection: Ellipsoidal-to-3D-Cartesian */
-    dist = proj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
     if (dist > 1e-12)
         return 9;
 
diff --git a/src/PJ_hgridshift.c b/src/PJ_hgridshift.c
index 1b65b1bd..0adc9e00 100644
--- a/src/PJ_hgridshift.c
+++ b/src/PJ_hgridshift.c
@@ -136,7 +136,7 @@ int pj_hgridshift_selftest (void) {
     a.coo.lpz.lam = PJ_TORAD(173);
     a.coo.lpz.phi = PJ_TORAD(-45);
     
-    dist = proj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
     if (dist > 0.00000001)
         return 1;
 
diff --git a/src/PJ_horner.c b/src/PJ_horner.c
index c29cdcd0..d6d2c51c 100644
--- a/src/PJ_horner.c
+++ b/src/PJ_horner.c
@@ -533,7 +533,7 @@ int pj_horner_selftest (void) {
     a.coo.uv.u =  878354.8539;
 
     /* Check roundtrip precision for 1 iteration each way, starting in forward direction */
-    dist = proj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
     if (dist > 0.01)
         return 1;
 
@@ -562,7 +562,7 @@ int pj_horner_selftest (void) {
         return 3;
 
     /* Check roundtrip precision for 1 iteration each way */
-    dist = proj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
     if (dist > 0.01)
         return 4;
 
diff --git a/src/PJ_molodensky.c b/src/PJ_molodensky.c
index fc68ee50..f09f07cd 100644
--- a/src/PJ_molodensky.c
+++ b/src/PJ_molodensky.c
@@ -341,7 +341,7 @@ int pj_molodensky_selftest (void) {
     }
 
     /* let's try a roundtrip */
-    if (proj_roundtrip(P, PJ_FWD, 100, in) > 1) {
+    if (proj_roundtrip(P, PJ_FWD, 100, in.coo) > 1) {
         proj_destroy(P);
         return 12;
     }
@@ -370,7 +370,7 @@ int pj_molodensky_selftest (void) {
     }
 
     /* let's try a roundtrip */
-    if (proj_roundtrip(P, PJ_FWD, 100, in) > 1) {
+    if (proj_roundtrip(P, PJ_FWD, 100, in.coo) > 1) {
         proj_destroy(P);
         return 22;
     }
diff --git a/src/PJ_vgridshift.c b/src/PJ_vgridshift.c
index d58972f4..ededd544 100644
--- a/src/PJ_vgridshift.c
+++ b/src/PJ_vgridshift.c
@@ -118,7 +118,7 @@ int pj_vgridshift_selftest (void) {
     a.coo.lpz.lam = PJ_TORAD(12.5);
     a.coo.lpz.phi = PJ_TORAD(55.5);
 
-    dist = proj_roundtrip (P, PJ_FWD, 1, a);
+    dist = proj_roundtrip (P, PJ_FWD, 1, a.coo);
     if (dist > 0.00000001)
         return 1;
 
diff --git a/src/pj_fwd.c b/src/pj_fwd.c
index 02f18d1c..e010f6ec 100644
--- a/src/pj_fwd.c
+++ b/src/pj_fwd.c
@@ -17,8 +17,8 @@ pj_fwd(LP lp, PJ *P) {
         return err;
     last_errno = proj_errno_reset (P);
 
-    /* Check input coordinates if angular */
-    if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+    /* Check validity of angular input coordinates */
+    if (P->left==PJ_IO_UNITS_RADIANS) {
 
         /* check for forward and latitude or longitude overange */
         t = fabs(lp.phi)-M_HALFPI;
diff --git a/src/pj_fwd3d.c b/src/pj_fwd3d.c
index be683745..d141178f 100644
--- a/src/pj_fwd3d.c
+++ b/src/pj_fwd3d.c
@@ -20,8 +20,8 @@ XYZ pj_fwd3d(LPZ lpz, PJ *P) {
 
     last_errno = proj_errno_reset(P);
 
-    /* Check input coordinates if angular */
-    if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+    /* Check validity of angular input coordinates */
+    if (P->left==PJ_IO_UNITS_RADIANS) {
 
         /* check for forward and latitude or longitude overange */
         t = fabs(lpz.phi)-M_HALFPI;
diff --git a/src/pj_inv.c b/src/pj_inv.c
index 55fc917f..68a5595b 100644
--- a/src/pj_inv.c
+++ b/src/pj_inv.c
@@ -42,7 +42,7 @@ LP pj_inv(XY xy, PJ *P) {
     if (P->ctx->last_errno)
         return err;
 
-    if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+    if (P->left==PJ_IO_UNITS_RADIANS) {
         /* reduce from del lp.lam */
         lp.lam += P->lam0;
 
diff --git a/src/pj_inv3d.c b/src/pj_inv3d.c
index a01cfa7e..53e39a76 100644
--- a/src/pj_inv3d.c
+++ b/src/pj_inv3d.c
@@ -42,7 +42,7 @@ LPZ pj_inv3d (XYZ xyz, PJ *P) {
     if (P->ctx->last_errno)
         return err;
 
-    if ((P->left==PJ_IO_UNITS_CLASSIC)||(P->left==PJ_IO_UNITS_RADIANS)) {
+    if (P->left==PJ_IO_UNITS_RADIANS) {
         /* reduce from del lp.lam */
         lpz.lam += P->lam0;
 
diff --git a/src/pj_obs_api.c b/src/pj_obs_api.c
index acdf0701..57b47341 100644
--- a/src/pj_obs_api.c
+++ b/src/pj_obs_api.c
@@ -90,9 +90,10 @@ 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_OBS obs) {
+double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo) {
     int i;
-    PJ_OBS o, u;
+    PJ_COORD o, u;
+    enum pj_io_units unit;
 
     if (0==P)
         return HUGE_VAL;
@@ -102,28 +103,46 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_OBS obs) {
         return HUGE_VAL;
     }
 
-    o.coo = obs.coo;
-
-    for (i = 0;  i < n;  i++) {
-        switch (direction) {
-            case PJ_FWD:
-                u  =  pj_fwdobs (o, P);
-                o  =  pj_invobs (u, P);
-                break;
-            case PJ_INV:
-                u  =  pj_invobs (o, P);
-                o  =  pj_fwdobs (u, P);
-                break;
-            default:
-                proj_errno_set (P, EINVAL);
-                return HUGE_VAL;
-        }
+    o = coo;
+
+    switch (direction) {
+        case PJ_FWD:
+            for (i = 0;  i < n;  i++) {
+                u  =  pj_fwdcoord (o, P);
+                o  =  pj_invcoord (u, P);
+            }
+            break;
+        case PJ_INV:
+            for (i = 0;  i < n;  i++) {
+                u  =  pj_invcoord (o, P);
+                o  =  pj_fwdcoord (u, P);
+            }
+            break;
+        default:
+            proj_errno_set (P, EINVAL);
+            return HUGE_VAL;
     }
 
-    return proj_xyz_dist (o.coo.xyz, obs.coo.xyz);
+    /* 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)
+        return hypot (proj_lp_dist (P, coo.lp, o.lp), coo.lpz.z - o.lpz.z);
+    
+    return proj_xyz_dist (coo.xyz, coo.xyz);
 }
 
 
+
+
+
+
+
+
+
+
+
+
+
 /* Apply the transformation P to the coordinate coo */
 PJ_OBS proj_trans_obs (PJ *P, PJ_DIRECTION direction, PJ_OBS obs) {
     if (0==P)
diff --git a/src/proj.c b/src/proj.c
index aa0f07c0..049593c9 100644
--- a/src/proj.c
+++ b/src/proj.c
@@ -145,7 +145,7 @@ 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 || Proj->left == PJ_IO_UNITS_CLASSIC) {
+                if (Proj->left == PJ_IO_UNITS_RADIANS) {
                     data.v *= RAD_TO_DEG;
                     data.u *= RAD_TO_DEG;
                 }
diff --git a/src/proj.h b/src/proj.h
index 600d5b30..167fb1eb 100644
--- a/src/proj.h
+++ b/src/proj.h
@@ -382,8 +382,8 @@ PJ_COORD proj_coord (double x, double y, double z, double t);
 PJ_OBS   proj_obs   (double x, double y, double z, double t, double o, double p, double k, int id, unsigned int flags);
 
 /* Measure internal consistency - in forward or inverse direction */
-double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_OBS obs);
-
+double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD coo);
+    
 /* Geodesic distance between two points with angular 2D coordinates */
 double proj_lp_dist (const PJ *P, LP a, LP b);
 
diff --git a/src/projects.h b/src/projects.h
index 8b20edd5..9391835d 100644
--- a/src/projects.h
+++ b/src/projects.h
@@ -589,6 +589,8 @@ C_NAMESPACE PJ *pj_##name (PJ *P) {                          \
         return 0;                                            \
     P->destructor = pj_default_destructor;                   \
     P->descr = des_##name;                                   \
+    P->left  = PJ_IO_UNITS_RADIANS;                          \
+    P->right = PJ_IO_UNITS_CLASSIC;                          \
     return P;                                                \
 }                                                            \
 PJ *pj_projection_specific_setup_##name (PJ *P)