Remove old geocent.h/.cpp code

Last user, apart from transform.cpp, was the SCH projection. Modify it
to use cart instead.

And move content of geocent.h and .cpp into transform.cpp directly, so
that it can be later wiped up easily.

1 files changed, 70 insertions, 69 deletions

diff --git a/src/projections/sch.cpp b/src/projections/sch.cpp
index 261555a7..6dcb7cac 100644
--- a/src/projections/sch.cpp
+++ b/src/projections/sch.cpp
@@ -37,7 +37,6 @@
 
 #include "proj.h"
 #include "proj_internal.h"
-#include "geocent.h"
 
 namespace { // anonymous namespace
 struct pj_opaque {
@@ -48,33 +47,29 @@ struct pj_opaque {
     double transMat[9];
     double xyzoff[3];
     double rcurv;
-    GeocentricInfo sph;
-    GeocentricInfo elp_0;
+    PJ* cart;
+    PJ* cart_sph;
 };
 } // anonymous namespace
 
 PROJ_HEAD(sch, "Spherical Cross-track Height") "\n\tMisc\n\tplat_0= plon_0= phdg_0= [h_0=]";
 
 static PJ_LPZ sch_inverse3d(PJ_XYZ xyz, PJ *P) {
-    PJ_LPZ lpz = {0.0, 0.0, 0.0};
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double temp[3];
-    /* Local lat,lon using radius */
-    double pxyz[] = {
-        xyz.y * P->a / Q->rcurv,
-        xyz.x * P->a / Q->rcurv,
-        xyz.z
-    };
 
-    if( pj_Convert_Geodetic_To_Geocentric( &(Q->sph), pxyz[0], pxyz[1], pxyz[2], temp, temp+1, temp+2) != 0) {
-            proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-            return lpz;
-    }
+    PJ_LPZ lpz;
+    lpz.lam = xyz.x * (P->a / Q->rcurv);
+    lpz.phi = xyz.y * (P->a / Q->rcurv);
+    lpz.z = xyz.z;
+    xyz  =  Q->cart_sph->fwd3d (lpz, Q->cart_sph);
+    double temp[] = {xyz.x, xyz.y, xyz.z};
 
     /* Apply rotation */
-    pxyz[0] = Q->transMat[0] * temp[0] + Q->transMat[1] * temp[1] + Q->transMat[2] * temp[2];
-    pxyz[1] = Q->transMat[3] * temp[0] + Q->transMat[4] * temp[1] + Q->transMat[5] * temp[2];
-    pxyz[2] = Q->transMat[6] * temp[0] + Q->transMat[7] * temp[1] + Q->transMat[8] * temp[2];
+    double pxyz[] = {
+        Q->transMat[0] * temp[0] + Q->transMat[1] * temp[1] + Q->transMat[2] * temp[2],
+        Q->transMat[3] * temp[0] + Q->transMat[4] * temp[1] + Q->transMat[5] * temp[2],
+        Q->transMat[6] * temp[0] + Q->transMat[7] * temp[1] + Q->transMat[8] * temp[2]
+    };
 
     /* Apply offset */
     pxyz[0] += Q->xyzoff[0];
@@ -82,28 +77,19 @@ static PJ_LPZ sch_inverse3d(PJ_XYZ xyz, PJ *P) {
     pxyz[2] += Q->xyzoff[2];
 
     /* Convert geocentric coordinates to lat lon */
-    pj_Convert_Geocentric_To_Geodetic( &(Q->elp_0), pxyz[0], pxyz[1], pxyz[2],
-            temp, temp+1, temp+2);
-
-
-    lpz.lam = temp[1] ;
-    lpz.phi = temp[0] ;
-    lpz.z = temp[2];
-
+    xyz.x = pxyz[0];
+    xyz.y = pxyz[1];
+    xyz.z = pxyz[2];
+    lpz  =  Q->cart->inv3d (xyz, Q->cart);
     return lpz;
 }
 
 static PJ_XYZ sch_forward3d(PJ_LPZ lpz, PJ *P) {
-    PJ_XYZ xyz = {0.0, 0.0, 0.0};
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double temp[3];
 
     /* Convert lat lon to geocentric coordinates */
-    if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), lpz.phi, lpz.lam, lpz.z, temp, temp+1, temp+2 ) != 0 ) {
-        proj_errno_set(P, PJD_ERR_TOLERANCE_CONDITION);
-        return xyz;
-    }
-
+    PJ_XYZ xyz  =  Q->cart->fwd3d (lpz, Q->cart);
+    double temp[] = {xyz.x, xyz.y, xyz.z};
 
     /* Adjust for offset */
     temp[0] -= Q->xyzoff[0];
@@ -119,52 +105,67 @@ static PJ_XYZ sch_forward3d(PJ_LPZ lpz, PJ *P) {
     };
 
     /* Convert to local lat,lon */
-    pj_Convert_Geocentric_To_Geodetic( &(Q->sph), pxyz[0], pxyz[1], pxyz[2],
-            temp, temp+1, temp+2);
-
+    xyz.x = pxyz[0];
+    xyz.y = pxyz[1];
+    xyz.z = pxyz[2];
+    lpz  =  Q->cart_sph->inv3d (xyz, Q->cart_sph);
 
     /* Scale by radius */
