GeodesicExact.hpp

Go to the documentation of this file.

#if !defined(GEOGRAPHICLIB_GEODESICEXACT_HPP)
#define GEOGRAPHICLIB_GEODESICEXACT_HPP 1

#include <GeographicLib/Constants.hpp>
#include <GeographicLib/EllipticFunction.hpp>

#if !defined(GEOGRAPHICLIB_GEODESICEXACT_ORDER)

#  define GEOGRAPHICLIB_GEODESICEXACT_ORDER 30
#endif

namespace GeographicLib {

  class GeodesicLineExact;

  class GEOGRAPHICLIB_EXPORT GeodesicExact {
  private:
    typedef Math::real real;
    friend class GeodesicLineExact;
    static const int nC4_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
    static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
    static const unsigned maxit1_ = 20;
    unsigned maxit2_;
    real tiny_, tol0_, tol1_, tol2_, tolb_, xthresh_;

    enum captype {
      CAP_NONE = 0U,
      CAP_E    = 1U<<0,
      // Skip 1U<<1 for compatibility with Geodesic (not required)
      CAP_D    = 1U<<2,
      CAP_H    = 1U<<3,
      CAP_C4   = 1U<<4,
      CAP_ALL  = 0x1FU,
      CAP_MASK = CAP_ALL,
      OUT_ALL  = 0x7F80U,
      OUT_MASK = 0xFF80U,       // Includes LONG_UNROLL
    };

    static real CosSeries(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 _C4x[nC4x_];

    void Lengths(const EllipticFunction& E,
                 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) const;
    real InverseStart(EllipticFunction& E,
                      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) 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,
                  EllipticFunction& E,
                  real& domg12, bool diffp, real& dlam12) 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 area.
    void C4coeff();
    void C4f(real k2, real c[]) const;
    // Large coefficients are split so that lo contains the low 52 bits and hi
    // the rest.  This choice avoids double rounding with doubles and higher
    // precision types.  float coefficients will suffer double rounding;
    // however the accuracy is already lousy for floats.
    static Math::real reale(long long hi, long long lo) {
      using std::ldexp;
      return ldexp(real(hi), 52) + lo;
    }

  public:

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


    GeodesicExact(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;


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

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

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

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

    GeodesicLineExact 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 GeodesicExact& WGS84();

  };

} // namespace GeographicLib

#endif  // GEOGRAPHICLIB_GEODESICEXACT_HPP