Geodesic.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_GEODESIC_HPP)
#define GEOGRAPHICLIB_GEODESIC_HPP 1

#include <GeographicLib/Constants.hpp>

#if !defined(GEOGRAPHICLIB_GEODESIC_ORDER)

#  define GEOGRAPHICLIB_GEODESIC_ORDER \
  (GEOGRAPHICLIB_PRECISION == 2 ? 6 : \
   (GEOGRAPHICLIB_PRECISION == 1 ? 3 : \
    (GEOGRAPHICLIB_PRECISION == 3 ? 7 : 8)))
#endif

namespace GeographicLib {

  class GeodesicLine;

  class GEOGRAPHICLIB_EXPORT Geodesic {
  private:
    typedef Math::real real;
    friend class GeodesicLine;
    static const int nA1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nC1_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nC1p_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nA2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nC2_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nA3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nA3x_ = nA3_;
    static const int nC3_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
    static const int nC4_ = GEOGRAPHICLIB_GEODESIC_ORDER;
    static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
    // Size for temporary array
    // nC = max(max(nC1_, nC1p_, nC2_) + 1, max(nC3_, nC4_))
    static const int nC_ = GEOGRAPHICLIB_GEODESIC_ORDER + 1;
    static const unsigned maxit1_ = 20;
    unsigned maxit2_;
    real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;

    enum captype {
      CAP_NONE = 0U,
      CAP_C1   = 1U<<0,
      CAP_C1p  = 1U<<1,
      CAP_C2   = 1U<<2,
      CAP_C3   = 1U<<3,
      CAP_C4   = 1U<<4,
      CAP_ALL  = 0x1FU,
      CAP_MASK = CAP_ALL,
      OUT_ALL  = 0x7F80U,
      OUT_MASK = 0xFF80U,       // Includes LONG_UNROLL
    };

    static real SinCosSeries(bool sinp,
                             real sinx, real cosx, const real c[], int n);
    static real Astroid(real x, real y);

    real _a, _f, _f1, _e2, _ep2, _n, _b, _c2, _etol2;
    real _A3x[nA3x_], _C3x[nC3x_], _C4x[nC4x_];

    void Lengths(real eps, real sig12,
                 real ssig1, real csig1, real dn1,
                 real ssig2, real csig2, real dn2,
                 real cbet1, real cbet2, unsigned outmask,
                 real& s12s, real& m12a, real& m0,
                 real& M12, real& M21, real Ca[]) const;
    real InverseStart(real sbet1, real cbet1, real dn1,
                      real sbet2, real cbet2, real dn2,
                      real lam12, real slam12, real clam12,
                      real& salp1, real& calp1,
                      real& salp2, real& calp2, real& dnm,
                      real Ca[]) const;
    real Lambda12(real sbet1, real cbet1, real dn1,
                  real sbet2, real cbet2, real dn2,
                  real salp1, real calp1, real slam120, real clam120,
                  real& salp2, real& calp2, real& sig12,
                  real& ssig1, real& csig1, real& ssig2, real& csig2,
                  real& eps, real& domg12,
                  bool diffp, real& dlam12, real Ca[]) const;
    real GenInverse(real lat1, real lon1, real lat2, real lon2,
                    unsigned outmask, real& s12,
                    real& salp1, real& calp1, real& salp2, real& calp2,
                    real& m12, real& M12, real& M21, real& S12) const;

    // These are Maxima generated functions to provide series approximations to
    // the integrals for the ellipsoidal geodesic.
    static real A1m1f(real eps);
    static void C1f(real eps, real c[]);
    static void C1pf(real eps, real c[]);
    static real A2m1f(real eps);
    static void C2f(real eps, real c[]);

    void A3coeff();
    real A3f(real eps) const;
    void C3coeff();
    void C3f(real eps, real c[]) const;
    void C4coeff();
    void C4f(real k2, real c[]) const;

  public:

