aboutsummaryrefslogtreecommitdiff
diff options
context:
space:
mode:
Diffstat
-rw-r--r--src/Makefile.am2
-rw-r--r--src/geocent.cpp456
-rw-r--r--src/geocent.h179
-rw-r--r--src/lib_proj.cmake2
-rw-r--r--src/projections/sch.cpp139
-rw-r--r--src/transform.cpp651
-rw-r--r--src/transformations/helmert.cpp1
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");
generated by cgit v1.2.3 (git 2.25.1) at 2025-12-16 04:55:42 +0000