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.

7 files changed, 721 insertions, 709 deletions

diff --git a/src/Makefile.am b/src/Makefile.am
index fe816af0..5f144332 100644
--- a/src/Makefile.am
+++ b/src/Makefile.am
@@ -208,7 +208,7 @@ libproj_la_SOURCES = \
 	fileapi.cpp \
 	\
 	datums.cpp datum_set.cpp transform.cpp \
-	geocent.cpp geocent.h utils.cpp \
+	utils.cpp \
 	mutex.cpp initcache.cpp geodesic.c \
 	strtod.cpp \
 	\
diff --git a/src/geocent.cpp b/src/geocent.cpp
deleted file mode 100644
index cbcc1df5..00000000
--- a/src/geocent.cpp
+++ /dev/null
@@ -1,456 +0,0 @@
-/***************************************************************************/
-/* RSC IDENTIFIER:  GEOCENTRIC
- *
- * ABSTRACT
- *
- *    This component provides conversions between Geodetic coordinates (latitude,
- *    longitude in radians and height in meters) and Geocentric coordinates
- *    (X, Y, Z) in meters.
- *
- * ERROR HANDLING
- *
- *    This component checks parameters for valid values.  If an invalid value
- *    is found, the error code is combined with the current error code using 
- *    the bitwise or.  This combining allows multiple error codes to be
- *    returned. The possible error codes are:
- *
- *      GEOCENT_NO_ERROR        : No errors occurred in function
- *      GEOCENT_LAT_ERROR       : Latitude out of valid range
- *                                 (-90 to 90 degrees)
- *      GEOCENT_LON_ERROR       : Longitude out of valid range
- *                                 (-180 to 360 degrees)
- *      GEOCENT_A_ERROR         : Semi-major axis lessthan or equal to zero
- *      GEOCENT_B_ERROR         : Semi-minor axis lessthan or equal to zero
- *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
- *
- *
- * REUSE NOTES
- *
- *    GEOCENTRIC is intended for reuse by any application that performs
- *    coordinate conversions between geodetic coordinates and geocentric
- *    coordinates.
- *    
- *
- * REFERENCES
- *    
- *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
- *    Ralph Toms, February 1996  UCRL-JC-123138.
- *    
- *    Further information on GEOCENTRIC can be found in the Reuse Manual.
- *
- *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
- *                                 Geospatial Information Division
- *                                 7701 Telegraph Road
- *                                 Alexandria, VA  22310-3864
- *
- * LICENSES
- *
- *    None apply to this component.
- *
- * RESTRICTIONS
- *
- *    GEOCENTRIC has no restrictions.
- *
- * ENVIRONMENT
- *
- *    GEOCENTRIC was tested and certified in the following environments:
- *
- *    1. Solaris 2.5 with GCC version 2.8.1
- *    2. Windows 95 with MS Visual C++ version 6
- *
- * MODIFICATIONS
- *
- *    Date              Description
- *    ----              -----------
- *    25-02-97          Original Code
- *
- */
-
-
-/***************************************************************************/
-/*
- *                               INCLUDES
- */
-#include <math.h>
-#include "geocent.h"
-/*
- *    math.h     - is needed for calls to sin, cos, tan and sqrt.
- *    geocent.h  - is needed for Error codes and prototype error checking.
- */
-
-
-/***************************************************************************/
-/*
- *                               DEFINES
- */
-#define PI         3.14159265358979323e0
-#define PI_OVER_2  (PI / 2.0e0)
-#define FALSE      0
-#define TRUE       1
-#define COS_67P5   0.38268343236508977  /* cosine of 67.5 degrees */
-#define AD_C       1.0026000            /* Toms region 1 constant */
-
-
-/***************************************************************************/
-/*
- *                              FUNCTIONS     
- */
-
-
-long pj_Set_Geocentric_Parameters (GeocentricInfo *gi, double a, double b) 
-
-{ /* BEGIN Set_Geocentric_Parameters */
-/*
- * The function Set_Geocentric_Parameters receives the ellipsoid parameters
- * as inputs and sets the corresponding state variables.
- *
- *    a  : Semi-major axis, in meters.          (input)
- *    b  : Semi-minor axis, in meters.          (input)
- */
-    long Error_Code = GEOCENT_NO_ERROR;
-
-    if (a <= 0.0)
-        Error_Code |= GEOCENT_A_ERROR;
-    if (b <= 0.0)
-        Error_Code |= GEOCENT_B_ERROR;
-    if (a < b)
-        Error_Code |= GEOCENT_A_LESS_B_ERROR;
-    if (!Error_Code)
-    {
-        gi->Geocent_a = a;
-        gi->Geocent_b = b;
-        gi->Geocent_a2 = a * a;
-        gi->Geocent_b2 = b * b;
-        gi->Geocent_e2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_a2;
-        gi->Geocent_ep2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_b2;
-    }
-    return (Error_Code);
-} /* END OF Set_Geocentric_Parameters */
-
-
-void pj_Get_Geocentric_Parameters (GeocentricInfo *gi,
-                                   double *a, 
-                                   double *b)
-{ /* BEGIN Get_Geocentric_Parameters */
-/*
- * The function Get_Geocentric_Parameters returns the ellipsoid parameters
- * to be used in geocentric coordinate conversions.
- *
- *    a  : Semi-major axis, in meters.          (output)
- *    b  : Semi-minor axis, in meters.          (output)
- */
-
-    *a = gi->Geocent_a;
-    *b = gi->Geocent_b;
-} /* END OF Get_Geocentric_Parameters */
-
-
-long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo *gi,
-                                        double Latitude,
-                                        double Longitude,
-                                        double Height,
-                                        double *X,
-                                        double *Y,
-                                        double *Z) 
-{ /* BEGIN Convert_Geodetic_To_Geocentric */
-/*
- * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
- * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
- * according to the current ellipsoid parameters.
- *
- *    Latitude  : Geodetic latitude in radians                     (input)
- *    Longitude : Geodetic longitude in radians                    (input)
- *    Height    : Geodetic height, in meters                       (input)
- *    X         : Calculated Geocentric X coordinate, in meters    (output)
- *    Y         : Calculated Geocentric Y coordinate, in meters    (output)
- *    Z         : Calculated Geocentric Z coordinate, in meters    (output)
- *
- */
-  long Error_Code = GEOCENT_NO_ERROR;
-  double Rn;            /*  Earth radius at location  */
-  double Sin_Lat;       /*  sin(Latitude)  */
-  double Sin2_Lat;      /*  Square of sin(Latitude)  */
-  double Cos_Lat;       /*  cos(Latitude)  */
-
-  /*
-  ** Don't blow up if Latitude is just a little out of the value
-  ** range as it may just be a rounding issue.  Also removed longitude
-  ** test, it should be wrapped by cos() and sin().  NFW for PROJ.4, Sep/2001.
-  */
-  if( Latitude < -PI_OVER_2 && Latitude > -1.001 * PI_OVER_2 )
-      Latitude = -PI_OVER_2;
-  else if( Latitude > PI_OVER_2 && Latitude < 1.001 * PI_OVER_2 )
-      Latitude = PI_OVER_2;
-  else if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
-  { /* Latitude out of range */
-    Error_Code |= GEOCENT_LAT_ERROR;
-  }
-
-  if (!Error_Code)
-  { /* no errors */
-    if (Longitude > PI)
-      Longitude -= (2*PI);
-    Sin_Lat = sin(Latitude);
-    Cos_Lat = cos(Latitude);
-    Sin2_Lat = Sin_Lat * Sin_Lat;
-    Rn = gi->Geocent_a / (sqrt(1.0e0 - gi->Geocent_e2 * Sin2_Lat));
-    *X = (Rn + Height) * Cos_Lat * cos(Longitude);
-    *Y = (Rn + Height) * Cos_Lat * sin(Longitude);
-    *Z = ((Rn * (1 - gi->Geocent_e2)) + Height) * Sin_Lat;
-
-  }
-  return (Error_Code);
-} /* END OF Convert_Geodetic_To_Geocentric */
-
-/*
- * The function Convert_Geocentric_To_Geodetic converts geocentric
- * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, 
- * and height), according to the current ellipsoid parameters.
- *
- *    X         : Geocentric X coordinate, in meters.         (input)
- *    Y         : Geocentric Y coordinate, in meters.         (input)
- *    Z         : Geocentric Z coordinate, in meters.         (input)
- *    Latitude  : Calculated latitude value in radians.       (output)
- *    Longitude : Calculated longitude value in radians.      (output)
- *    Height    : Calculated height value, in meters.         (output)
- */
-
-#define USE_ITERATIVE_METHOD
-
-void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
-                                        double X,
-                                        double Y, 
-                                        double Z,
-                                        double *Latitude,
-                                        double *Longitude,
-                                        double *Height)
-{ /* BEGIN Convert_Geocentric_To_Geodetic */
-#if !defined(USE_ITERATIVE_METHOD)
-/*
- * The method used here is derived from 'An Improved Algorithm for
- * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
- */
-
-/* Note: Variable names follow the notation used in Toms, Feb 1996 */
-
-    double W;        /* distance from Z axis */
-    double W2;       /* square of distance from Z axis */
-    double T0;       /* initial estimate of vertical component */
-    double T1;       /* corrected estimate of vertical component */
-    double S0;       /* initial estimate of horizontal component */
-    double S1;       /* corrected estimate of horizontal component */
-    double Sin_B0;   /* sin(B0), B0 is estimate of Bowring aux variable */
-    double Sin3_B0;  /* cube of sin(B0) */
-    double Cos_B0;   /* cos(B0) */
-    double Sin_p1;   /* sin(phi1), phi1 is estimated latitude */
-    double Cos_p1;   /* cos(phi1) */
-    double Rn;       /* Earth radius at location */
-    double Sum;      /* numerator of cos(phi1) */
-    int At_Pole;     /* indicates location is in polar region */
-
-    At_Pole = FALSE;
-    if (X != 0.0)
-    {
-        *Longitude = atan2(Y,X);
-    }
-    else
-    {
-        if (Y > 0)
-        {
-            *Longitude = PI_OVER_2;
-        }
-        else if (Y < 0)
-        {
-            *Longitude = -PI_OVER_2;
-        }
-        else
-        {
-            At_Pole = TRUE;
-            *Longitude = 0.0;
-            if (Z > 0.0)
-            {  /* north pole */
-                *Latitude = PI_OVER_2;
-            }
-            else if (Z < 0.0)
-            {  /* south pole */
-                *Latitude = -PI_OVER_2;
-            }
-            else
-            {  /* center of earth */
-                *Latitude = PI_OVER_2;
-                *Height = -Geocent_b;
-                return;
-            } 
-        }
-    }
-    W2 = X*X + Y*Y;
-    W = sqrt(W2);
-    T0 = Z * AD_C;
-    S0 = sqrt(T0 * T0 + W2);
-    Sin_B0 = T0 / S0;
-    Cos_B0 = W / S0;
-    Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
-    T1 = Z + gi->Geocent_b * gi->Geocent_ep2 * Sin3_B0;
-    Sum = W - gi->Geocent_a * gi->Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0;
-    S1 = sqrt(T1*T1 + Sum * Sum);
-    Sin_p1 = T1 / S1;
-    Cos_p1 = Sum / S1;
-    Rn = gi->Geocent_a / sqrt(1.0 - gi->Geocent_e2 * Sin_p1 * Sin_p1);
-    if (Cos_p1 >= COS_67P5)
-    {
-        *Height = W / Cos_p1 - Rn;
-    }
-    else if (Cos_p1 <= -COS_67P5)
-    {
-        *Height = W / -Cos_p1 - Rn;
-    }
-    else
-    {
-        *Height = Z / Sin_p1 + Rn * (gi->Geocent_e2 - 1.0);
-    }
-    if (At_Pole == FALSE)
-    {
-        *Latitude = atan(Sin_p1 / Cos_p1);
-    }
-#else /* defined(USE_ITERATIVE_METHOD) */
-/*
-* Reference...
-* ============
-* Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für
-* das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover
-* Nr. 137, p. 130-131.
-
-* Programmed by GGA- Leibniz-Institute of Applied Geophysics
-*               Stilleweg 2
-*               D-30655 Hannover
-*               Federal Republic of Germany
-*               Internet: www.gga-hannover.de
-*
-*               Hannover, March 1999, April 2004.
-*               see also: comments in statements
-* remarks:
-* Mathematically exact and because of symmetry of rotation-ellipsoid,
-* each point (X,Y,Z) has at least two solutions (Latitude1,Longitude1,Height1) and
-* (Latitude2,Longitude2,Height2). Is point=(0.,0.,Z) (P=0.), so you get even
-* four solutions,	every two symmetrical to the semi-minor axis.
-* Here Height1 and Height2 have at least a difference in order of
-* radius of curvature (e.g. (0,0,b)=> (90.,0.,0.) or (-90.,0.,-2b);
-* (a+100.)*(sqrt(2.)/2.,sqrt(2.)/2.,0.) => (0.,45.,100.) or
-* (0.,225.,-(2a+100.))).
-* The algorithm always computes (Latitude,Longitude) with smallest |Height|.
-* For normal computations, that means |Height|<10000.m, algorithm normally
-* converges after to 2-3 steps!!!
-* But if |Height| has the amount of length of ellipsoid's axis
-* (e.g. -6300000.m),	algorithm needs about 15 steps.
-*/
-
-/* local definitions and variables */
-/* end-criterium of loop, accuracy of sin(Latitude) */
-#define genau   1.E-12
-#define genau2  (genau*genau)
-#define maxiter 30
-
-    double P;        /* distance between semi-minor axis and location */
-    double RR;       /* distance between center and location */
-    double CT;       /* sin of geocentric latitude */
-    double ST;       /* cos of geocentric latitude */
-    double RX;
-    double RK;
-    double RN;       /* Earth radius at location */
-    double CPHI0;    /* cos of start or old geodetic latitude in iterations */
-    double SPHI0;    /* sin of start or old geodetic latitude in iterations */
-    double CPHI;     /* cos of searched geodetic latitude */
-    double SPHI;     /* sin of searched geodetic latitude */
-    double SDPHI;    /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
-    int iter;        /* # of continuous iteration, max. 30 is always enough (s.a.) */
-
-    P = sqrt(X*X+Y*Y);
-    RR = sqrt(X*X+Y*Y+Z*Z);
-
-/*	special cases for latitude and longitude */
-    if (P/gi->Geocent_a < genau) {
-
-/*  special case, if P=0. (X=0., Y=0.) */
-	*Longitude = 0.;
-
-/*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
- *  of ellipsoid (=center of mass), Latitude becomes PI/2 */
-        if (RR/gi->Geocent_a < genau) {
-            *Latitude = PI_OVER_2;
-            *Height   = -gi->Geocent_b;
-            return ;
-
-        }
-    }
-    else {
-/*  ellipsoidal (geodetic) longitude
- *  interval: -PI < Longitude <= +PI */
-        *Longitude=atan2(Y,X);
-    }
-
-/* --------------------------------------------------------------
- * Following iterative algorithm was developed by
- * "Institut für Erdmessung", University of Hannover, July 1988.
- * Internet: www.ife.uni-hannover.de
- * Iterative computation of CPHI,SPHI and Height.
- * Iteration of CPHI and SPHI to 10**-12 radian resp.
- * 2*10**-7 arcsec.
- * --------------------------------------------------------------
- */
-    CT = Z/RR;
-    ST = P/RR;
-    {
-        const double denominator = 1.0-gi->Geocent_e2*(2.0-gi->Geocent_e2)*ST*ST;
-        if( denominator == 0 )
-        {
-            *Latitude = HUGE_VAL;
-            *Longitude = HUGE_VAL;
-            *Height = HUGE_VAL;
-            return;
-        }
-        RX = 1.0/sqrt(denominator);
-    }
-    CPHI0 = ST*(1.0-gi->Geocent_e2)*RX;
-    SPHI0 = CT*RX;
-    iter = 0;
-
-/* loop to find sin(Latitude) resp. Latitude
- * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
-    do
-    {
-        iter++;
-        RN = gi->Geocent_a/sqrt(1.0-gi->Geocent_e2*SPHI0*SPHI0);
-
-/*  ellipsoidal (geodetic) height */
-        *Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi->Geocent_e2*SPHI0*SPHI0);
-
-        /* avoid zero division */
-        if (RN+*Height==0.0) {
-            *Latitude = 0.0;
-            return;
-        }
-        RK = gi->Geocent_e2*RN/(RN+*Height);
-        {
-            const double denominator = 1.0-RK*(2.0-RK)*ST*ST;
-            if( denominator == 0 )
-            {
-                *Latitude = HUGE_VAL;
-                *Longitude = HUGE_VAL;
-                *Height = HUGE_VAL;
-                return;
-            }
-            RX = 1.0/sqrt(denominator);
-        }
-        CPHI = ST*(1.0-RK)*RX;
-        SPHI = CT*RX;
-        SDPHI = SPHI*CPHI0-CPHI*SPHI0;
-        CPHI0 = CPHI;
-        SPHI0 = SPHI;
-    }
-    while (SDPHI*SDPHI > genau2 && iter < maxiter);
-
-/*	ellipsoidal (geodetic) latitude */
-    *Latitude=atan2(SPHI, fabs(CPHI));
-
-#endif /* defined(USE_ITERATIVE_METHOD) */
-} /* END OF Convert_Geocentric_To_Geodetic */
diff --git a/src/geocent.h b/src/geocent.h
deleted file mode 100644
index d6e90d23..00000000
--- a/src/geocent.h
+++ /dev/null
@@ -1,179 +0,0 @@
-#ifndef GEOCENT_H
-#define GEOCENT_H
-
-/***************************************************************************/
-/* RSC IDENTIFIER:  GEOCENTRIC
- *
- * ABSTRACT
- *
- *    This component provides conversions between Geodetic coordinates (latitude,
- *    longitude in radians and height in meters) and Geocentric coordinates
- *    (X, Y, Z) in meters.
- *
- * ERROR HANDLING
- *
- *    This component checks parameters for valid values.  If an invalid value
- *    is found, the error code is combined with the current error code using 
- *    the bitwise or.  This combining allows multiple error codes to be
- *    returned. The possible error codes are:
- *
- *      GEOCENT_NO_ERROR        : No errors occurred in function
- *      GEOCENT_LAT_ERROR       : Latitude out of valid range
- *                                 (-90 to 90 degrees)
- *      GEOCENT_LON_ERROR       : Longitude out of valid range
- *                                 (-180 to 360 degrees)
- *      GEOCENT_A_ERROR         : Semi-major axis less than or equal to zero
- *      GEOCENT_B_ERROR         : Semi-minor axis less than or equal to zero
- *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
- *
- *
- * REUSE NOTES
- *
- *    GEOCENTRIC is intended for reuse by any application that performs
- *    coordinate conversions between geodetic coordinates and geocentric
- *    coordinates.
- *    
- *
- * REFERENCES
- *    
- *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
- *    Ralph Toms, February 1996  UCRL-JC-123138.
- *    
- *    Further information on GEOCENTRIC can be found in the Reuse Manual.
- *
- *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
- *                                 Geospatial Information Division
- *                                 7701 Telegraph Road
- *                                 Alexandria, VA  22310-3864
- *
- * LICENSES
- *
- *    None apply to this component.
- *
- * RESTRICTIONS
- *
- *    GEOCENTRIC has no restrictions.
- *
- * ENVIRONMENT
- *
- *    GEOCENTRIC was tested and certified in the following environments:
- *
- *    1. Solaris 2.5 with GCC version 2.8.1
- *    2. Windows 95 with MS Visual C++ version 6
- *
- * MODIFICATIONS
- *
- *    Date              Description
- *    ----              -----------
- *
- *
- */
-
-
-/***************************************************************************/
-/*
- *                              DEFINES
- */
-#define GEOCENT_NO_ERROR        0x0000
-#define GEOCENT_LAT_ERROR       0x0001
-#define GEOCENT_LON_ERROR       0x0002
-#define GEOCENT_A_ERROR         0x0004
-#define GEOCENT_B_ERROR         0x0008
-#define GEOCENT_A_LESS_B_ERROR  0x0010
-
-
-/***************************************************************************/
-/*
- *                              FUNCTION PROTOTYPES
- */
-
-/* ensure proper linkage to c++ programs */
-#ifdef __cplusplus
-extern "C" {
-#endif
-
-typedef struct 
-{
-    double Geocent_a;        /* Semi-major axis of ellipsoid in meters */
-    double Geocent_b;        /* Semi-minor axis of ellipsoid           */
-    double Geocent_a2;       /* Square of semi-major axis */
-    double Geocent_b2;       /* Square of semi-minor axis */
-    double Geocent_e2;       /* Eccentricity squared  */
-    double Geocent_ep2;      /* 2nd eccentricity squared */
-} GeocentricInfo;
-
-void pj_Init_Geocentric( GeocentricInfo *gi );
-long pj_Set_Geocentric_Parameters( GeocentricInfo *gi, 
-                                   double a, 
-                                   double b);
-
-/*
- * The function Set_Geocentric_Parameters receives the ellipsoid parameters
- * as inputs and sets the corresponding state variables.
- *
- *    a  : Semi-major axis, in meters.          (input)
- *    b  : Semi-minor axis, in meters.          (input)
- */
-
-
-void pj_Get_Geocentric_Parameters ( GeocentricInfo *gi,
-                                    double *a, 
-                                    double *b);
-
-/*
- * The function Get_Geocentric_Parameters returns the ellipsoid parameters
- * to be used in geocentric coordinate conversions.
- *
- *    a  : Semi-major axis, in meters.          (output)
- *    b  : Semi-minor axis, in meters.          (output)
- */
-
-
-long pj_Convert_Geodetic_To_Geocentric ( GeocentricInfo *gi,
-                                         double Latitude,
-                                         double Longitude,
-                                         double Height,
-                                         double *X,
-                                         double *Y,
-                                         double *Z);
-/*
- * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
- * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
- * according to the current ellipsoid parameters.
- *
- *    Latitude  : Geodetic latitude in radians                     (input)
- *    Longitude : Geodetic longitude in radians                    (input)
- *    Height    : Geodetic height, in meters                       (input)
- *    X         : Calculated Geocentric X coordinate, in meters.   (output)
- *    Y         : Calculated Geocentric Y coordinate, in meters.   (output)
- *    Z         : Calculated Geocentric Z coordinate, in meters.   (output)
- *
- */
-
-
-void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
-                                        double X,
-                                        double Y, 
-                                        double Z,
-                                        double *Latitude,
-                                        double *Longitude,
-                                        double *Height);
-/*
- * The function Convert_Geocentric_To_Geodetic converts geocentric
- * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, 
- * and height), according to the current ellipsoid parameters.
- *
- *    X         : Geocentric X coordinate, in meters.         (input)
- *    Y         : Geocentric Y coordinate, in meters.         (input)
- *    Z         : Geocentric Z coordinate, in meters.         (input)
- *    Latitude  : Calculated latitude value in radians.       (output)
- *    Longitude : Calculated longitude value in radians.      (output)
- *    Height    : Calculated height value, in meters.         (output)
- */
-
-
-#ifdef __cplusplus
-}
-#endif
-
-#endif /* GEOCENT_H */
diff --git a/src/lib_proj.cmake b/src/lib_proj.cmake
index 7c93b205..96cef9f4 100644
--- a/src/lib_proj.cmake
+++ b/src/lib_proj.cmake
@@ -217,8 +217,6 @@ set(SRC_LIBPROJ_CORE
   fileapi.cpp
   fwd.cpp
   gauss.cpp
-  geocent.cpp
-  geocent.h
   geodesic.c
   init.cpp
   initcache.cpp
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;
 
diff --git a/src/transform.cpp b/src/transform.cpp
index ea3d9ae2..cff89232 100644
--- a/src/transform.cpp
+++ b/src/transform.cpp
@@ -33,11 +33,660 @@
 
 #include "proj.h"
 #include "proj_internal.h"
-#include "geocent.h"
 #include "grids.hpp"
 
 using namespace NS_PROJ;
 
+
+///////////////////////////////////////////////////////////////////////////////
+/// From older geocent.h
+///////////////////////////////////////////////////////////////////////////////
+
+
+/***************************************************************************/
+/* RSC IDENTIFIER:  GEOCENTRIC
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between Geodetic coordinates (latitude,
+ *    longitude in radians and height in meters) and Geocentric coordinates
+ *    (X, Y, Z) in meters.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found, the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *      GEOCENT_NO_ERROR        : No errors occurred in function
+ *      GEOCENT_LAT_ERROR       : Latitude out of valid range
+ *                                 (-90 to 90 degrees)
+ *      GEOCENT_LON_ERROR       : Longitude out of valid range
+ *                                 (-180 to 360 degrees)
+ *      GEOCENT_A_ERROR         : Semi-major axis less than or equal to zero
+ *      GEOCENT_B_ERROR         : Semi-minor axis less than or equal to zero
+ *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
+ *
+ *
+ * REUSE NOTES
+ *
+ *    GEOCENTRIC is intended for reuse by any application that performs
+ *    coordinate conversions between geodetic coordinates and geocentric
+ *    coordinates.
+ *    
+ *
+ * REFERENCES
+ *    
+ *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
+ *    Ralph Toms, February 1996  UCRL-JC-123138.
+ *    
+ *    Further information on GEOCENTRIC can be found in the Reuse Manual.
+ *
+ *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
+ *                                 Geospatial Information Division
+ *                                 7701 Telegraph Road
+ *                                 Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    GEOCENTRIC has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    GEOCENTRIC was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC version 2.8.1
+ *    2. Windows 95 with MS Visual C++ version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                              DEFINES
+ */
+#define GEOCENT_NO_ERROR        0x0000
+#define GEOCENT_LAT_ERROR       0x0001
+#define GEOCENT_LON_ERROR       0x0002
+#define GEOCENT_A_ERROR         0x0004
+#define GEOCENT_B_ERROR         0x0008
+#define GEOCENT_A_LESS_B_ERROR  0x0010
+
+
+/***************************************************************************/
+/*
+ *                              FUNCTION PROTOTYPES
+ */
+
+/* ensure proper linkage to c++ programs */
+#ifdef __cplusplus
+extern "C" {
+#endif
+
+typedef struct 
+{
+    double Geocent_a;        /* Semi-major axis of ellipsoid in meters */
+    double Geocent_b;        /* Semi-minor axis of ellipsoid           */
+    double Geocent_a2;       /* Square of semi-major axis */
+    double Geocent_b2;       /* Square of semi-minor axis */
+    double Geocent_e2;       /* Eccentricity squared  */
+    double Geocent_ep2;      /* 2nd eccentricity squared */
+} GeocentricInfo;
+
+void pj_Init_Geocentric( GeocentricInfo *gi );
+long pj_Set_Geocentric_Parameters( GeocentricInfo *gi, 
+                                   double a, 
+                                   double b);
+
+/*
+ * The function Set_Geocentric_Parameters receives the ellipsoid parameters
+ * as inputs and sets the corresponding state variables.
+ *
+ *    a  : Semi-major axis, in meters.          (input)
+ *    b  : Semi-minor axis, in meters.          (input)
+ */
+
+
+void pj_Get_Geocentric_Parameters ( GeocentricInfo *gi,
+                                    double *a, 
+                                    double *b);
+
+/*
+ * The function Get_Geocentric_Parameters returns the ellipsoid parameters
+ * to be used in geocentric coordinate conversions.
+ *
+ *    a  : Semi-major axis, in meters.          (output)
+ *    b  : Semi-minor axis, in meters.          (output)
+ */
+
+
+long pj_Convert_Geodetic_To_Geocentric ( GeocentricInfo *gi,
+                                         double Latitude,
+                                         double Longitude,
+                                         double Height,
+                                         double *X,
+                                         double *Y,
+                                         double *Z);
+/*
+ * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
+ * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
+ * according to the current ellipsoid parameters.
+ *
+ *    Latitude  : Geodetic latitude in radians                     (input)
+ *    Longitude : Geodetic longitude in radians                    (input)
+ *    Height    : Geodetic height, in meters                       (input)
+ *    X         : Calculated Geocentric X coordinate, in meters.   (output)
+ *    Y         : Calculated Geocentric Y coordinate, in meters.   (output)
+ *    Z         : Calculated Geocentric Z coordinate, in meters.   (output)
+ *
+ */
+
+
+void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
+                                        double X,
+                                        double Y, 
+                                        double Z,
+                                        double *Latitude,
+                                        double *Longitude,
+                                        double *Height);
+/*
+ * The function Convert_Geocentric_To_Geodetic converts geocentric
+ * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, 
+ * and height), according to the current ellipsoid parameters.
+ *
+ *    X         : Geocentric X coordinate, in meters.         (input)
+ *    Y         : Geocentric Y coordinate, in meters.         (input)
+ *    Z         : Geocentric Z coordinate, in meters.         (input)
+ *    Latitude  : Calculated latitude value in radians.       (output)
+ *    Longitude : Calculated longitude value in radians.      (output)
+ *    Height    : Calculated height value, in meters.         (output)
+ */
+
+
+#ifdef __cplusplus
+}
+#endif
+
+
+///////////////////////////////////////////////////////////////////////////////
+/// From older geocent.cpp
+///////////////////////////////////////////////////////////////////////////////
+
+
+/***************************************************************************/
+/* RSC IDENTIFIER:  GEOCENTRIC
+ *
+ * ABSTRACT
+ *
+ *    This component provides conversions between Geodetic coordinates (latitude,
+ *    longitude in radians and height in meters) and Geocentric coordinates
+ *    (X, Y, Z) in meters.
+ *
+ * ERROR HANDLING
+ *
+ *    This component checks parameters for valid values.  If an invalid value
+ *    is found, the error code is combined with the current error code using 
+ *    the bitwise or.  This combining allows multiple error codes to be
+ *    returned. The possible error codes are:
+ *
+ *      GEOCENT_NO_ERROR        : No errors occurred in function
+ *      GEOCENT_LAT_ERROR       : Latitude out of valid range
+ *                                 (-90 to 90 degrees)
+ *      GEOCENT_LON_ERROR       : Longitude out of valid range
+ *                                 (-180 to 360 degrees)
+ *      GEOCENT_A_ERROR         : Semi-major axis lessthan or equal to zero
+ *      GEOCENT_B_ERROR         : Semi-minor axis lessthan or equal to zero
+ *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
+ *
+ *
+ * REUSE NOTES
+ *
+ *    GEOCENTRIC is intended for reuse by any application that performs
+ *    coordinate conversions between geodetic coordinates and geocentric
+ *    coordinates.
+ *    
+ *
+ * REFERENCES
+ *    
+ *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
+ *    Ralph Toms, February 1996  UCRL-JC-123138.
+ *    
+ *    Further information on GEOCENTRIC can be found in the Reuse Manual.
+ *
+ *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
+ *                                 Geospatial Information Division
+ *                                 7701 Telegraph Road
+ *                                 Alexandria, VA  22310-3864
+ *
+ * LICENSES
+ *
+ *    None apply to this component.
+ *
+ * RESTRICTIONS
+ *
+ *    GEOCENTRIC has no restrictions.
+ *
+ * ENVIRONMENT
+ *
+ *    GEOCENTRIC was tested and certified in the following environments:
+ *
+ *    1. Solaris 2.5 with GCC version 2.8.1
+ *    2. Windows 95 with MS Visual C++ version 6
+ *
+ * MODIFICATIONS
+ *
+ *    Date              Description
+ *    ----              -----------
+ *    25-02-97          Original Code
+ *
+ */
+
+
+/***************************************************************************/
+/*
+ *                               INCLUDES
+ */
+#include <math.h>
+//#include "geocent.h"
+/*
+ *    math.h     - is needed for calls to sin, cos, tan and sqrt.
+ *    geocent.h  - is needed for Error codes and prototype error checking.
+ */
+
+
+/***************************************************************************/
+/*
+ *                               DEFINES
+ */
+#define PI         3.14159265358979323e0
+#define PI_OVER_2  (PI / 2.0e0)
+#define FALSE      0
+#define TRUE       1
+#define COS_67P5   0.38268343236508977  /* cosine of 67.5 degrees */
+#define AD_C       1.0026000            /* Toms region 1 constant */
+
+
+/***************************************************************************/
+/*
+ *                              FUNCTIONS     
+ */
+
+
+long pj_Set_Geocentric_Parameters (GeocentricInfo *gi, double a, double b) 
+
+{ /* BEGIN Set_Geocentric_Parameters */
+/*
+ * The function Set_Geocentric_Parameters receives the ellipsoid parameters
+ * as inputs and sets the corresponding state variables.
+ *
+ *    a  : Semi-major axis, in meters.          (input)
+ *    b  : Semi-minor axis, in meters.          (input)
+ */
+    long Error_Code = GEOCENT_NO_ERROR;
+
+    if (a <= 0.0)
+        Error_Code |= GEOCENT_A_ERROR;
+    if (b <= 0.0)
+        Error_Code |= GEOCENT_B_ERROR;
+    if (a < b)
+        Error_Code |= GEOCENT_A_LESS_B_ERROR;
+    if (!Error_Code)
+    {
+        gi->Geocent_a = a;
+        gi->Geocent_b = b;
+        gi->Geocent_a2 = a * a;
+        gi->Geocent_b2 = b * b;
+        gi->Geocent_e2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_a2;
+        gi->Geocent_ep2 = (gi->Geocent_a2 - gi->Geocent_b2) / gi->Geocent_b2;
+    }
+    return (Error_Code);
+} /* END OF Set_Geocentric_Parameters */
+
+
+void pj_Get_Geocentric_Parameters (GeocentricInfo *gi,
+                                   double *a, 
+                                   double *b)
+{ /* BEGIN Get_Geocentric_Parameters */
+/*
+ * The function Get_Geocentric_Parameters returns the ellipsoid parameters
+ * to be used in geocentric coordinate conversions.
+ *
+ *    a  : Semi-major axis, in meters.          (output)
+ *    b  : Semi-minor axis, in meters.          (output)
+ */
+
+    *a = gi->Geocent_a;
+    *b = gi->Geocent_b;
+} /* END OF Get_Geocentric_Parameters */
+
+
+long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo *gi,
+                                        double Latitude,
+                                        double Longitude,
+                                        double Height,
+                                        double *X,
+                                        double *Y,
+                                        double *Z) 
+{ /* BEGIN Convert_Geodetic_To_Geocentric */
+/*
+ * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
+ * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
+ * according to the current ellipsoid parameters.
+ *
+ *    Latitude  : Geodetic latitude in radians                     (input)
+ *    Longitude : Geodetic longitude in radians                    (input)
+ *    Height    : Geodetic height, in meters                       (input)
+ *    X         : Calculated Geocentric X coordinate, in meters    (output)
+ *    Y         : Calculated Geocentric Y coordinate, in meters    (output)
+ *    Z         : Calculated Geocentric Z coordinate, in meters    (output)
+ *
+ */
+  long Error_Code = GEOCENT_NO_ERROR;
+  double Rn;            /*  Earth radius at location  */
+  double Sin_Lat;       /*  sin(Latitude)  */
+  double Sin2_Lat;      /*  Square of sin(Latitude)  */
+  double Cos_Lat;       /*  cos(Latitude)  */
+
+  /*
+  ** Don't blow up if Latitude is just a little out of the value
+  ** range as it may just be a rounding issue.  Also removed longitude
+  ** test, it should be wrapped by cos() and sin().  NFW for PROJ.4, Sep/2001.
+  */
+  if( Latitude < -PI_OVER_2 && Latitude > -1.001 * PI_OVER_2 )
+      Latitude = -PI_OVER_2;
+  else if( Latitude > PI_OVER_2 && Latitude < 1.001 * PI_OVER_2 )
+      Latitude = PI_OVER_2;
+  else if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
+  { /* Latitude out of range */
+    Error_Code |= GEOCENT_LAT_ERROR;
+  }
+
+  if (!Error_Code)
+  { /* no errors */
+    if (Longitude > PI)
+      Longitude -= (2*PI);
+    Sin_Lat = sin(Latitude);
+    Cos_Lat = cos(Latitude);
+    Sin2_Lat = Sin_Lat * Sin_Lat;
+    Rn = gi->Geocent_a / (sqrt(1.0e0 - gi->Geocent_e2 * Sin2_Lat));
+    *X = (Rn + Height) * Cos_Lat * cos(Longitude);
+    *Y = (Rn + Height) * Cos_Lat * sin(Longitude);
+    *Z = ((Rn * (1 - gi->Geocent_e2)) + Height) * Sin_Lat;
+
+  }
+  return (Error_Code);
+} /* END OF Convert_Geodetic_To_Geocentric */
+
+/*
+ * The function Convert_Geocentric_To_Geodetic converts geocentric
+ * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude, 
+ * and height), according to the current ellipsoid parameters.
+ *
+ *    X         : Geocentric X coordinate, in meters.         (input)
+ *    Y         : Geocentric Y coordinate, in meters.         (input)
+ *    Z         : Geocentric Z coordinate, in meters.         (input)
+ *    Latitude  : Calculated latitude value in radians.       (output)
+ *    Longitude : Calculated longitude value in radians.      (output)
+ *    Height    : Calculated height value, in meters.         (output)
+ */
+
+#define USE_ITERATIVE_METHOD
+
+void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo *gi,
+                                        double X,
+                                        double Y, 
+                                        double Z,
+                                        double *Latitude,
+                                        double *Longitude,
+                                        double *Height)
+{ /* BEGIN Convert_Geocentric_To_Geodetic */
+#if !defined(USE_ITERATIVE_METHOD)
+/*
+ * The method used here is derived from 'An Improved Algorithm for
+ * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
+ */
+
+/* Note: Variable names follow the notation used in Toms, Feb 1996 */
+
+    double W;        /* distance from Z axis */
+    double W2;       /* square of distance from Z axis */
+    double T0;       /* initial estimate of vertical component */
+    double T1;       /* corrected estimate of vertical component */
+    double S0;       /* initial estimate of horizontal component */
+    double S1;       /* corrected estimate of horizontal component */
+    double Sin_B0;   /* sin(B0), B0 is estimate of Bowring aux variable */
+    double Sin3_B0;  /* cube of sin(B0) */
+    double Cos_B0;   /* cos(B0) */
+    double Sin_p1;   /* sin(phi1), phi1 is estimated latitude */
+    double Cos_p1;   /* cos(phi1) */
+    double Rn;       /* Earth radius at location */
+    double Sum;      /* numerator of cos(phi1) */
+    int At_Pole;     /* indicates location is in polar region */
+
+    At_Pole = FALSE;
+    if (X != 0.0)
+    {
+        *Longitude = atan2(Y,X);
+    }
+    else
+    {
+        if (Y > 0)
+        {
+            *Longitude = PI_OVER_2;
+        }
+        else if (Y < 0)
+        {
+            *Longitude = -PI_OVER_2;
+        }
+        else
+        {
+            At_Pole = TRUE;
+            *Longitude = 0.0;
+            if (Z > 0.0)
+            {  /* north pole */
+                *Latitude = PI_OVER_2;
+            }
+            else if (Z < 0.0)
+            {  /* south pole */
+                *Latitude = -PI_OVER_2;
+            }
+            else
+            {  /* center of earth */
+                *Latitude = PI_OVER_2;
+                *Height = -Geocent_b;
+                return;
+            } 
+        }
+    }
+    W2 = X*X + Y*Y;
+    W = sqrt(W2);
+    T0 = Z * AD_C;
+    S0 = sqrt(T0 * T0 + W2);
+    Sin_B0 = T0 / S0;
+    Cos_B0 = W / S0;
+    Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
+    T1 = Z + gi->Geocent_b * gi->Geocent_ep2 * Sin3_B0;
+    Sum = W - gi->Geocent_a * gi->Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0;
+    S1 = sqrt(T1*T1 + Sum * Sum);
+    Sin_p1 = T1 / S1;
+    Cos_p1 = Sum / S1;
+    Rn = gi->Geocent_a / sqrt(1.0 - gi->Geocent_e2 * Sin_p1 * Sin_p1);
+    if (Cos_p1 >= COS_67P5)
+    {
+        *Height = W / Cos_p1 - Rn;
+    }
+    else if (Cos_p1 <= -COS_67P5)
+    {
+        *Height = W / -Cos_p1 - Rn;
+    }
+    else
+    {
+        *Height = Z / Sin_p1 + Rn * (gi->Geocent_e2 - 1.0);
+    }
+    if (At_Pole == FALSE)
+    {
+        *Latitude = atan(Sin_p1 / Cos_p1);
+    }
+#else /* defined(USE_ITERATIVE_METHOD) */
+/*
+* Reference...
+* ============
+* Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für
+* das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover
+* Nr. 137, p. 130-131.
+
+* Programmed by GGA- Leibniz-Institute of Applied Geophysics
+*               Stilleweg 2
+*               D-30655 Hannover
+*               Federal Republic of Germany
+*               Internet: www.gga-hannover.de
+*
+*               Hannover, March 1999, April 2004.
+*               see also: comments in statements
+* remarks:
+* Mathematically exact and because of symmetry of rotation-ellipsoid,
+* each point (X,Y,Z) has at least two solutions (Latitude1,Longitude1,Height1) and
+* (Latitude2,Longitude2,Height2). Is point=(0.,0.,Z) (P=0.), so you get even
+* four solutions,	every two symmetrical to the semi-minor axis.
+* Here Height1 and Height2 have at least a difference in order of
+* radius of curvature (e.g. (0,0,b)=> (90.,0.,0.) or (-90.,0.,-2b);
+* (a+100.)*(sqrt(2.)/2.,sqrt(2.)/2.,0.) => (0.,45.,100.) or
+* (0.,225.,-(2a+100.))).
+* The algorithm always computes (Latitude,Longitude) with smallest |Height|.
+* For normal computations, that means |Height|<10000.m, algorithm normally
+* converges after to 2-3 steps!!!
+* But if |Height| has the amount of length of ellipsoid's axis
+* (e.g. -6300000.m),	algorithm needs about 15 steps.
+*/
+
+/* local definitions and variables */
+/* end-criterium of loop, accuracy of sin(Latitude) */
+#define genau   1.E-12
+#define genau2  (genau*genau)
+#define maxiter 30
+
+    double P;        /* distance between semi-minor axis and location */
+    double RR;       /* distance between center and location */
+    double CT;       /* sin of geocentric latitude */
+    double ST;       /* cos of geocentric latitude */
+    double RX;
+    double RK;
+    double RN;       /* Earth radius at location */
+    double CPHI0;    /* cos of start or old geodetic latitude in iterations */
+    double SPHI0;    /* sin of start or old geodetic latitude in iterations */
+    double CPHI;     /* cos of searched geodetic latitude */
+    double SPHI;     /* sin of searched geodetic latitude */
+    double SDPHI;    /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
+    int iter;        /* # of continuous iteration, max. 30 is always enough (s.a.) */
+
+    P = sqrt(X*X+Y*Y);
+    RR = sqrt(X*X+Y*Y+Z*Z);
+
+/*	special cases for latitude and longitude */
+    if (P/gi->Geocent_a < genau) {
+
+/*  special case, if P=0. (X=0., Y=0.) */
+	*Longitude = 0.;
+
+/*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
+ *  of ellipsoid (=center of mass), Latitude becomes PI/2 */
+        if (RR/gi->Geocent_a < genau) {
+            *Latitude = PI_OVER_2;
+            *Height   = -gi->Geocent_b;
+            return ;
+
+        }
+    }
+    else {
+/*  ellipsoidal (geodetic) longitude
+ *  interval: -PI < Longitude <= +PI */
+        *Longitude=atan2(Y,X);
+    }
+
+/* --------------------------------------------------------------
+ * Following iterative algorithm was developed by
+ * "Institut für Erdmessung", University of Hannover, July 1988.
+ * Internet: www.ife.uni-hannover.de
+ * Iterative computation of CPHI,SPHI and Height.
+ * Iteration of CPHI and SPHI to 10**-12 radian resp.
+ * 2*10**-7 arcsec.
+ * --------------------------------------------------------------
+ */
+    CT = Z/RR;
+    ST = P/RR;
+    {
+        const double denominator = 1.0-gi->Geocent_e2*(2.0-gi->Geocent_e2)*ST*ST;
+        if( denominator == 0 )
+        {
+            *Latitude = HUGE_VAL;
+            *Longitude = HUGE_VAL;
+            *Height = HUGE_VAL;
+            return;
+        }
+        RX = 1.0/sqrt(denominator);
+    }
+    CPHI0 = ST*(1.0-gi->Geocent_e2)*RX;
+    SPHI0 = CT*RX;
+    iter = 0;
+
+/* loop to find sin(Latitude) resp. Latitude
+ * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
+    do
+    {
+        iter++;
+        RN = gi->Geocent_a/sqrt(1.0-gi->Geocent_e2*SPHI0*SPHI0);
+
+/*  ellipsoidal (geodetic) height */
+        *Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi->Geocent_e2*SPHI0*SPHI0);
+
+        /* avoid zero division */
+        if (RN+*Height==0.0) {
+            *Latitude = 0.0;
+            return;
+        }
+        RK = gi->Geocent_e2*RN/(RN+*Height);
+        {
+            const double denominator = 1.0-RK*(2.0-RK)*ST*ST;
+            if( denominator == 0 )
+            {
+                *Latitude = HUGE_VAL;
+                *Longitude = HUGE_VAL;
+                *Height = HUGE_VAL;
+                return;
+            }
+            RX = 1.0/sqrt(denominator);
+        }
+        CPHI = ST*(1.0-RK)*RX;
+        SPHI = CT*RX;
+        SDPHI = SPHI*CPHI0-CPHI*SPHI0;
+        CPHI0 = CPHI;
+        SPHI0 = SPHI;
+    }
+    while (SDPHI*SDPHI > genau2 && iter < maxiter);
+
+/*	ellipsoidal (geodetic) latitude */
+    *Latitude=atan2(SPHI, fabs(CPHI));
+
+#endif /* defined(USE_ITERATIVE_METHOD) */
+} /* END OF Convert_Geocentric_To_Geodetic */
+
+
+///////////////////////////////////////////////////////////////////////////////
+/// Main of transform.cpp
+///////////////////////////////////////////////////////////////////////////////
+
+
 static int adjust_axis( projCtx ctx, const char *axis, int denormalize_flag,
                            long point_count, int point_offset,
                            double *x, double *y, double *z );
diff --git a/src/transformations/helmert.cpp b/src/transformations/helmert.cpp
index 2045dea2..33a703bb 100644
--- a/src/transformations/helmert.cpp
+++ b/src/transformations/helmert.cpp
@@ -53,7 +53,6 @@ Last update: 2018-10-26
 #include <math.h>
 
 #include "proj_internal.h"
-#include "geocent.h"
 
 PROJ_HEAD(helmert, "3(6)-, 4(8)- and 7(14)-parameter Helmert shift");
 PROJ_HEAD(molobadekas, "Molodensky-Badekas transformation");