    enum mask {
      NONE          = 0U,
      LATITUDE      = 1U<<7  | CAP_NONE,
      LONGITUDE     = 1U<<8  | CAP_C3,
      AZIMUTH       = 1U<<9  | CAP_NONE,
      DISTANCE      = 1U<<10 | CAP_C1,
      DISTANCE_IN   = 1U<<11 | CAP_C1 | CAP_C1p,
      REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
      GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
      AREA          = 1U<<14 | CAP_C4,
      LONG_UNROLL   = 1U<<15,
      ALL           = OUT_ALL| CAP_ALL,
    };


    Geodesic(real a, real f);


    Math::real Direct(real lat1, real lon1, real azi1, real s12,
                      real& lat2, real& lon2, real& azi2,
                      real& m12, real& M12, real& M21, real& S12)
      const {
      real t;
      return GenDirect(lat1, lon1, azi1, false, s12,
                       LATITUDE | LONGITUDE | AZIMUTH |
                       REDUCEDLENGTH | GEODESICSCALE | AREA,
                       lat2, lon2, azi2, t, m12, M12, M21, S12);
    }

    Math::real Direct(real lat1, real lon1, real azi1, real s12,
                      real& lat2, real& lon2)
      const {
      real t;
      return GenDirect(lat1, lon1, azi1, false, s12,
                       LATITUDE | LONGITUDE,
                       lat2, lon2, t, t, t, t, t, t);
    }

    Math::real Direct(real lat1, real lon1, real azi1, real s12,
                      real& lat2, real& lon2, real& azi2)
      const {
      real t;
      return GenDirect(lat1, lon1, azi1, false, s12,
                       LATITUDE | LONGITUDE | AZIMUTH,
                       lat2, lon2, azi2, t, t, t, t, t);
    }

    Math::real Direct(real lat1, real lon1, real azi1, real s12,
                      real& lat2, real& lon2, real& azi2, real& m12)
      const {
      real t;
      return GenDirect(lat1, lon1, azi1, false, s12,
                       LATITUDE | LONGITUDE | AZIMUTH | REDUCEDLENGTH,
                       lat2, lon2, azi2, t, m12, t, t, t);
    }

    Math::real Direct(real lat1, real lon1, real azi1, real s12,
                      real& lat2, real& lon2, real& azi2,
                      real& M12, real& M21)
      const {
      real t;
      return GenDirect(lat1, lon1, azi1, false, s12,
                       LATITUDE | LONGITUDE | AZIMUTH | GEODESICSCALE,
                       lat2, lon2, azi2, t, t, M12, M21, t);
    }

