Gnomonic.cpp

/**
 * \file Gnomonic.cpp
 * \brief Implementation for GeographicLib::Gnomonic class
 *
 * Copyright (c) Charles Karney (2010-2019) <charles@karney.com> and licensed
 * under the MIT/X11 License.  For more information, see
 * https://geographiclib.sourceforge.io/
 **********************************************************************/

#include <GeographicLib/Gnomonic.hpp>

#if defined(_MSC_VER)
// Squelch warnings about potentially uninitialized local variables and
// constant conditional expressions
#  pragma warning (disable: 4701 4127)
#endif

namespace GeographicLib {

  using namespace std;

  Gnomonic::Gnomonic(const Geodesic& earth)
    : eps0_(numeric_limits<real>::epsilon())
    , eps_(real(0.01) * sqrt(eps0_))
    , _earth(earth)
    , _a(_earth.EquatorialRadius())
    , _f(_earth.Flattening())
  {}

  void Gnomonic::Forward(real lat0, real lon0, real lat, real lon,
                         real& x, real& y, real& azi, real& rk) const {
    real azi0, m, M, t;
    _earth.GenInverse(lat0, lon0, lat, lon,
                      Geodesic::AZIMUTH | Geodesic::REDUCEDLENGTH |
                      Geodesic::GEODESICSCALE,
                      t, azi0, azi, m, M, t, t);
    rk = M;
    if (M <= 0)
      x = y = Math::NaN();
    else {
      real rho = m/M;
      Math::sincosd(azi0, x, y);
      x *= rho; y *= rho;
    }
  }

  void Gnomonic::Reverse(real lat0, real lon0, real x, real y,
                         real& lat, real& lon, real& azi, real& rk) const {
    real
      azi0 = Math::atan2d(x, y),
      rho = Math::hypot(x, y),
      s = _a * atan(rho/_a);
    bool little = rho <= _a;
    if (!little)
      rho = 1/rho;
    GeodesicLine line(_earth.Line(lat0, lon0, azi0,
                                  Geodesic::LATITUDE | Geodesic::LONGITUDE |
                                  Geodesic::AZIMUTH | Geodesic::DISTANCE_IN |
                                  Geodesic::REDUCEDLENGTH |
                                  Geodesic::GEODESICSCALE));
    int count = numit_, trip = 0;
    real lat1, lon1, azi1, M;
    while (count-- || GEOGRAPHICLIB_PANIC) {
      real m, t;
      line.Position(s, lat1, lon1, azi1, m, M, t);
      if (trip)
        break;
      // If little, solve rho(s) = rho with drho(s)/ds = 1/M^2
      // else solve 1/rho(s) = 1/rho with d(1/rho(s))/ds = -1/m^2
      real ds = little ? (m - rho * M) * M : (rho * m - M) * m;
      s -= ds;
      // Reversed test to allow escape with NaNs
      if (!(abs(ds) >= eps_ * _a))
        ++trip;
    }
    if (trip) {
      lat = lat1; lon = lon1; azi = azi1; rk = M;
    } else
      lat = lon = azi = rk = Math::NaN();
    return;
  }

} // namespace GeographicLib