include/ecef-wgs84.h

   1 /*
   2   Earth Centered Earth Fixed (ECEF) WGS84 Navigation Utility Copyright © 2013-2017  Infinite Delta Corp
   3
   4   This program is free software: you can redistribute it and/or modify
   5   it under the terms of the GNU Lesser General Public License as published by
   6   the Free Software Foundation, either version 2 of the License, or
   7   (at your option) any later version.
   8
   9   This program is distributed in the hope that it will be useful,
  10   but WITHOUT ANY WARRANTY; without even the implied warranty of
  11   MERCHANTABILITY or FITNESS FOR A PARTICULAR PURPOSE.  See the
  12   GNU General Public License for more details.
  13
  14   You should have received a copy of the GNU Lesser General Public License
  15   along with this program.  If not, see <http://www.gnu.org/licenses/>.
  16
  17 */
  18
  19 #ifndef ECEF_WGS84_H
  20 #define ECEF_WGS84_H
  21
  22 #include <ecef-vector.h>
  23
  24 // http://en.wikipedia.org/wiki/Geodetic_system
  25 // wgs84fin.pdf -> TR8350.2
  26 // http://proceedings.esri.com/library/userconf/proc14/papers/1034_755.pdf
  27
  28 // http://www.ngs.noaa.gov/AERO/aero.html
  29 // http://geodesy.noaa.gov/AERO/UDDFdat.htm
  30 // http://www.ngs.noaa.gov/AERO/uddf/GREAT-LAKES/MICHIGAN/GRR__01B.F77
  31
  32 namespace ecef
  33 {
  34   namespace wgs84Const
  35     {
  36     const double a                      = 6378137.0;            // Semi-major Axis Meters
  37     const double fr                     = 298.257223563;        // Flattening recipical
  38
  39     const double f                      = 1/fr;                 // Flattening
  40     const double c                      = a*(1-f);              // Minor Axis Meters
  41     const double c2                     = c*c;                  // Square Minor Axis M^2
  42     const double e2                     = 2*f-f*f;              // First Eccentricity Squared
  43     const double EarthCenterScale       = 1-a*a/c2;             // WGS84 Earth Center Z scaler
  44     const double maxEast                = M_PIl;
  45     const double maxWest                = -M_PIl;
  46     const double maxNorth               = M_PI_2l;
  47     const double maxSouth               = -M_PI_2l;
  48     }
  49
  50 struct wgs84 : public vec
  51   {
  52   wgs84() : vec(wgs84Const::a,0,0) {}
  53   wgs84(const vec v) : vec(v) {}
  54   wgs84(const double lat, const double lon, const double h=0.0)
  55     {
  56     double slat = ::sin(lat);
  57     double clat = ::cos(lat);
  58     double n = wgs84Const::a/sqrt(1-wgs84Const::e2*sqr(slat));
  59     x = (n+h)*clat*::cos(lon);
  60     y = (n+h)*clat*::sin(lon);
  61     z = (n*(1-wgs84Const::e2)+h)*slat;
  62     }
  63   // wgs84& operator=(const vec v) {x=v.x; y=v.y; z=v.z; return *this;}
  64
  65   double lon() const {double r=distZ(); return tolerance(r) ? 0 : ::acos(x/r)*sign(y);}
  66   };
  67
  68 struct wgs84Line : public line
  69   {
  70   wgs84Line() : line(wgs84()) {}
  71   wgs84Line(const vec end) : line(end) {}
  72   wgs84Line(const wgs84 end) : line(end) {}
  73   wgs84Line(const wgs84 start, const wgs84 end) : line(start, end) {}
  74   int trimToEllipsoid();
  75   };
  76
  77 // Fast: +/- 0.0004% @0m, or 0.2% @10,000m altitude
  78 //   Assuming altitude is zero (0) or insignificant, everything can be based on WGS84 earth center.
  79
  80 struct wgs84Fast : public wgs84
  81   {
  82   wgs84Fast() : wgs84() {}
  83   wgs84Fast(const wgs84 p) : wgs84(p) {}
  84   wgs84Fast(const double lat, const double lon, const double h=0.0) : wgs84(lat, lon, h) {}
  85
  86   wgs84 center() const {return vec(0,0,z*wgs84Const::EarthCenterScale);} // gravity earth center
  87   double tZ() const {return z * (1 - wgs84Const::EarthCenterScale);}
  88   wgs84 up() const {return vec(x,y,tZ());}
  89   wgs84 north() const {return zAxis.doubleCross(up());} // North tangent vector
  90   double sinLat() const {return tZ() / sqrt(sqr(x) + sqr(y) + sqr(tZ()));}
  91   double lat() const {return ::asin(sinLat());}
  92   };
  93 }
  94
  95 #endif