projections/impl/geocent.hpp

0001 // Boost.Geometry
0002 // This file is manually converted from PROJ4
0003
0004 // This file was modified by Oracle on 2017.
0005 // Modifications copyright (c) 2017, Oracle and/or its affiliates.
0006 // Contributed and/or modified by Adam Wulkiewicz, on behalf of Oracle
0007
0008 // Use, modification and distribution is subject to the Boost Software License,
0009 // Version 1.0. (See accompanying file LICENSE_1_0.txt or copy at
0010 // http://www.boost.org/LICENSE_1_0.txt)
0011
0012 // This file is converted from PROJ4, http://trac.osgeo.org/proj
0013 // PROJ4 is originally written by Gerald Evenden (then of the USGS)
0014 // PROJ4 is maintained by Frank Warmerdam
0015 // This file was converted to Geometry Library by Adam Wulkiewicz
0016
0017 // Original copyright notice:
0018
0019 /***************************************************************************/
0020 /* RSC IDENTIFIER:  GEOCENTRIC
0021  *
0022  * ABSTRACT
0023  *
0024  *    This component provides conversions between Geodetic coordinates (latitude,
0025  *    longitude in radians and height in meters) and Geocentric coordinates
0026  *    (X, Y, Z) in meters.
0027  *
0028  * ERROR HANDLING
0029  *
0030  *    This component checks parameters for valid values.  If an invalid value
0031  *    is found, the error code is combined with the current error code using
0032  *    the bitwise or.  This combining allows multiple error codes to be
0033  *    returned. The possible error codes are:
0034  *
0035  *      GEOCENT_NO_ERROR        : No errors occurred in function
0036  *      GEOCENT_LAT_ERROR       : Latitude out of valid range
0037  *                                 (-90 to 90 degrees)
0038  *      GEOCENT_LON_ERROR       : Longitude out of valid range
0039  *                                 (-180 to 360 degrees)
0040  *      GEOCENT_A_ERROR         : Semi-major axis lessthan or equal to zero
0041  *      GEOCENT_B_ERROR         : Semi-minor axis lessthan or equal to zero
0042  *      GEOCENT_A_LESS_B_ERROR  : Semi-major axis less than semi-minor axis
0043  *
0044  *
0045  * REUSE NOTES
0046  *
0047  *    GEOCENTRIC is intended for reuse by any application that performs
0048  *    coordinate conversions between geodetic coordinates and geocentric
0049  *    coordinates.
0050  *
0051  *
0052  * REFERENCES
0053  *
0054  *    An Improved Algorithm for Geocentric to Geodetic Coordinate Conversion,
0055  *    Ralph Toms, February 1996  UCRL-JC-123138.
0056  *
0057  *    Further information on GEOCENTRIC can be found in the Reuse Manual.
0058  *
0059  *    GEOCENTRIC originated from : U.S. Army Topographic Engineering Center
0060  *                                 Geospatial Information Division
0061  *                                 7701 Telegraph Road
0062  *                                 Alexandria, VA  22310-3864
0063  *
0064  * LICENSES
0065  *
0066  *    None apply to this component.
0067  *
0068  * RESTRICTIONS
0069  *
0070  *    GEOCENTRIC has no restrictions.
0071  *
0072  * ENVIRONMENT
0073  *
0074  *    GEOCENTRIC was tested and certified in the following environments:
0075  *
0076  *    1. Solaris 2.5 with GCC version 2.8.1
0077  *    2. Windows 95 with MS Visual C++ version 6
0078  *
0079  * MODIFICATIONS
0080  *
0081  *    Date              Description
0082  *    ----              -----------
0083  *    25-02-97          Original Code
0084  *
0085  */
0086
0087
0088 #ifndef BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_GEOCENT_HPP
0089 #define BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_GEOCENT_HPP
0090
0091
0092 #include <boost/geometry/util/math.hpp>
0093
0094
0095 namespace boost { namespace geometry { namespace projections
0096 {
0097
0098 namespace detail
0099 {
0100
0101 /***************************************************************************/
0102 /*
0103  *                               DEFINES
0104  */
0105 static const long GEOCENT_NO_ERROR       = 0x0000;
0106 static const long GEOCENT_LAT_ERROR      = 0x0001;
0107 static const long GEOCENT_LON_ERROR      = 0x0002;
0108 static const long GEOCENT_A_ERROR        = 0x0004;
0109 static const long GEOCENT_B_ERROR        = 0x0008;
0110 static const long GEOCENT_A_LESS_B_ERROR = 0x0010;
0111
0112 template <typename T>
0113 struct GeocentricInfo
0114 {
0115     T Geocent_a;        /* Semi-major axis of ellipsoid in meters */
0116     T Geocent_b;        /* Semi-minor axis of ellipsoid           */
0117     T Geocent_a2;       /* Square of semi-major axis */
0118     T Geocent_b2;       /* Square of semi-minor axis */
0119     T Geocent_e2;       /* Eccentricity squared  */
0120     T Geocent_ep2;      /* 2nd eccentricity squared */
0121 };
0122
0123 template <typename T>
0124 inline T COS_67P5()
0125 {
0126     /*return 0.38268343236508977*/;
0127     return cos(T(67.5) * math::d2r<T>());  /* cosine of 67.5 degrees */
0128 }
0129 template <typename T>
0130 inline T AD_C()
0131 {
0132     return 1.0026000;            /* Toms region 1 constant */
0133 }
0134
0135
0136 /***************************************************************************/
0137 /*
0138  *                              FUNCTIONS
0139  */
0140
0141 template <typename T>
0142 inline long pj_Set_Geocentric_Parameters (GeocentricInfo<T> & gi, T const& a, T const& b)
0143
0144 { /* BEGIN Set_Geocentric_Parameters */
0145 /*
0146  * The function Set_Geocentric_Parameters receives the ellipsoid parameters
0147  * as inputs and sets the corresponding state variables.
0148  *
0149  *    a  : Semi-major axis, in meters.          (input)
0150  *    b  : Semi-minor axis, in meters.          (input)
0151  */
0152     long Error_Code = GEOCENT_NO_ERROR;
0153
0154     if (a <= 0.0)
0155         Error_Code |= GEOCENT_A_ERROR;
0156     if (b <= 0.0)
0157         Error_Code |= GEOCENT_B_ERROR;
0158     if (a < b)
0159         Error_Code |= GEOCENT_A_LESS_B_ERROR;
0160     if (!Error_Code)
0161     {
0162         gi.Geocent_a = a;
0163         gi.Geocent_b = b;
0164         gi.Geocent_a2 = a * a;
0165         gi.Geocent_b2 = b * b;
0166         gi.Geocent_e2 = (gi.Geocent_a2 - gi.Geocent_b2) / gi.Geocent_a2;
0167         gi.Geocent_ep2 = (gi.Geocent_a2 - gi.Geocent_b2) / gi.Geocent_b2;
0168     }
0169     return (Error_Code);
0170 } /* END OF Set_Geocentric_Parameters */
0171
0172
0173 template <typename T>
0174 inline void pj_Get_Geocentric_Parameters (GeocentricInfo<T> const& gi,
0175                                           T & a,
0176                                           T & b)
0177 { /* BEGIN Get_Geocentric_Parameters */
0178 /*
0179  * The function Get_Geocentric_Parameters returns the ellipsoid parameters
0180  * to be used in geocentric coordinate conversions.
0181  *
0182  *    a  : Semi-major axis, in meters.          (output)
0183  *    b  : Semi-minor axis, in meters.          (output)
0184  */
0185
0186     a = gi.Geocent_a;
0187     b = gi.Geocent_b;
0188 } /* END OF Get_Geocentric_Parameters */
0189
0190
0191 template <typename T>
0192 inline long pj_Convert_Geodetic_To_Geocentric (GeocentricInfo<T> const& gi,
0193                                                T Longitude, T Latitude, T Height,
0194                                                T & X, T & Y, T & Z)
0195 { /* BEGIN Convert_Geodetic_To_Geocentric */
0196 /*
0197  * The function Convert_Geodetic_To_Geocentric converts geodetic coordinates
0198  * (latitude, longitude, and height) to geocentric coordinates (X, Y, Z),
0199  * according to the current ellipsoid parameters.
0200  *
0201  *    Latitude  : Geodetic latitude in radians                     (input)
0202  *    Longitude : Geodetic longitude in radians                    (input)
0203  *    Height    : Geodetic height, in meters                       (input)
0204  *    X         : Calculated Geocentric X coordinate, in meters    (output)
0205  *    Y         : Calculated Geocentric Y coordinate, in meters    (output)
0206  *    Z         : Calculated Geocentric Z coordinate, in meters    (output)
0207  *
0208  */
0209   long Error_Code = GEOCENT_NO_ERROR;
0210   T Rn;            /*  Earth radius at location  */
0211   T Sin_Lat;       /*  sin(Latitude)  */
0212   T Sin2_Lat;      /*  Square of sin(Latitude)  */
0213   T Cos_Lat;       /*  cos(Latitude)  */
0214
0215   static const T PI = math::pi<T>();
0216   static const T PI_OVER_2 = math::half_pi<T>();
0217
0218   /*
0219   ** Don't blow up if Latitude is just a little out of the value
0220   ** range as it may just be a rounding issue.  Also removed longitude
0221   ** test, it should be wrapped by cos() and sin().  NFW for PROJ.4, Sep/2001.
0222   */
0223   if( Latitude < -PI_OVER_2 && Latitude > -1.001 * PI_OVER_2 )
0224       Latitude = -PI_OVER_2;
0225   else if( Latitude > PI_OVER_2 && Latitude < 1.001 * PI_OVER_2 )
0226       Latitude = PI_OVER_2;
0227   else if ((Latitude < -PI_OVER_2) || (Latitude > PI_OVER_2))
0228   { /* Latitude out of range */
0229     Error_Code |= GEOCENT_LAT_ERROR;
0230   }
0231
0232   if (!Error_Code)
0233   { /* no errors */
0234     if (Longitude > PI)
0235       Longitude -= (2*PI);
0236     Sin_Lat = sin(Latitude);
0237     Cos_Lat = cos(Latitude);
0238     Sin2_Lat = Sin_Lat * Sin_Lat;
0239     Rn = gi.Geocent_a / (sqrt(1.0e0 - gi.Geocent_e2 * Sin2_Lat));
0240     X = (Rn + Height) * Cos_Lat * cos(Longitude);
0241     Y = (Rn + Height) * Cos_Lat * sin(Longitude);
0242     Z = ((Rn * (1 - gi.Geocent_e2)) + Height) * Sin_Lat;
0243   }
0244   return (Error_Code);
0245 } /* END OF Convert_Geodetic_To_Geocentric */
0246
0247 /*
0248  * The function Convert_Geocentric_To_Geodetic converts geocentric
0249  * coordinates (X, Y, Z) to geodetic coordinates (latitude, longitude,
0250  * and height), according to the current ellipsoid parameters.
0251  *
0252  *    X         : Geocentric X coordinate, in meters.         (input)
0253  *    Y         : Geocentric Y coordinate, in meters.         (input)
0254  *    Z         : Geocentric Z coordinate, in meters.         (input)
0255  *    Latitude  : Calculated latitude value in radians.       (output)
0256  *    Longitude : Calculated longitude value in radians.      (output)
0257  *    Height    : Calculated height value, in meters.         (output)
0258  */
0259
0260 #define BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD
0261
0262 template <typename T>
0263 inline void pj_Convert_Geocentric_To_Geodetic (GeocentricInfo<T> const& gi,
0264                                                T X, T Y, T Z,
0265                                                T & Longitude, T & Latitude, T & Height)
0266 { /* BEGIN Convert_Geocentric_To_Geodetic */
0267
0268     static const T PI_OVER_2 = math::half_pi<T>();
0269
0270 #if !defined(BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD)
0271
0272     static const T COS_67P5 = detail::COS_67P5<T>();
0273     static const T AD_C = detail::AD_C<T>();
0274
0275 /*
0276  * The method used here is derived from 'An Improved Algorithm for
0277  * Geocentric to Geodetic Coordinate Conversion', by Ralph Toms, Feb 1996
0278  */
0279
0280 /* Note: Variable names follow the notation used in Toms, Feb 1996 */
0281
0282     T W;        /* distance from Z axis */
0283     T W2;       /* square of distance from Z axis */
0284     T T0;       /* initial estimate of vertical component */
0285     T T1;       /* corrected estimate of vertical component */
0286     T S0;       /* initial estimate of horizontal component */
0287     T S1;       /* corrected estimate of horizontal component */
0288     T Sin_B0;   /* sin(B0), B0 is estimate of Bowring aux variable */
0289     T Sin3_B0;  /* cube of sin(B0) */
0290     T Cos_B0;   /* cos(B0) */
0291     T Sin_p1;   /* sin(phi1), phi1 is estimated latitude */
0292     T Cos_p1;   /* cos(phi1) */
0293     T Rn;       /* Earth radius at location */
0294     T Sum;      /* numerator of cos(phi1) */
0295     bool At_Pole;     /* indicates location is in polar region */
0296
0297     At_Pole = false;
0298     if (X != 0.0)
0299     {
0300         Longitude = atan2(Y,X);
0301     }
0302     else
0303     {
0304         if (Y > 0)
0305         {
0306             Longitude = PI_OVER_2;
0307         }
0308         else if (Y < 0)
0309         {
0310             Longitude = -PI_OVER_2;
0311         }
0312         else
0313         {
0314             At_Pole = true;
0315             Longitude = 0.0;
0316             if (Z > 0.0)
0317             {  /* north pole */
0318                 Latitude = PI_OVER_2;
0319             }
0320             else if (Z < 0.0)
0321             {  /* south pole */
0322                 Latitude = -PI_OVER_2;
0323             }
0324             else
0325             {  /* center of earth */
0326                 Latitude = PI_OVER_2;
0327                 Height = -Geocent_b;
0328                 return;
0329             }
0330         }
0331     }
0332     W2 = X*X + Y*Y;
0333     W = sqrt(W2);
0334     T0 = Z * AD_C;
0335     S0 = sqrt(T0 * T0 + W2);
0336     Sin_B0 = T0 / S0;
0337     Cos_B0 = W / S0;
0338     Sin3_B0 = Sin_B0 * Sin_B0 * Sin_B0;
0339     T1 = Z + gi.Geocent_b * gi.Geocent_ep2 * Sin3_B0;
0340     Sum = W - gi.Geocent_a * gi.Geocent_e2 * Cos_B0 * Cos_B0 * Cos_B0;
0341     S1 = sqrt(T1*T1 + Sum * Sum);
0342     Sin_p1 = T1 / S1;
0343     Cos_p1 = Sum / S1;
0344     Rn = gi.Geocent_a / sqrt(1.0 - gi.Geocent_e2 * Sin_p1 * Sin_p1);
0345     if (Cos_p1 >= COS_67P5)
0346     {
0347         Height = W / Cos_p1 - Rn;
0348     }
0349     else if (Cos_p1 <= -COS_67P5)
0350     {
0351         Height = W / -Cos_p1 - Rn;
0352     }
0353     else
0354     {
0355         Height = Z / Sin_p1 + Rn * (gi.Geocent_e2 - 1.0);
0356     }
0357     if (At_Pole == false)
0358     {
0359         Latitude = atan(Sin_p1 / Cos_p1);
0360     }
0361 #else /* defined(BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD) */
0362 /*
0363 * Reference...
0364 * ============
0365 * Wenzel, H.-G.(1985): Hochauflösende Kugelfunktionsmodelle für
0366 * das Gravitationspotential der Erde. Wiss. Arb. Univ. Hannover
0367 * Nr. 137, p. 130-131.
0368
0369 * Programmed by GGA- Leibniz-Institute of Applied Geophysics
0370 *               Stilleweg 2
0371 *               D-30655 Hannover
0372 *               Federal Republic of Germany
0373 *               Internet: www.gga-hannover.de
0374 *
0375 *               Hannover, March 1999, April 2004.
0376 *               see also: comments in statements
0377 * remarks:
0378 * Mathematically exact and because of symmetry of rotation-ellipsoid,
0379 * each point (X,Y,Z) has at least two solutions (Latitude1,Longitude1,Height1) and
0380 * (Latitude2,Longitude2,Height2). Is point=(0.,0.,Z) (P=0.), so you get even
0381 * four solutions,   every two symmetrical to the semi-minor axis.
0382 * Here Height1 and Height2 have at least a difference in order of
0383 * radius of curvature (e.g. (0,0,b)=> (90.,0.,0.) or (-90.,0.,-2b);
0384 * (a+100.)*(sqrt(2.)/2.,sqrt(2.)/2.,0.) => (0.,45.,100.) or
0385 * (0.,225.,-(2a+100.))).
0386 * The algorithm always computes (Latitude,Longitude) with smallest |Height|.
0387 * For normal computations, that means |Height|<10000.m, algorithm normally
0388 * converges after to 2-3 steps!!!
0389 * But if |Height| has the amount of length of ellipsoid's axis
0390 * (e.g. -6300000.m),    algorithm needs about 15 steps.
0391 */
0392
0393 /* local definitions and variables */
0394 /* end-criterium of loop, accuracy of sin(Latitude) */
0395 static const T genau   = 1.E-12;
0396 static const T genau2  = (genau*genau);
0397 static const int maxiter = 30;
0398
0399     T P;        /* distance between semi-minor axis and location */
0400     T RR;       /* distance between center and location */
0401     T CT;       /* sin of geocentric latitude */
0402     T ST;       /* cos of geocentric latitude */
0403     T RX;
0404     T RK;
0405     T RN;       /* Earth radius at location */
0406     T CPHI0;    /* cos of start or old geodetic latitude in iterations */
0407     T SPHI0;    /* sin of start or old geodetic latitude in iterations */
0408     T CPHI;     /* cos of searched geodetic latitude */
0409     T SPHI;     /* sin of searched geodetic latitude */
0410     T SDPHI;    /* end-criterium: addition-theorem of sin(Latitude(iter)-Latitude(iter-1)) */
0411     int iter;   /* # of continuous iteration, max. 30 is always enough (s.a.) */
0412
0413     P = sqrt(X*X+Y*Y);
0414     RR = sqrt(X*X+Y*Y+Z*Z);
0415
0416 /*  special cases for latitude and longitude */
0417     if (P/gi.Geocent_a < genau) {
0418
0419 /*  special case, if P=0. (X=0., Y=0.) */
0420     Longitude = 0.;
0421
0422 /*  if (X,Y,Z)=(0.,0.,0.) then Height becomes semi-minor axis
0423  *  of ellipsoid (=center of mass), Latitude becomes PI/2 */
0424         if (RR/gi.Geocent_a < genau) {
0425             Latitude = PI_OVER_2;
0426             Height   = -gi.Geocent_b;
0427             return ;
0428
0429         }
0430     }
0431     else {
0432 /*  ellipsoidal (geodetic) longitude
0433  *  interval: -PI < Longitude <= +PI */
0434         Longitude=atan2(Y,X);
0435     }
0436
0437 /* --------------------------------------------------------------
0438  * Following iterative algorithm was developed by
0439  * "Institut für Erdmessung", University of Hannover, July 1988.
0440  * Internet: www.ife.uni-hannover.de
0441  * Iterative computation of CPHI,SPHI and Height.
0442  * Iteration of CPHI and SPHI to 10**-12 radian resp.
0443  * 2*10**-7 arcsec.
0444  * --------------------------------------------------------------
0445  */
0446     CT = Z/RR;
0447     ST = P/RR;
0448     RX = 1.0/sqrt(1.0-gi.Geocent_e2*(2.0-gi.Geocent_e2)*ST*ST);
0449     CPHI0 = ST*(1.0-gi.Geocent_e2)*RX;
0450     SPHI0 = CT*RX;
0451     iter = 0;
0452
0453 /* loop to find sin(Latitude) resp. Latitude
0454  * until |sin(Latitude(iter)-Latitude(iter-1))| < genau */
0455     do
0456     {
0457         iter++;
0458         RN = gi.Geocent_a/sqrt(1.0-gi.Geocent_e2*SPHI0*SPHI0);
0459
0460 /*  ellipsoidal (geodetic) height */
0461         Height = P*CPHI0+Z*SPHI0-RN*(1.0-gi.Geocent_e2*SPHI0*SPHI0);
0462
0463         RK = gi.Geocent_e2*RN/(RN+Height);
0464         RX = 1.0/sqrt(1.0-RK*(2.0-RK)*ST*ST);
0465         CPHI = ST*(1.0-RK)*RX;
0466         SPHI = CT*RX;
0467         SDPHI = SPHI*CPHI0-CPHI*SPHI0;
0468         CPHI0 = CPHI;
0469         SPHI0 = SPHI;
0470     }
0471     while (SDPHI*SDPHI > genau2 && iter < maxiter);
0472
0473 /*  ellipsoidal (geodetic) latitude */
0474     Latitude=atan(SPHI/fabs(CPHI));
0475
0476     return;
0477 #endif /* defined(BOOST_GEOMETRY_PROJECTIONS_USE_ITERATIVE_METHOD) */
0478 } /* END OF Convert_Geocentric_To_Geodetic */
0479
0480
0481 } // namespace detail
0482
0483
0484 }}} // namespace boost::geometry::projections
0485
0486
0487 #endif // BOOST_GEOMETRY_SRS_PROJECTIONS_IMPL_GEOCENT_HPP