src/common/greatcircle.cpp

   1 #include "greatcircle.h"
   2
   3 #define DELTA 1e-5
   4
   5 static bool antipodes(const Coordinates &c1, const Coordinates &c2)
   6 {
   7         return ((qAbs(c1.lat() + c2.lat()) < DELTA)
   8           && (qAbs(180.0 - qAbs(c1.lon() - c2.lon())) < DELTA));
   9 }
  10
  11 GreatCircle::GreatCircle(const Coordinates &c1, const Coordinates &c2)
  12 {
  13         double lat1, lon1, lat2, lon2;
  14
  15         if (antipodes(c1, c2)) {
  16                 /* In case of antipodes (which would lead to garbage output without
  17                    this hack), move the points DELTA degrees closer to each other in
  18                    a way that the route never crosses the antimeridian. */
  19                 if (c1.lon() < c2.lon()) {
  20                         lon1 = deg2rad(c1.lon() + DELTA);
  21                         lon2 = deg2rad(c2.lon() - DELTA);
  22                 } else {
  23                         lon1 = deg2rad(c1.lon() - DELTA);
  24                         lon2 = deg2rad(c2.lon() + DELTA);
  25                 }
  26                 if (c1.lat() < c2.lat()) {
  27                         lat1 = deg2rad(c1.lat() + DELTA);
  28                         lat2 = deg2rad(c2.lat() - DELTA);
  29                 } else {
  30                         lat1 = deg2rad(c1.lat() - DELTA);
  31                         lat2 = deg2rad(c2.lat() + DELTA);
  32                 }
  33         } else {
  34                 lat1 = deg2rad(c1.lat());
  35                 lon1 = deg2rad(c1.lon());
  36                 lat2 = deg2rad(c2.lat());
  37                 lon2 = deg2rad(c2.lon());
  38         }
  39
  40         double cosLat1 = cos(lat1);
  41         double cosLat2 = cos(lat2);
  42         _sinLat1 = sin(lat1);
  43         _sinLat2 = sin(lat2);
  44
  45         _xA = cosLat1 * cos(lon1);
  46         _xB = cosLat2 * cos(lon2);
  47         _yA = cosLat1 * sin(lon1);
  48         _yB = cosLat2 * sin(lon2);
  49
  50         _d = 2 * asin(sqrt(pow((sin((lat1 - lat2) / 2)), 2) + cosLat1 * cosLat2
  51           * pow(sin((lon1 - lon2) / 2), 2)));
  52         _sinD = sin(_d);
  53 }
  54
  55 Coordinates GreatCircle::pointAt(double f) const
  56 {
  57         double A = sin((1.0 - f) * _d) / _sinD;
  58         double B = sin(f * _d) / _sinD;
  59         double x = A * _xA + B * _xB;
  60         double y = A * _yA + B * _yB;
  61         double z = A * _sinLat1 + B * _sinLat2;
  62
  63         return Coordinates(rad2deg(atan2(y, x)), rad2deg(atan2(z, sqrt(x*x + y*y))));
  64 }