src/FLoc_EllipsoidModel.cpp

   1 //
   2 // Open Service Platform
   3 // Copyright (c) 2012 Samsung Electronics Co., Ltd.
   4 //
   5 // Licensed under the Apache License, Version 2.0 (the License);
   6 // you may not use this file except in compliance with the License.
   7 // You may obtain a copy of the License at
   8 //
   9 //     http://www.apache.org/licenses/LICENSE-2.0
  10 //
  11 // Unless required by applicable law or agreed to in writing, software
  12 // distributed under the License is distributed on an "AS IS" BASIS,
  13 // WITHOUT WARRANTIES OR CONDITIONS OF ANY KIND, either express or implied.
  14 // See the License for the specific language governing permissions and
  15 // limitations under the License.
  16 //
  17
  18 /**
  19  * @file        FLoc_EllipsoidModel.h
  20  * @brief       This is the implementation file for the _EllipsoidModel class.
  21  *
  22  * This header file contains implementation of the _EllipsoidModel class.
  23  */
  24
  25 #include <math.h>
  26 #include <FBaseDouble.h>
  27 #include "FLoc_EllipsoidModel.h"
  28 #include "FLoc_MathUtils.h"
  29
  30 using namespace Tizen::Base;
  31
  32 namespace Tizen { namespace Locations
  33 {
  34
  35 const double _EllipsoidModel::SEMI_MAJOR_AXIS = 6378137.0;         //in meter
  36 const double _EllipsoidModel::SEMI_MINOR_AXIS = 6356752.3142;      // in meter
  37 const double _EllipsoidModel::FLATTENING = 298.257223563;
  38 const double _EllipsoidModel::EARTH_RADIUS = _EllipsoidModel::SEMI_MAJOR_AXIS;
  39
  40 float
  41 _EllipsoidModel::GetDistance(double lat1, double lon1, double lat2, double lon2)
  42 {
  43         return SolveInverseProblem(lat1, lon1, lat2, lon2, RT_DISTANCE);
  44 }
  45
  46 float
  47 _EllipsoidModel::GetAzimuth(double lat1, double lon1, double lat2, double lon2)
  48 {
  49         return SolveInverseProblem(lat1, lon1, lat2, lon2, RT_FORWARD_AZIMUTH);
  50 }
  51
  52 float
  53 _EllipsoidModel::SolveInverseProblem(double lat1, double lon1, double lat2, double lon2, ResultType type)
  54 {
  55         // References: Direct and inverse solutions of geodesics on the ellipsoid with application of nested equations
  56         // The section 4, Inverse Formula, is used.
  57         // Instead of the equations (3), (4), and (6), the simplified ones (3a), (4a), and (6a) are adopted.
  58
  59         const int ITERATION_MAX = 10;
  60         const double DELTA_TOLERANCE = 1.0e-12;
  61         const double a = EARTH_RADIUS; // WGS-84 major axis (6378137.0)
  62         const double b = SEMI_MINOR_AXIS; // WGS-84 semi-minor axis
  63         const double f = 1.0 / FLATTENING; // inverse flattening
  64
  65         double L = (_MathUtils::DEG2RAD* lon2) -(_MathUtils::DEG2RAD* lon1);     // difference in longitude, positive east.
  66         double U1 = atan((1.0 - f) * tan(_MathUtils::DEG2RAD * lat1)); // reduced latitude
  67         double U2 = atan((1.0 - f) * tan(_MathUtils::DEG2RAD * lat2)); // reduced latitude
  68
  69         double cosU1 = cos(U1);
  70         double cosU2 = cos(U2);
  71         double sinU1 = sin(U1);
  72         double sinU2 = sin(U2);
  73         double cosU1cosU2 = cosU1 * cosU2;
  74         double sinU1sinU2 = sinU1 * sinU2;
  75
  76         double sigma = 0.0; // angular distance P1, P2 on the sphere
  77         double lambda = L; // difference in longitude on an auxiliary sphere. use L for first approximation.
  78         double A = 0.0, B = 0.0, C = 0.0, sinSigma = 0.0, cosSigma = 0.0, cos2SM = 0.0;
  79         int i = 0;
  80
  81         while (i < ITERATION_MAX)
  82         {
  83                 double sinAlpha = 0.0, cosSqAlpha = 0.0, uSq = 0.0, delta = 0.0;
  84                 double lambdaOld = lambda;
  85
  86                 double cosLambda = cos(lambda);
  87                 double sinLambda = sin(lambda);
  88
  89                 double t1 = cosU2 * sinLambda;
  90                 double t2 = cosU1 * sinU2 - sinU1 * cosU2 * cosLambda;
  91                 sinSigma = sqrt(t1 * t1 + t2 * t2); // (14)
  92                 cosSigma = sinU1sinU2 + cosU1cosU2 * cosLambda; // (15)
  93
  94                 sigma = atan2(sinSigma, cosSigma); // (16)
  95
  96                 sinAlpha = Double::Compare(sinSigma, 0.0) == 0 ? 0.0 : cosU1cosU2 * sinLambda / sinSigma; // (17)
  97                 cosSqAlpha = 1.0 - (sinAlpha * sinAlpha);
  98                 cos2SM = Double::Compare(cosSqAlpha, 0.0) == 0 ? 0.0 : cosSigma - (2.0 * sinU1sinU2 / cosSqAlpha); // (18)
  99
 100                 uSq = cosSqAlpha * ((a * a - b * b) / (b * b)); // u^2 = cos(alpha)^2 * (a^2 = b^2) / b^2
 101
 102                 A = 1.0 + (uSq / 256.0) * (64.0 + uSq * (-12.0 - (5.0 * uSq))); // (3a)
 103                 B = (uSq / 512.0) * (128.0 + uSq * (-64.0 + (37.0 * uSq))); // (4a)
 104                 C = (f / 16.0) * cosSqAlpha * (4.0 + f * (4.0 - 3.0 * cosSqAlpha)); // (10)
 105
 106                 lambda = L + ((1.0 - C) * f * sinAlpha * (sigma + C * sinSigma * (cos2SM + (C * cosSigma * (-1.0 + 2.0 * cos2SM * cos2SM))))); // from (11)
 107
 108                 delta = (lambda - lambdaOld) / lambda;
 109                 if (fabs(delta) < DELTA_TOLERANCE)
 110                 {
 111                         break;
 112                 }
 113                 ++i;
 114         }
 115
 116         float ret = 0.0F;
 117         if (RT_DISTANCE == type)
 118         {
 119                 double deltaSigma = B * sinSigma * (cos2SM + (B / 4.0) * (cosSigma * (-1.0 + 2.0 * cos2SM * cos2SM))); // (6a)
 120
 121                 ret = b * A * (sigma - deltaSigma);
 122
 123         }
 124         else if (RT_FORWARD_AZIMUTH == type)
 125         {
 126                 double cosLambda = cos(lambda);
 127                 double sinLambda = sin(lambda);
 128
 129                 ret = atan2(cosU2 * sinLambda, cosU1 * sinU2 - sinU1 * cosU2 * cosLambda);  // (20)
 130
 131                 if (ret < 0.0F)
 132                 {
 133                         ret += _MathUtils::PI2;
 134                 }
 135                 ret *= _MathUtils::RAD2DEG;
 136         }
 137
 138         return ret;
 139 }
 140
 141 } } // Tizen::Locations