    Math::real Direct(real lat1, real lon1, real azi1, real s12,
                      real& lat2, real& lon2, real& azi2,
                      real& m12, real& M12, real& M21)
      const {
      real t;
      return GenDirect(lat1, lon1, azi1, false, s12,
                       LATITUDE | LONGITUDE | AZIMUTH |
                       REDUCEDLENGTH | GEODESICSCALE,
                       lat2, lon2, azi2, t, m12, M12, M21, t);
    }


    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2, real& azi2, real& s12,
                   real& m12, real& M12, real& M21, real& S12)
      const {
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
                REDUCEDLENGTH | GEODESICSCALE | AREA,
                lat2, lon2, azi2, s12, m12, M12, M21, S12);
    }

    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2) const {
      real t;
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE,
                lat2, lon2, t, t, t, t, t, t);
    }

    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2, real& azi2) const {
      real t;
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE | AZIMUTH,
                lat2, lon2, azi2, t, t, t, t, t);
    }

    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2, real& azi2, real& s12)
      const {
      real t;
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
                lat2, lon2, azi2, s12, t, t, t, t);
    }

    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2, real& azi2,
                   real& s12, real& m12) const {
      real t;
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
                REDUCEDLENGTH,
                lat2, lon2, azi2, s12, m12, t, t, t);
    }

    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2, real& azi2, real& s12,
                   real& M12, real& M21) const {
      real t;
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
                GEODESICSCALE,
                lat2, lon2, azi2, s12, t, M12, M21, t);
    }

    void ArcDirect(real lat1, real lon1, real azi1, real a12,
                   real& lat2, real& lon2, real& azi2, real& s12,
                   real& m12, real& M12, real& M21) const {
      real t;
      GenDirect(lat1, lon1, azi1, true, a12,
                LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
                REDUCEDLENGTH | GEODESICSCALE,
                lat2, lon2, azi2, s12, m12, M12, M21, t);
    }


    Math::real GenDirect(real lat1, real lon1, real azi1,
                         bool arcmode, real s12_a12, unsigned outmask,
                         real& lat2, real& lon2, real& azi2,
                         real& s12, real& m12, real& M12, real& M21,
                         real& S12) const;


    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& s12, real& azi1, real& azi2, real& m12,
                       real& M12, real& M21, real& S12) const {
      return GenInverse(lat1, lon1, lat2, lon2,
                        DISTANCE | AZIMUTH |
                        REDUCEDLENGTH | GEODESICSCALE | AREA,
                        s12, azi1, azi2, m12, M12, M21, S12);
    }

    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& s12) const {
      real t;
      return GenInverse(lat1, lon1, lat2, lon2,
                        DISTANCE,
                        s12, t, t, t, t, t, t);
    }

    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& azi1, real& azi2) const {
      real t;
      return GenInverse(lat1, lon1, lat2, lon2,
                        AZIMUTH,
                        t, azi1, azi2, t, t, t, t);
    }

    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& s12, real& azi1, real& azi2)
      const {
      real t;
      return GenInverse(lat1, lon1, lat2, lon2,
                        DISTANCE | AZIMUTH,
                        s12, azi1, azi2, t, t, t, t);
    }

    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& s12, real& azi1, real& azi2, real& m12)
      const {
      real t;
      return GenInverse(lat1, lon1, lat2, lon2,
                        DISTANCE | AZIMUTH | REDUCEDLENGTH,
                        s12, azi1, azi2, m12, t, t, t);
    }

    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& s12, real& azi1, real& azi2,
                       real& M12, real& M21) const {
      real t;
      return GenInverse(lat1, lon1, lat2, lon2,
                        DISTANCE | AZIMUTH | GEODESICSCALE,
                        s12, azi1, azi2, t, M12, M21, t);
    }

    Math::real Inverse(real lat1, real lon1, real lat2, real lon2,
                       real& s12, real& azi1, real& azi2, real& m12,
                       real& M12, real& M21) const {
      real t;
      return GenInverse(lat1, lon1, lat2, lon2,
                        DISTANCE | AZIMUTH |
                        REDUCEDLENGTH | GEODESICSCALE,
                        s12, azi1, azi2, m12, M12, M21, t);
    }


    Math::real GenInverse(real lat1, real lon1, real lat2, real lon2,
                          unsigned outmask,
                          real& s12, real& azi1, real& azi2,
                          real& m12, real& M12, real& M21, real& S12) const;


    GeodesicLine Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
      const;

    GeodesicLine InverseLine(real lat1, real lon1, real lat2, real lon2,
                             unsigned caps = ALL) const;

    GeodesicLine DirectLine(real lat1, real lon1, real azi1, real s12,
                            unsigned caps = ALL) const;

    GeodesicLine ArcDirectLine(real lat1, real lon1, real azi1, real a12,
                               unsigned caps = ALL) const;

    GeodesicLine GenDirectLine(real lat1, real lon1, real azi1,
                               bool arcmode, real s12_a12,
                               unsigned caps = ALL) const;


    Math::real MajorRadius() const { return _a; }

    Math::real Flattening() const { return _f; }

    Math::real EllipsoidArea() const
    { return 4 * Math::pi() * _c2; }

    static const Geodesic& WGS84();

  };

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_GEODESIC_HPP