-    xyz.x = temp[1] * Q->rcurv / P->a;
-    xyz.y = temp[0] * Q->rcurv / P->a;
-    xyz.z = temp[2];
+    xyz.x = lpz.lam * (Q->rcurv / P->a);
+    xyz.y = lpz.phi * (Q->rcurv / P->a);
+    xyz.z = lpz.z;
 
     return xyz;
 }
 
+static PJ *destructor (PJ *P, int errlev) {
+    if (nullptr==P)
+        return nullptr;
+
+    auto Q = static_cast<struct pj_opaque*>(P->opaque);
+    if( Q )
+    {
+        if (Q->cart)
+            Q->cart->destructor (Q->cart, errlev);
+        if (Q->cart_sph)
+            Q->cart_sph->destructor (Q->cart_sph, errlev);
+    }
+
+    return pj_default_destructor(P, errlev);
+}
 
 static PJ *setup(PJ *P) { /* general initialization */
     struct pj_opaque *Q = static_cast<struct pj_opaque*>(P->opaque);
-    double reast, rnorth;
-    double chdg, shdg;
-    double clt, slt;
-    double clo, slo;
-    double temp;
-    double pxyz[3];
-
-    temp = P->a * sqrt(1.0 - P->es);
 
     /* Setup original geocentric system */
-    if ( pj_Set_Geocentric_Parameters(&(Q->elp_0), P->a, temp) != 0)
-        return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+    // Pass a dummy ellipsoid definition that will be overridden just afterwards
+    Q->cart = proj_create(P->ctx, "+proj=cart +a=1");
+    if (Q->cart == nullptr)
+        return destructor(P, ENOMEM);
+
+    /* inherit ellipsoid definition from P to Q->cart */
+    pj_inherit_ellipsoid_def (P, Q->cart);
 
-    clt = cos(Q->plat);
-    slt = sin(Q->plat);
-    clo = cos(Q->plon);
-    slo = sin(Q->plon);
+    const double clt = cos(Q->plat);
+    const double slt = sin(Q->plat);
+    const double clo = cos(Q->plon);
+    const double slo = sin(Q->plon);
 
     /* Estimate the radius of curvature for given peg */
-    temp = sqrt(1.0 - (P->es) * slt * slt);
-    reast = (P->a)/temp;
-    rnorth = (P->a) * (1.0 - (P->es))/pow(temp,3);
+    const double temp = sqrt(1.0 - (P->es) * slt * slt);
+    const double reast = (P->a)/temp;
+    const double rnorth = (P->a) * (1.0 - (P->es))/pow(temp,3);
 
-    chdg = cos(Q->phdg);
-    shdg = sin(Q->phdg);
+    const double chdg = cos(Q->phdg);
+    const double shdg = sin(Q->phdg);
 
     Q->rcurv = Q->h0 + (reast*rnorth)/(reast * chdg * chdg + rnorth * shdg * shdg);
 
     /* Set up local sphere at the given peg point */
-    if ( pj_Set_Geocentric_Parameters(&(Q->sph), Q->rcurv, Q->rcurv) != 0)
-            return pj_default_destructor(P, PJD_ERR_FAILED_TO_FIND_PROJ);
+    Q->cart_sph = proj_create(P->ctx, "+proj=cart +a=1");
+    if (Q->cart_sph == nullptr)
+        return destructor(P, ENOMEM);
+    pj_calc_ellipsoid_params(Q->cart_sph, Q->rcurv, 0);
 
     /* Set up the transformation matrices */
     Q->transMat[0] = clt * clo;
@@ -177,15 +178,14 @@ static PJ *setup(PJ *P) { /* general initialization */
     Q->transMat[7] =  clt*chdg;
     Q->transMat[8] =  clt*shdg;
 
-
-    if( pj_Convert_Geodetic_To_Geocentric( &(Q->elp_0), Q->plat, Q->plon, Q->h0,
-                                           pxyz, pxyz+1, pxyz+2 ) != 0 )
-        return pj_default_destructor(P, PJD_ERR_LAT_OR_LON_EXCEED_LIMIT);
-
-
-    Q->xyzoff[0] = pxyz[0] - (Q->rcurv) * clt * clo;
-    Q->xyzoff[1] = pxyz[1] - (Q->rcurv) * clt * slo;
-    Q->xyzoff[2] = pxyz[2] - (Q->rcurv) * slt;
+    PJ_LPZ lpz;
+    lpz.lam = Q->plon;
+    lpz.phi = Q->plat;
+    lpz.z = Q->h0;
+    PJ_XYZ xyz  =  Q->cart->fwd3d (lpz, Q->cart);
+    Q->xyzoff[0] = xyz.x - (Q->rcurv) * clt * clo;
+    Q->xyzoff[1] = xyz.y- (Q->rcurv) * clt * slo;
+    Q->xyzoff[2] = xyz.z - (Q->rcurv) * slt;
 
     P->fwd3d = sch_forward3d;
     P->inv3d = sch_inverse3d;
@@ -198,6 +198,7 @@ PJ *PROJECTION(sch) {
     if (nullptr==Q)
         return pj_default_destructor(P, ENOMEM);
     P->opaque = Q;
+    P->destructor = destructor;
 
     Q->h0 = 0.0;