Make the 4D API fully 4D (#788)

Make 4D API fully 4D:

Remove PJ_XY, PJ_LP, PJ_XYZ, PJ_LPZ etc. from the API surface and
make all formal parameters of the API fully 4D PJ_COORD.

This operation primarily influences the proj_XXX_dist functions,
which mostly work by calling Charles Karney's geodesic
subsystem, keeping the distance, and throwing away the start and
end azimuths for the geodesic computed.

Also a PJ_GEOD(esic) persona is introduced for the PJ_COORD type.
The proj_geod function returns a PJ_GEOD, representing the geodesic
between the points represented by its PJ_COORD arguments.

Finally, the proj_factors functions had its lp argument changed
from PJ_LP to PJ_COORD.

1 files changed, 31 insertions, 16 deletions

diff --git a/src/proj_4D_api.c b/src/proj_4D_api.c
index 8f6305a4..2f689085 100644
--- a/src/proj_4D_api.c
+++ b/src/proj_4D_api.c
@@ -67,32 +67,47 @@ int proj_angular_output (PJ *P, enum PJ_DIRECTION dir) {
 }
 
 
+/* Geodesic distance (in meter) + fwd and rev azimuth between two points on the ellipsoid */
+PJ_COORD proj_geod (const PJ *P, PJ_COORD a, PJ_COORD b) {
+    PJ_COORD c;
+    /* Note: the geodesic code takes arguments in degrees */
+    geod_inverse (P->geod,
+        PJ_TODEG(a.lpz.phi), PJ_TODEG(a.lpz.lam),
+        PJ_TODEG(b.lpz.phi), PJ_TODEG(b.lpz.lam),
+        c.v, c.v+1, c.v+2
+    );
+
+    return c;
+}
+
+
 /* Geodesic distance (in meter) between two points with angular 2D coordinates */
-double proj_lp_dist (const PJ *P, LP a, LP b) {
+double proj_lp_dist (const PJ *P, PJ_COORD a, PJ_COORD b) {
     double s12, azi1, azi2;
     /* Note: the geodesic code takes arguments in degrees */
-    geod_inverse (P->geod, PJ_TODEG(a.phi), PJ_TODEG(a.lam), PJ_TODEG(b.phi), PJ_TODEG(b.lam), &s12, &azi1, &azi2);
+    geod_inverse (P->geod,
+        PJ_TODEG(a.lpz.phi), PJ_TODEG(a.lpz.lam),
+        PJ_TODEG(b.lpz.phi), PJ_TODEG(b.lpz.lam),
+        &s12, &azi1, &azi2
+    );
     return s12;
 }
 
 /* The geodesic distance AND the vertical offset */
-double proj_lpz_dist (const PJ *P, LPZ a, LPZ b) {
-    PJ_COORD aa, bb;
-    aa.lpz = a;
-    bb.lpz = b;
-    if (HUGE_VAL==a.lam || HUGE_VAL==b.lam)
+double proj_lpz_dist (const PJ *P, PJ_COORD a, PJ_COORD b) {
+    if (HUGE_VAL==a.lpz.lam || HUGE_VAL==b.lpz.lam)
         return HUGE_VAL;
-    return hypot (proj_lp_dist (P, aa.lp, bb.lp), a.z - b.z);
+    return hypot (proj_lp_dist (P, a, b), a.lpz.z - b.lpz.z);
 }
 
 /* Euclidean distance between two points with linear 2D coordinates */
-double proj_xy_dist (XY a, XY b) {
-    return hypot (a.x - b.x, a.y - b.y);
+double proj_xy_dist (PJ_COORD a, PJ_COORD b) {
+    return hypot (a.xy.x - b.xy.x, a.xy.y - b.xy.y);
 }
 
 /* Euclidean distance between two points with linear 3D coordinates */
-double proj_xyz_dist (XYZ a, XYZ b) {
-    return hypot (hypot (a.x - b.x, a.y - b.y), a.z - b.z);
+double proj_xyz_dist (PJ_COORD a, PJ_COORD b) {
+    return hypot (proj_xy_dist (a, b), a.xyz.z - b.xyz.z);
 }
 
 
@@ -125,9 +140,9 @@ double proj_roundtrip (PJ *P, PJ_DIRECTION direction, int n, PJ_COORD *coo) {
 
     /* 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, t.lpz);
+        return proj_lpz_dist (P, org, t);
 
-    return proj_xyz_dist (org.xyz, t.xyz);
+    return proj_xyz_dist (org, t);
 }
 
 
@@ -996,7 +1011,7 @@ PJ_INIT_INFO proj_init_info(const char *initname){
 
 
 /*****************************************************************************/
-PJ_FACTORS proj_factors(PJ *P, LP lp) {
+PJ_FACTORS proj_factors(PJ *P, PJ_COORD lp) {
 /******************************************************************************
     Cartographic characteristics at point lp.
 
@@ -1012,7 +1027,7 @@ PJ_FACTORS proj_factors(PJ *P, LP lp) {
     if (0==P)
         return factors;
 
-    if (pj_factors(lp, P, 0.0, &f))
+    if (pj_factors(lp.lp, P, 0.0, &f))
         return factors;
 
     factors.meridional_scale  =  f.h;