Switch proj_roundtrip to accept PJ_COORD, rather than PJ_OBS, and make it do proper geodesic distances for forward roundtrips

7 files changed, 50 insertions, 31 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_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.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);