ADD: new track message, Entity class and Position class

libs/geographiclib/develop/CMakeLists.txt

@@ -0,0 +1,95 @@
+# Build the test programs...
+
+# Use fftw for AreaEst?  Because CMake support for installed fftw
+# package is broken, using FFTW required building local copies of the
+# library (a separate cmake configure and build for each precision).  So
+# default to not using fftw and fall back on kissfftw.  For these
+# reasons + the overall complexity of adding a dependency to
+# GeographicLib, the DST class uses kissfft (despite it being somewhat
+# slower that fftw).
+set (USE_FFTW OFF)
+if (USE_FFTW)
+  if (GEOGRAPHICLIB_PRECISION EQUAL 1)
+    set (FFTW_PACKAGE FFTW3f)
+  elseif (GEOGRAPHICLIB_PRECISION EQUAL 2)
+    set (FFTW_PACKAGE FFTW3)
+  elseif (GEOGRAPHICLIB_PRECISION EQUAL 3)
+    set (FFTW_PACKAGE FFTW3l)
+  elseif (GEOGRAPHICLIB_PRECISION EQUAL 4)
+    set (FFTW_PACKAGE FFTW3q)
+  else () # GEOGRAPHICLIB_PRECISION EQUAL 5
+    set (FFTW_PACKAGE OFF)
+  endif ()
+  if (FFTW_PACKAGE)
+    find_package (${FFTW_PACKAGE})
+    if (${FFTW_PACKAGE}_FOUND)
+      set (FFTW_FOUND ON)
+      set (FFTW_LIBRARIES ${${FFTW_PACKAGE}_LIBRARIES})
+      set (FFTW_LIBRARY_DIRS ${${FFTW_PACKAGE}_LIBRARY_DIRS})
+      set (FFTW_INCLUDE_DIRS ${${FFTW_PACKAGE}_INCLUDE_DIRS})
+    endif ()
+  endif ()
+endif ()
+
+set (DEVELPROGRAMS
+  ProjTest TMTest GeodTest ConicTest NaNTester HarmTest EllipticTest intersect
+  ClosestApproach M12zero GeodShort NormalTest)
+
+if (Boost_FOUND AND NOT GEOGRAPHICLIB_PRECISION EQUAL 4)
+  # Skip LevelEllipsoid for quad precision because of compiler errors
+  # with boost 1.69 and g++ 9.2.1 (Fedora 30).  Problem reported as
+  # https://github.com/boostorg/odeint/issues/40
+  set (DEVELPROGRAMS ${DEVELPROGRAMS} LevelEllipsoid)
+  include_directories ("${Boost_INCLUDE_DIRS}")
+  if (APPLE)
+    # Suppress warnings from Boost library
+    # warnings with Mac OS X and boost 1.63
+    # no warnings with Linux and boost 1.60
+    set (CMAKE_CXX_FLAGS
+      "${CMAKE_CXX_FLAGS} -Wno-unused-parameter -Wno-unused-variable -Wno-unused-local-typedef")
+  endif ()
+endif ()
+
+# Loop over all the tools, specifying the source and library.
+add_custom_target (develprograms)
+foreach (DEVELPROGRAM ${DEVELPROGRAMS})
+
+  add_executable (${DEVELPROGRAM} EXCLUDE_FROM_ALL ${DEVELPROGRAM}.cpp)
+  add_dependencies (develprograms ${DEVELPROGRAM})
+  target_link_libraries (${DEVELPROGRAM} ${PROJECT_LIBRARIES}
+    ${HIGHPREC_LIBRARIES})
+
+endforeach ()
+
+add_executable (GeodExact EXCLUDE_FROM_ALL GeodExact.cpp
+  Geodesic30.cpp GeodesicLine30.cpp
+  Geodesic30.hpp GeodesicLine30.hpp)
+add_dependencies (develprograms GeodExact)
+target_link_libraries (GeodExact ${PROJECT_LIBRARIES} ${HIGHPREC_LIBRARIES})
+set (DEVELPROGRAMS ${DEVELPROGRAMS} GeodExact)
+
+add_executable (AreaEst EXCLUDE_FROM_ALL AreaEst.cpp)
+add_dependencies (develprograms AreaEst)
+target_link_libraries (AreaEst ${PROJECT_LIBRARIES} ${FFTW_LIBRARIES}
+  ${HIGHPREC_LIBRARIES})
+if (FFTW_FOUND)
+  target_compile_definitions (AreaEst PUBLIC HAVE_FFTW=1)
+  target_include_directories (AreaEst PUBLIC ${FFTW_INCLUDE_DIRS})
+endif ()
+
+set (DEVELPROGRAMS ${DEVELPROGRAMS} AreaEst)
+
+add_executable (reformat EXCLUDE_FROM_ALL reformat.cpp)
+add_dependencies (develprograms reformat)
+set (DEVELPROGRAMS ${DEVELPROGRAMS} reformat)
+
+if (MSVC OR CMAKE_CONFIGURATION_TYPES)
+  # Add _d suffix for your debug versions of the tools
+  set_target_properties (${DEVELPROGRAMS} PROPERTIES
+    DEBUG_POSTFIX ${CMAKE_DEBUG_POSTFIX})
+endif ()
+
+# Put all the programs into a folder in the IDE
+set_property (TARGET develprograms ${DEVELPROGRAMS} PROPERTY FOLDER develop)
+
+# Don't install develop programs

libs/geographiclib/develop/ClosestApproach.cpp

@@ -0,0 +1,152 @@
+// Determine time and position of closest approach for two vessels traveling
+// along geodesics.  Thanks to
+//    Jorge D. Taramona
+//    Inspector de Navegación Aérea
+//    Dirección General de Aviación Civil del Perú
+// for suggesting this problem.  See discussion entry from 2014-08-19:
+// https://sourceforge.net/p/geographiclib/discussion/1026620/thread/33ce09e0/
+
+// Compile with
+// g++ -O3  -std=c++11 -o ClosestApproach ClosestApproach.cpp -lGeographic -L/usr/local/lib -Wl,-rpath=/usr/local/lib
+
+// Example: Planes leave
+
+// (1) Istanbul 42N 29E bearing 51W traveling at 900 km/hr
+// (2) Reyjkavik 64N 22W bearing 154E traveling at 800 km/hr
+
+// at the same time; compute time of closest approach.  Guess that the time of
+// closest approach is after 1 hr.  So run
+
+// echo 42 29 -51 900e3 64 -22 154 800e3 1 | ./ClosestApproach
+
+// 0 1 2696504 -3.532249e+12 46.74926 19.83775 57.4055 -16.17089
+// 1 2.751696 1083726 6.786464e+10 52.80517 -0.05892146 45.43817 -9.817806
+// 2 2.71894 1082700 -2657490 52.72565 0.3576737 45.66487 -9.90999
+// 3 2.718941 1082700 -0.005293798 52.72565 0.3576574 45.66486 -9.909987
+// 4 2.718941 1082700 0 52.72565 0.3576574 45.66486 -9.909987
+// 5 2.718941 1082700 0 52.72565 0.3576574 45.66486 -9.909987
+// 6 2.718941 1082700 0 52.72565 0.3576574 45.66486 -9.909987
+// 7 2.718941 1082700 0 52.72565 0.3576574 45.66486 -9.909987
+// 8 2.718941 1082700 0 52.72565 0.3576574 45.66486 -9.909987
+// 9 2.718941 1082700 0 52.72565 0.3576574 45.66486 -9.909987
+
+// At closest approach
+//   time = 2.72 hr
+//   dist = 1083 km
+//   position of first plane = 52.7N 0.3E
+//   position of second plane = 45.7N 9.9W
+
+// CAVEAT: if the initial guess is bad, this program may return the time at
+// which the distance is a maximum.
+
+// Explanation:
+
+// Consider f(t) = s12^2 / 2
+// Find zeros of g(t) = f'(t) = s12' * s12 by Newton's method
+// using g'(t) = (s12')^2 + s12'' * s12
+
+// Using s12^2 / 2 as the governing function (as opposed, for example, to
+// |s12|) means that in the planar limit f(t) is a quadratic function and g(t)
+// is a linear function.  So Newton's method just needs a single iteration.
+
+// Let gam{1,2} be angle between s12 and v{1,2}.
+
+// s12' = v2 * cos(gam2) - v1 * cos(gam1)
+// g(t) = s12 * (v2 * cos(gam2) - v1 * cos(gam1))
+// s12'' =  - v2*sin(gam2) * dgam2/dt + v1*sin(gam1) * dgam1/dt
+// dgam1/dt = (+ M12*v1*sin(gam1) - v2*sin(gam2)) / m12
+// dgam2/dt = (- M21*v2*sin(gam2) + v1*sin(gam1)) / m12
+// s12'' = (M12*v1^2*sin(gam1)^2 + M21*v2^2*sin(gam2)^2
+//          - 2*v1*v2*sin(gam1)*sin(gam2)) / m12
+// g'(t) = v1^2 * (cos(gam1)^2 + M12*s12/m12 * sin(gam1)^2)
+//       + v2^2 * (cos(gam2)^2 + M21*s12/m12 * sin(gam2)^2)
+//       - 2*v1*v2 * (cos(gam1)*cos(gam2) + s12/m12 * sin(gam1)*sin(gam2))
+
+// Planar limit (m12 = s12, M12 = M21 = 1):
+// g(t) = (v2-v1) . (r2-r1)
+// g'(t) = v1^2 + v2^2 - 2*v1*v2 * cos(gam1-gam2)
+//       = |v1-v2|^2 = const
+
+// Special case when v2 = 0 (simple interception), v1 = 1
+
+// s12' = - cos(gam1)
+// g(t) = - s12 * cos(gam1)
+// s12'' =  sin(gam1) * dgam1/dt
+// dgam1/dt = M12*sin(gam1) / m12
+// s12'' = M12*sin(gam1)^2 / m12
+// g'(t) =  (cos(gam1)^2 + M12*s12/m12 * sin(gam1)^2)
+
+// Planar limit (m12 = s12, M12 = M21 = 1):
+// g(t) = -v1 . (r2-r1)
+// g'(t) = 1 = const
+
+#include <iostream>
+#include <iomanip>
+#include <cmath>
+
+#include <GeographicLib/Geodesic.hpp>
+#include <GeographicLib/GeodesicExact.hpp>
+#include <GeographicLib/Utility.hpp>
+
+using namespace std;
+using namespace GeographicLib;
+
+int main() {
+  try {
+    typedef Math::real real;
+    Utility::set_digits();
+    const Geodesic/*Exact*/ g = Geodesic/*Exact*/(Constants::WGS84_a(),
+                                          Constants::WGS84_f());
+    real lat1, lon1, azi1, v1;
+    real lat2, lon2, azi2, v2;
+    real t0;
+    cout << setprecision(7);
+    while (cin
+           >> lat1 >> lon1 >> azi1 >> v1
+           >> lat2 >> lon2 >> azi2 >> v2
+           >> t0) {
+      real t = t0;
+      for (int i = 0; i < 10; ++i) {
+        // Compute positions at time t
+        real lat1a, lon1a, azi1a;
+        real lat2a, lon2a, azi2a;
+        g.Direct(lat1, lon1, azi1, t*v1, lat1a, lon1a, azi1a);
+        g.Direct(lat2, lon2, azi2, t*v2, lat2a, lon2a, azi2a);
+        // Compute distance at time t
+        real s12, azi1c, azi2c, m12, M12, M21;
+        g.Inverse(lat1a, lon1a, lat2a, lon2a,
+                  s12, azi1c, azi2c, m12, M12, M21);
+        real
+          r12 = s12 / m12,
+          gam1 = (azi1c - azi1a) * Math::degree(),
+          gam2 = (azi2c - azi2a) * Math::degree(),
+          x1 = v1 * cos(gam1), y1 = v1 * sin(gam1),
+          x2 = v2 * cos(gam2), y2 = v2 * sin(gam2);
+        // Find g = 0 using Newton's method where
+        // g = d( (1/2)*s12^2 ) / dt = s12 * d(s12)/dt
+        real g = s12 * (x2 - x1);
+        real dg =             // dg = dg/dt
+          (Math::sq(x1) + M12 * r12 * Math::sq(y1)) +
+          (Math::sq(x2) + M21 * r12 * Math::sq(y2)) -
+          2 * (x1 * x2 + r12 * y1 * y2);
+        cout << i << " "
+             << fixed << setprecision(12) << t << " " << s12 << " "
+             << defaultfloat << setprecision(6) << g << " "
+             << fixed << setprecision(20)
+             << lat1a << " " << lon1a << " "
+             << fixed << setprecision(6)
+             << lat2a << " " << lon2a << "\n";
+        t -= g/dg;              // Newton iteration
+      }
+    }
+  }
+  catch (const exception& e) {
+    cerr << "Caught exception: " << e.what() << "\n";
+    return 1;
+  }
+  catch (...) {
+    cerr << "Caught unknown exception\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/ConicTest.cpp

@@ -0,0 +1,273 @@
+/**
+ * \file ConicTest.cpp
+ * \brief Command line utility for testing transverse Mercator projections
+ *
+ * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ **********************************************************************/
+
+#include "GeographicLib/LambertConformalConic.hpp"
+#include "GeographicLib/AlbersEqualArea.hpp"
+#include "GeographicLib/Constants.hpp"
+#include "GeographicLib/Geodesic.hpp"
+#include "GeographicLib/DMS.hpp"
+#include <string>
+#include <limits>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <iomanip>
+#include <stdexcept>
+
+#if defined(_MSC_VER)
+// Squelch warnings about constant conditional expressions
+#  pragma warning (disable: 4127)
+#endif
+
+GeographicLib::Math::real
+dist(GeographicLib::Math::real a, GeographicLib::Math::real f,
+     GeographicLib::Math::real lat0, GeographicLib::Math::real lon0,
+     GeographicLib::Math::real lat1, GeographicLib::Math::real lon1) {
+  using namespace GeographicLib;
+  typedef Math::real real;
+  using std::cos; using std::sin; using std::sqrt; using std::hypot;
+  real
+    phi = lat0 * Math::degree(),
+    e2 = f * (2 - f),
+    sinphi = sin(phi),
+    n = 1/sqrt(1 - e2 * sinphi * sinphi),
+      // See Wikipedia article on latitude
+    hlon = cos(phi) * n,
+    hlat = (1 - e2) * n * n * n,
+    dlon = lon1 - lon0;
+  if (dlon >= 180) dlon -= 360;
+  else if (dlon < -180) dlon += 360;
+  return a * Math::degree() *
+    hypot((lat1 - lat0) * hlat, dlon * hlon);
+}
+
+class TestData {
+  // Read test data with one line of buffering
+public:
+  typedef GeographicLib::Math::real real;
+private:
+  std::istream& _is;
+  bool _usesaved;               // Should GetNext use saved values?
+  bool _valid;                  // Are there saved values?
+  real _lat0, _lat, _lon, _x, _y, _k;
+  TestData& operator=(const TestData&);
+public:
+  TestData(std::istream& is) : _is(is), _usesaved(false), _valid(false) {}
+  bool GetNext(real& lat0, real& lat, real& lon,
+               real& x, real& y, real& k) {
+    if (_usesaved)
+      _usesaved = false;
+    else {
+      // Avoid a warning about void* changed to bool
+      _valid = (_is >> _lat0 >> _lat >> _lon >> _x >> _y >> _k) ? true : false;
+      if (!_valid)
+        return false;
+    }
+    lat0 = _lat0; lat = _lat; lon = _lon; x = _x; y = _y; k = _k;
+    return true;
+  }
+  bool BackUp() {
+    if (!_valid || _usesaved)
+      return false;             // Can't backup up again
+    else {
+      _usesaved = true;         // Set flag for GetNext
+      return true;
+    }
+  }
+};
+
+int usage(int retval) {
+  ( retval ? std::cerr : std::cout ) <<
+"ConicTest -l -s\n\
+\n\
+Checks conic projections\n";
+  return retval;
+}
+
+int main(int argc, char* argv[]) {
+  using namespace GeographicLib;
+  using namespace std;
+  typedef Math::real real;
+  bool lambert = true;
+  bool albers = false;
+  bool checkstdlats = false;
+  real a = Constants::WGS84_a(), f = Constants::WGS84_f();
+  for (int m = 1; m < argc; ++m) {
+    string arg(argv[m]);
+    if (arg == "-l") {
+      lambert = true;
+      albers = false;
+    } else if (arg == "-a") {
+      lambert = false;
+      albers = true;
+    } else if (arg == "-s") {
+      checkstdlats = true;
+    } else if (arg == "-e") {
+      if (m + 2 >= argc) return usage(1);
+      try {
+        a = Utility::val<real>(string(argv[m + 1]));
+        f = Utility::fract<real>(string(argv[m + 2]));
+      }
+      catch (const exception& e) {
+        cerr << "Error decoding arguments of -e: " << e.what() << "\n";
+        return 1;
+      }
+      m += 2;
+      if (f > 1) f = 1/f;
+    } else
+      return usage(arg != "-h");
+  }
+
+  try {
+    if (checkstdlats) {         // stdin contains lat1 lat2 lat0 k0
+      cout << setprecision(17);
+      real quant = real(1e12);
+      while (true) {
+        real lat1, lat2, lat0, k0;
+        if (!(cin >> lat1 >> lat2 >> lat0 >> k0))
+          break;
+        int
+          sign1 = lat1 < 0 ? -1 : 1,
+          sign2 = lat2 < 0 ? -1 : 1;
+        lat1 = real(floor(sign1 * lat1 * quant + 0.5L));
+        lat2 = real(floor(sign2 * lat2 * quant + 0.5L));
+        real
+          colat1 = (90 * quant - lat1) / quant,
+          colat2 = (90 * quant - lat2) / quant;
+        lat1 /= quant;
+        lat2 /= quant;
+        real
+          sinlat1 = sign1 * (lat1 < 45 ? sin(lat1 * Math::degree())
+                             : cos(colat1 * Math::degree())),
+          sinlat2 = sign2 * (lat2 < 45 ? sin(lat2 * Math::degree())
+                             : cos(colat2 * Math::degree())),
+          coslat1 = (lat1 < 45 ? cos(lat1 * Math::degree())
+                     : sin(colat1 * Math::degree())),
+          coslat2 = (lat2 < 45 ? cos(lat2 * Math::degree())
+                     : sin(colat2 * Math::degree()));
+        lat1 *= sign1;
+        lat2 *= sign2;
+        const LambertConformalConic lam(a, f, /* real(lat1), real(lat2), */
+                                        real(sinlat1), real(coslat1),
+                                        real(sinlat2), real(coslat2),
+                                        real(1));
+        const AlbersEqualArea alb(a, f, /* real(lat1), real(lat2), */
+                                        real(sinlat1), real(coslat1),
+                                        real(sinlat2), real(coslat2),
+                                        real(1));
+        real
+          lat0a = albers ? alb.OriginLatitude() : lam.OriginLatitude();
+          //k0a = albers ? alb.CentralScale() : lam.CentralScale();
+        if (!(abs(lat0a-lat0) <= 4.5e-14))
+          cout << lat1 << " " << lat2 << " " << lat0
+               << " " << lat0a << " " << lat0a - lat0 << "\n";
+      }
+    } else { // Check projection
+      // stdin contains lat0 lat lon x y k
+      TestData testset(cin);
+      cout << setprecision(8);
+      while (true) {
+        real lat0, lat, lon, x, y, k;
+        if (!testset.GetNext(lat0, lat, lon, x, y, k))
+          break;
+        if (!testset.BackUp())
+          break;
+        int
+          sign0 = lat0 < 0 ? -1 : 1;
+        real quant = real(1e12);
+        real
+          lat00 = real(floor(sign0 * lat0 * quant + 0.5L)),
+          colat00 = (90 * quant - lat00) / quant;
+        lat00 /= quant;
+        real
+          sinlat0 = real(sign0 * (lat00 < 45 ?
+                                  sin(lat00 * Math::degree()) :
+                                  cos(colat00 * Math::degree()))),
+          coslat0 = real(lat00 < 45 ? cos(lat00 * Math::degree())
+                         : sin(colat00 * Math::degree()));
+        const LambertConformalConic lcc(a, f,
+                                        sinlat0, coslat0, sinlat0, coslat0,
+                                        real(1));
+        const AlbersEqualArea alb(a, f,
+                                  sinlat0, coslat0, sinlat0, coslat0,
+                                  real(1));
+        real maxerrf = 0, maxerrr = 0, maxerrkf = 0, maxerrkr = 0;
+        real latf = 0, lonf = 0, latr = 0, lonr = 0,
+          latkf = 0, lonkf = 0, latkr = 0, lonkr = 0;
+        // cout << "New lat0: " << lat0 << "\n";
+        while (true) {
+          real lat0x;
+          if (!testset.GetNext(lat0x, lat, lon, x, y, k))
+            break;
+          if (lat0 != lat0x) {
+            testset.BackUp();
+            break;
+          }
+          real latb, lonb, xa, ya, gammaa, gammab, ka, kb;
+          if (albers)
+            alb.Forward(real(0), real(lat), real(lon), xa, ya, gammaa, ka);
+          else
+            lcc.Forward(real(0), real(lat), real(lon), xa, ya, gammaa, ka);
+          real errf = real(hypot(real(xa) - x, real(ya) - y));
+          if (lambert)
+            errf /= real(k);
+          real errkf = real(abs(real(ka) - k)/k);
+          if (albers)
+            alb.Reverse(real(0), real(x), real(y), latb, lonb, gammab, kb);
+          else
+            lcc.Reverse(real(0), real(x), real(y), latb, lonb, gammab, kb);
+          real errr = real(dist(real(a), real(f),
+                                lat, lon, real(latb), real(lonb)));
+          /*
+          cout << latb << " " << lonb << " " << xa << " " << ya << " "
+                    << ka << " " << kb << " "
+                    << gammaa << " " << gammab << "\n";
+          */
+          real errkr = real(abs(real(kb) - k)/k);
+          if (!(errf <= maxerrf)) {
+            maxerrf = errf;
+            latf = real(lat);
+            lonf = real(lon);
+          }
+          if (!(errr <= maxerrr)) {
+            maxerrr = errr;
+            latr = real(lat);
+            lonr = real(lon);
+          }
+          if (!(errkf <= maxerrkf || abs(lat) >= 89)) {
+            maxerrkf = errkf;
+            latkf = real(lat);
+            lonkf = real(lon);
+          }
+          if (!(errkr <= maxerrkr || abs(lat) >= 89)) {
+            maxerrkr = errkr;
+            latkr = real(lat);
+            lonkr = real(lon);
+          }
+          cout << lat0 << " " << lat << " " << lon << " "
+               << errf << " " << errr << " "
+               << errkf << " " << errkr << "\n";
+        }
+        cout << "Max errf: " << lat0 << " "
+             << maxerrf << " " << latf << " " << lonf << "\n";
+        cout << "Max errr: " << lat0 << " "
+             << maxerrr << " " << latr << " " << lonr << "\n";
+        cout << "Max errkf: " << lat0 << " "
+             << maxerrkf << " " << latkf << " " << lonkf << "\n";
+        cout << "Max errkr: " << lat0 << " "
+             << maxerrkr << " " << latkr << " " << lonkr << "\n";
+      }
+    }
+  }
+  catch (const exception& e) {
+    cout << "ERROR: " << e.what() << "\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/EllipticTest.cpp

@@ -0,0 +1,163 @@
+#include <iostream>
+#include <iomanip>
+#include <GeographicLib/EllipticFunction.hpp>
+#include <GeographicLib/Ellipsoid.hpp>
+
+#if defined(_MSC_VER)
+// Squelch warnings about constant conditional expressions
+#  pragma warning (disable: 4127)
+#endif
+
+using namespace GeographicLib;
+using namespace std;
+int main() {
+  typedef GeographicLib::Math::real real;
+  try {
+    if (true) {
+      for (int i = -7; i <= 7; ++i) {
+        if (i == 0)
+          continue;
+        real a = 1, b = a * real(pow(real(2), i)), f = (a-b)/a;
+        Ellipsoid e(a, f);
+        real m = e.QuarterMeridian();
+        real na = a * real(1e7)/m, nb = na * (1-f);
+        Ellipsoid ne(na, f);
+        cout << b/a << " " << na << " " << nb << " "
+             << ne.QuarterMeridian() << "\n";
+      }
+      for (int i = -1; i <= 1; ++i) {
+        if (i == 0)
+          continue;
+        real a = 1, b = a * real(pow(real(3), i)), f = (a-b)/a;
+        Ellipsoid e(a, f);
+        real m = e.QuarterMeridian();
+        real na = a * real(1e7)/m, nb = na * (1-f);
+        Ellipsoid ne(na, f);
+        cout << b/a << " " << na << " " << nb << " "
+             << ne.QuarterMeridian() << "\n";
+      }
+      return 0;
+    }
+    if (false) {
+      // Longitude check
+      real a = real(1.5), b = 1, f = (a - b)/a,
+        e2 = (a*a - b*b)/(a*a), ep2 = (a*a - b*b)/(b*b),
+        alp0 = real(0.85), calp0 = cos(alp0), salp0 = sin(alp0),
+        k2 = ep2 * Math::sq(calp0), kp2 = k2/(1 + k2),
+        sigma = real(1.2), theta = atan(sqrt(1+k2) * tan(sigma)),
+        kap2 = Math::sq(calp0)/(1-Math::sq(salp0)*e2),
+        omg = atan(salp0 * tan(sigma)),
+        ups = atan(sqrt((1+ep2)/(1+k2*Math::sq(sin(sigma)))) * tan(omg)),
+        // psi = atan(sqrt((1+k2*Math::sq(sin(sigma)))/(1+ep2)) * tan(omg)),
+        psi = atan(sqrt((1-e2)/(1+k2*Math::sq(cos(theta)))) *
+                   salp0*tan(theta));
+      EllipticFunction ella1(-k2, Math::sq(calp0));
+      EllipticFunction ella2(-k2, -ep2);
+      EllipticFunction ellb1(kp2, kap2);
+      EllipticFunction ellb2(kp2, e2);
+      EllipticFunction ellc(kp2, kp2);
+      cout << setprecision(15);
+      cout << (1 - f) * salp0 * ella1.G(sigma) << " "
+           << ups - ep2/sqrt(1+ep2) * salp0 * ella2.H(sigma) << " "
+           << (1-f) * sqrt(1-kp2) * salp0 * ellb1.Pi(theta) << " "
+           << psi + (1-f) * sqrt(1-kp2) * salp0 *
+        (ellb2.F(theta) - ellb2.Pi(theta)) << "\n";
+      cout << b*ella1.E(sigma) << " "
+           << b * sqrt(1-kp2) * ellc.Pi(theta) << " "
+           << b / sqrt(1-kp2) *
+        (ellc.E(theta) - kp2 * cos(theta) * sin(theta)/
+         sqrt(1 - kp2*Math::sq(sin(theta)))) << "\n";
+      return 0;
+    }
+    if (false) {
+      real b = 1, a = 10,
+        e2 = (a*a - b*b)/(a*a), ep2 = (a*a - b*b)/(b*b);
+      EllipticFunction elle(e2, e2);
+      EllipticFunction ellep(-ep2, -ep2);
+      cout << fixed << setprecision(8);
+      for (int i = 0; i <= 90; i += 5) {
+        real
+          beta = i * Math::degree(),
+          phi = atan(sqrt(1 + ep2) * tan(beta)),
+          u = elle.F(phi),
+          y = b * ellep.E(beta),
+          M = a*elle.E();
+        cout << (y / M) * 90 << " "
+                  << i << " "
+                  << phi / Math::degree() << " "
+                  << (u / elle.K()) * 90 << "\n";
+      }
+      /* Create plot with
+t=load('data.txt');
+plot(t(:,1),t(:,1),'-k',...
+     t(:,1),t(:,2),'-kx',...
+     t(:,1),t(:,4),'-ko',...
+     t(:,1),t(:,3),'-k+')
+axis equal;axis([0,90,0,90]);
+title('meridian measures for a/b = 10');
+xlabel('meridian distance (\circ)');
+ylabel('measure (\circ)');
+legend('distance',...
+       'parametric latitude',...
+       'elliptic variable',...
+       'geographic latitude',...
+       'Location','SouthEast');
+set(gcf,'PaperSize',[4.5,4.5]);
+set(gcf,'PaperPosition',[0,0,4.5,4.5]);
+print meridian-measures.png -dpng
+
+       */
+      return 0;
+    }
+
+    if (false) {
+      real alpha2 = real(0.8), k2 = real(-0.4);
+      EllipticFunction ellG(k2,alpha2);
+      EllipticFunction ellH(k2,k2/alpha2);
+
+      cout << setprecision(10);
+      for (int i = -179; i <= 180; i += 10) {
+        real
+          phi = i * Math::degree(),
+          sn = sin(phi), cn = cos(phi), dn = ellG.Delta(sn, cn),
+          g = ellG.G(phi),
+          h = (k2/alpha2)*ellH.H(phi) + sqrt(1-k2/alpha2)/sqrt(1-alpha2)*
+          atan2(sqrt(1-alpha2)*sqrt(1-k2/alpha2)*sn, dn*cn);
+
+        cout << i << " " << g << " " << h << " " << h-g << "\n";
+      }
+      return 0;
+    }
+    if (false) {
+      // For tabulated values in A+S
+      real
+        ASalpha = 30*Math::degree<real>(),
+        k2 = Math::sq(sin(ASalpha)),
+        alpha2 = real(0.3);
+
+      EllipticFunction ell(k2, alpha2);
+      real dphi = Math::degree();
+      cout << fixed << setprecision(10);
+      for (int i = 0; i <= 90; i += 15) {
+        real phi = i * dphi;
+        cout << i << " "
+                  << ell.F(phi) << " "
+                  << ell.E(phi) << " "
+                  << ell.D(phi) << " "
+                  << ell.Pi(phi) << " "
+                  << ell.G(phi) << " "
+                  << ell.H(phi) << "\n";
+      }
+      return 0;
+    }
+  }
+  catch (const exception& e) {
+    cerr << "Caught exception: " << e.what() << "\n";
+    return 1;
+  }
+  catch (...) {
+    cerr << "Caught unknown exception\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/GeodExact.cpp

@@ -0,0 +1,34 @@
+#include <iostream>
+#include "Geodesic30.hpp"
+
+int main() {
+  {
+    GeographicLib::Geodesic30<double> g(6.4e6, 0.01);
+    {
+      double lat2, lon2, azi2;
+      g.Direct(10, 30, 20, 1e6, lat2, lon2, azi2);
+      std::cout << lat2 << " " << lon2 << " " << azi2 << "\n";
+    }
+    {
+      double s12, azi1, azi2;
+      g.Inverse(-20, 0, 10, 70, s12, azi1, azi2);
+      std::cout << s12 << " " << azi1 << " " << azi2 << "\n";
+    }
+  }
+  /*
+  {
+    GeographicLib::Geodesic30<long double> g(6.4e6L, 0.01L);
+    {
+      long double lat2, lon2, azi2;
+      g.Direct(10, 30, 20, 1e6L, lat2, lon2, azi2);
+      std::cout << lat2 << " " << lon2 << " " << azi2 << "\n";
+    }
+    {
+      long double s12, azi1, azi2;
+      g.Inverse(-20, 0, 10, 70, s12, azi1, azi2);
+      std::cout << s12 << " " << azi1 << " " << azi2 << "\n";
+    }
+  }
+  */
+  return 0;
+}

libs/geographiclib/develop/GeodShort.cpp

@@ -0,0 +1,337 @@
+/**
+ * \file GeodShort.cpp
+ *
+ * Test various solutions to the inverse problem in the limit of short lines
+ * (following Bowring 1981).
+ **********************************************************************/
+
+#include <ctime>                // for std::time
+#include <functional>           // for std::function and std::bind
+#include <random>               // for C++11 random numbers
+#include <limits>               // for digits
+#include <iostream>
+#include <iomanip>
+#include <algorithm>
+#include <GeographicLib/Geodesic.hpp>
+#include <GeographicLib/GeodesicExact.hpp>
+#include <GeographicLib/Utility.hpp>
+#include <GeographicLib/EllipticFunction.hpp>
+
+using namespace GeographicLib;
+using namespace std;
+
+class GeodShort {
+private:
+  typedef Math::real real;
+  real _a, _f, _f1, _b, _e2, _ep2, _e;
+  EllipticFunction _E;
+  inline real eatanhe(real x) const {
+    return _f >= 0 ? _e * atanh(_e * x) : - _e * atan(_e * x);
+  }
+  static inline real psi0f(real phi) { return asinh(tan(phi)); }
+  static inline real invpsi0f(real psi) { return atan(sinh(psi)); }
+  inline real psif(real phi) { return psi0f(phi) - eatanhe(sin(phi)); }
+  static inline void SinCosNorm(real& sinx, real& cosx) {
+    real r = hypot(sinx, cosx);
+    sinx /= r;
+    cosx /= r;
+  }
+  static real GreatCircle(real sbet1, real cbet1,
+                          real sbet2, real cbet2,
+                          real omg12, real& azi1, real& azi2) {
+    real
+      // bet12 = bet2 - bet1 in [0, pi)
+      sbet12 = sbet2 * cbet1 - cbet2 * sbet1,
+      sbet12a = sbet2 * cbet1 + cbet2 * sbet1,
+      somg12 = sin(omg12), comg12 = cos(omg12);
+    real
+      salp1 = cbet2 * somg12,
+      calp1 = (comg12 >= 0 ?
+               sbet12 + cbet2 * sbet1 * Math::sq(somg12) / (1 + comg12) :
+               sbet12a - cbet2 * sbet1 * Math::sq(somg12) / (1 - comg12));
+    real
+      ssig12 = hypot(salp1, calp1),
+      csig12 = sbet1 * sbet2 + cbet1 * cbet2 * comg12;
+
+    real
+      salp2 = cbet1 * somg12,
+      calp2 = sbet12 - cbet1 * sbet2 * Math::sq(somg12) / (1 + comg12),
+      sig12 = atan2(ssig12, csig12);
+
+    azi1 = atan2(salp1, calp1) / Math::degree();
+    azi2 = atan2(salp2, calp2) / Math::degree();
+    return sig12;
+  }
+public:
+  GeodShort(real a, real f)
+    : _a(a)
+    , _f(f)
+    , _f1(1 - _f)
+    , _b(_a * _f1)
+    , _e2(_f * (2 - _f))
+    , _ep2(_e2/Math::sq(_f1))
+    , _e(sqrt(abs(_e2)))
+    , _E(-_ep2) {}
+  real Inverse(real lat1, real lon1, real lat2, real lon2,
+               real& azi1, real& azi2) {
+    int mode = 1;
+    real
+      phi1 = Math::degree() * lat1,
+      phi2 = Math::degree() * lat2,
+      lam12 = Math::degree() * Math::AngNormalize(lon2 - lon1),
+      sbet1 = _f1 * sin(phi1), cbet1 = cos(phi1),
+      sbet2 = _f1 * sin(phi2), cbet2 = cos(phi2);
+    SinCosNorm(sbet1, cbet1); SinCosNorm(sbet2, cbet2);
+    real
+      dn1 = sqrt(1 + _ep2 * Math::sq(sbet1)),
+      dn2 = sqrt(1 + _ep2 * Math::sq(sbet2)),
+      dnm;
+    if (mode == 0)
+       dnm = (dn1 + dn2) / 2;
+    else {
+      real sbetm2 = Math::sq(sbet1 + sbet2);
+      sbetm2 = sbetm2 / (sbetm2 + Math::sq(cbet1 + cbet2));
+      dnm = sqrt(1 + _ep2 * sbetm2);
+    }
+    return _b * dnm * GreatCircle(sbet1, cbet1, sbet2, cbet2,
+                                  lam12 / (_f1 * dnm),
+                                  azi1, azi2);
+  }
+  real Inverse2(real lat1, real lon1, real lat2, real lon2,
+               real& azi1, real& azi2) {
+    real
+      phi1 = Math::degree() * lat1,
+      phi2 = Math::degree() * lat2,
+      lam12 = Math::degree() * Math::AngNormalize(lon2 - lon1),
+      sbet1 = _f1 * sin(phi1), cbet1 = cos(phi1),
+      sbet2 = _f1 * sin(phi2), cbet2 = cos(phi2);
+    SinCosNorm(sbet1, cbet1); SinCosNorm(sbet2, cbet2);
+    real dnm;
+    real sbetm2 = Math::sq(sbet1 + sbet2);
+    sbetm2 = sbetm2 / (sbetm2 + Math::sq(cbet1 + cbet2));
+    dnm = sqrt(1 + _ep2 * sbetm2);
+    // Adjust bet1 and bet2 via conformal map
+    real
+      A = 1/(dnm * _f1),
+      phim = atan((sbet1+sbet2)/(_f1*(cbet1+cbet2))),
+      betm = atan((sbet1+sbet2)/(cbet1+cbet2)),
+      psim = psif(phim),
+      psipm = psi0f(betm),
+      K = psipm - A * psim,
+      bet1 = invpsi0f( A*psif(phi1) + K ),
+      bet2 = invpsi0f( A*psif(phi2) + K );
+    sbet1 = sin(bet1); cbet1 = cos(bet1);
+    sbet2 = sin(bet2); cbet2 = cos(bet2);
+    return _b * dnm * GreatCircle(sbet1, cbet1, sbet2, cbet2,
+                                  lam12 / (_f1 * dnm),
+                                  azi1, azi2);
+  }
+  real Bowring0(real lat1, real lon1, real lat2, real lon2,
+                real& azi1, real& azi2) {
+    int mode = 2;
+    real
+      phi1 = Math::degree() * lat1,
+      phi2 = Math::degree() * lat2,
+      lam12 = Math::degree() * Math::AngNormalize(lon2 - lon1),
+      m = 0.5,
+      phim = ( abs(phi1) >= abs(phi2) ?
+               (1 - m) * phi1 + m * phi2 :
+               (1 - m) * phi2 + m * phi1 ),
+      dphi1 = phi1 - phim,
+      dphi2 = phi2 - phim,
+      cosm = cos(phim),
+      sinm = sin(phim),
+      A = sqrt(1 + _ep2 * Math::sq(Math::sq(cosm))),
+      B = sqrt(1 + _ep2 * Math::sq(cosm)),
+      C = sqrt(1 + _ep2),
+      omg12 = A * lam12,
+      phipm = atan(tan(phim)/B),
+      R = _a * C/(B*B),
+      phip1, phip2;
+    if (mode == 0) {
+      phip1 = phipm + (dphi1/B) * (1 + (3*_ep2*dphi1)/(4*B*B) *
+                                   sin(2 * (phim + dphi1/3)));
+      phip2 = phipm + (dphi2/B) * (1 + (3*_ep2*dphi2)/(4*B*B) *
+                                   sin(2 * (phim + dphi2/3)));
+    } else if (mode == 1) {
+      real
+        psi1 = psif(phi1), psi2 = psif(phi2),
+        psim = psif(phim), psipm = psi0f(phipm);
+      phip1 = invpsi0f(A * (psi1-psim) + psipm);
+      phip2 = invpsi0f(A * (psi2-psim) + psipm);
+    } else {
+      phip1 = phipm + (dphi1/B) *
+        ( 1 + _ep2*dphi1/(2*B*B) *
+          (3*cosm*sinm + dphi1*(1-2*sinm*sinm +
+                                _ep2 * (cosm*cosm * (1 + 3*sinm*sinm)))/B*B));
+      phip2 = phipm + (dphi2/B) *
+        ( 1 + _ep2*dphi2/(2*B*B) *
+          (3*cosm*sinm + dphi2*(1-2*sinm*sinm +
+                                _ep2 * (cosm*cosm * (1 + 3*sinm*sinm)))/B*B));
+    }
+    return R * GreatCircle(sin(phip1), cos(phip1), sin(phip2), cos(phip2),
+                           omg12, azi1, azi2);
+  }
+  real Bowring1(real lat1, real lon1, real lat2, real lon2,
+                real& azi1, real& azi2) {
+    real
+      phi1 = Math::degree() * lat1,
+      phi2 = Math::degree() * lat2,
+      lam12 = Math::degree() * Math::AngNormalize(lon2 - lon1),
+      bet1 = atan(_f1 * tan(phi1)),
+      bet2 = atan(_f1 * tan(phi2)),
+      betm = (bet1 + bet2)/2,
+      phim = atan(tan(betm) / _f1),
+      cosm = cos(phim),
+      A = sqrt(1 + _ep2 * Math::sq(Math::sq(cosm))),
+      B = sqrt(1 + _ep2 * Math::sq(cosm)),
+      C = sqrt(1 + _ep2),
+      omg12 = A * lam12,
+      phipm = atan(tan(phim)/B),
+      R = _a * C/(B*B),
+      m1 = _E.E(sin(bet1), cos(bet1), sqrt(1 + _ep2 * Math::sq(sin(bet1)))),
+      m2 = _E.E(sin(bet2), cos(bet2), sqrt(1 + _ep2 * Math::sq(sin(bet2)))),
+      mm = _E.E(sin(betm), cos(betm), sqrt(1 + _ep2 * Math::sq(sin(betm)))),
+      phip1 = phipm + (m1 - mm) * _a * _f1 / R,
+      phip2 = phipm + (m2 - mm) * _a * _f1 / R;
+    return R * GreatCircle(sin(phip1), cos(phip1), sin(phip2), cos(phip2),
+                           omg12, azi1, azi2);
+  }
+  real Bowring2(real lat1, real lon1, real lat2, real lon2,
+                real& azi1, real& azi2) {
+    real highfact = 1;
+    real
+      phi1 = Math::degree() * lat1,
+      phi2 = Math::degree() * lat2,
+      lam12 = Math::degree() * Math::AngNormalize(lon2 - lon1),
+      bet1 = atan(_f1 * tan(phi1)),
+      bet2 = atan(_f1 * tan(phi2)),
+      betm = (bet1 + bet2)/2,
+      sbetm = sin(betm),
+      cbetm = cos(betm),
+      phim = atan(tan(betm) / _f1),
+      cosm = cos(phim),
+      A = sqrt(1 + _ep2 * Math::sq(Math::sq(cosm))),
+      B = sqrt(1 + _ep2 * Math::sq(cosm)),
+      C = sqrt(1 + _ep2),
+      omg12 = A * lam12,
+      phipm = atan(tan(phim)/B),
+      R = _a * C/(B*B),
+      dbet1 = bet1-betm,
+      dbet2 = bet2-betm,
+      dnm2 = 1+_ep2*Math::sq(sbetm),
+      dnm = sqrt(dnm2),
+      phip1 = phipm + dbet1/dnm *
+      (1 + dbet1 * _ep2/(2*dnm2) *
+       ( cbetm*sbetm + highfact * dbet1 * ( Math::sq(cbetm) -
+                                            Math::sq(sbetm)*dnm2)/(3*dnm2) )),
+      phip2 = phipm + dbet2/dnm *
+      (1 + dbet2 * _ep2/(2*dnm2) *
+       ( cbetm*sbetm + highfact * dbet2 * ( Math::sq(cbetm) -
+                                            Math::sq(sbetm)*dnm2)/(3*dnm2) ));
+    return R * GreatCircle(sin(phip1), cos(phip1), sin(phip2), cos(phip2),
+                           omg12, azi1, azi2);
+  }
+};
+
+int main(int argc, char* argv[]) {
+  try {
+    // Args are f and sig
+    if (argc != 3) {
+      cerr << "Usage: GeodShort f sig\n";
+      return 1;
+    }
+    typedef Math::real real;
+    real
+      f = Utility::fract<real>(string(argv[1])),
+      sig = Utility::val<real>(string(argv[2]));
+    Geodesic g(1, f);
+    GeodesicExact ge(1, f);
+    GeodShort s(1, f);
+    bool exact = abs(f) > 0.02;
+    real norm, consist;
+    {
+      real m;
+      if (exact)
+        ge.Inverse(0, 0, 90, 0, m);
+      else
+        g.Inverse(0, 0, 90, 0, m);
+      norm = max(m, Math::pi()/2 * g.EquatorialRadius());
+      consist = min(m, Math::pi()/2 * g.EquatorialRadius()) /
+        (Math::pi()/2);
+    }
+    unsigned seed = random_device()(); // Set seed from random_device
+    mt19937 r(seed);                   // Initialize URNG
+    uniform_real_distribution<double> U;
+    cout << norm << " " << consist << " "
+         << f << " " << sig << " " << seed << endl;
+    real maxerr1 = -1, maxerr2 = -1, maxerr3 = -1;
+    for (unsigned k = 0; k < 10000000; ++k) {
+      real
+        lat1 = 90*real(U(r)),
+        lon1 = 0,
+        azi1 = 180*real(U(r)),
+        lat2, lon2, s12, azi2;
+      if (exact)
+        ge.ArcDirect(lat1, lon1, azi1, sig, lat2, lon2, azi2, s12);
+      else
+        g.ArcDirect(lat1, lon1, azi1, sig, lat2, lon2, azi2, s12);
+      real
+        s12a = s.Bowring2(lat1, lon1, lat2, lon2, azi1, azi2),
+        err1 = abs(s12a - s12) / norm;
+      if (err1 > maxerr1) {
+        maxerr1 = err1;
+        cout << "A " << k << " "
+                  << lat1 << " " << azi1 << " " << maxerr1 << endl;
+      }
+      real lat, lon, azi1a, azi2a;
+      if (exact)
+        ge.Direct(lat1, lon1, azi1, s12a, lat, lon);
+      else
+        g.Direct(lat1, lon1, azi1, s12a, lat, lon);
+      real err2 = s.Inverse(lat2, lon2, lat, lon, azi1a, azi2a);
+      if (exact)
+        ge.Direct(lat2, lon2, azi2, -s12a, lat, lon);
+      else
+        g.Direct(lat2, lon2, azi2, -s12a, lat, lon);
+      err2 = max(err2, s.Inverse(lat1, lon1, lat, lon, azi1a, azi2a)) / norm;
+      if (err2 > maxerr2) {
+        maxerr2 = err2;
+        cout << "B " << k << " "
+                  << lat1 << " " << azi1 << " " << maxerr2 << endl;
+      }
+      real latx, lonx;
+      if (exact) {
+        ge.Direct(lat1, lon1, azi1, s12a + consist, lat, lon);
+        ge.Direct(lat2, lon2, azi2, consist, latx, lonx);
+      } else {
+        g.Direct(lat1, lon1, azi1, s12a + consist, lat, lon);
+        g.Direct(lat2, lon2, azi2, consist, latx, lonx);
+      }
+      real err3 = s.Inverse(lat, lon, latx, lonx, azi1a, azi2a);
+      if (exact) {
+        ge.Direct(lat1, lon1, azi1, -consist, lat, lon);
+        ge.Direct(lat2, lon2, azi2, -s12a - consist, latx, lonx);
+      } else {
+        g.Direct(lat1, lon1, azi1, -consist, lat, lon);
+        g.Direct(lat2, lon2, azi2, -s12a - consist, latx, lonx);
+      }
+      err3 = max(err3, s.Inverse(lat, lon, latx, lonx, azi1a, azi2a)) / norm;
+      if (err3 > maxerr3) {
+        maxerr3 = err3;
+        cout << "C " << k << " "
+                  << lat1 << " " << azi1 << " " << maxerr3 << endl;
+      }
+    }
+  }
+  catch (const exception& e) {
+    cerr << "Caught exception: " << e.what() << "\n";
+    return 1;
+  }
+  catch (...) {
+    cerr << "Caught unknown exception\n";
+    return 1;
+  }
+  cout << "DONE\n";
+  return 0;
+}

libs/geographiclib/develop/GeodTest.cpp

@@ -0,0 +1,498 @@
+/**
+ * \file GeodTest.cpp
+ **********************************************************************/
+
+#include <string>
+#include <iostream>
+#include <iomanip>
+#include <sstream>
+
+#include "GeographicLib/Geodesic.hpp"
+#include "GeographicLib/GeodesicLine.hpp"
+#include "GeographicLib/GeodesicExact.hpp"
+#include "GeographicLib/GeodesicLineExact.hpp"
+#include "GeographicLib/Constants.hpp"
+#include <GeographicLib/Utility.hpp>
+
+#include <cmath>
+#include <vector>
+#include <utility>
+#include <algorithm>
+#include <limits>
+
+using namespace std;
+using namespace GeographicLib;
+
+int usage(int retval) {
+  ( retval ? cerr : cout ) <<
+"GeodTest [ -a | -c | -t0 | -t1 | -t2 | -t3 | -h ]\n\
+\n\
+Check GeographicLib::Geodesic class.\n\
+-a (default) accuracy test (reads test data on standard input)\n\
+-E accuracy test with GeodesicExact (reads test data on standard input)\n\
+-F accuracy test with GeodesicExact (reads test data on standard input\n\
+   first line gives a and f)\n\
+-c coverage test (reads test data on standard input)\n\
+-t0 time GeodecicLine with distances using synthetic data\n\
+-t1 time GeodecicLine with angles using synthetic data\n\
+-t2 time Geodecic::Direct using synthetic data\n\
+-t3 time Geodecic::Inverse with synthetic data\n\
+-T0 time GeodecicLineExact with distances using synthetic data\n\
+-T1 time GeodecicLineExact with angles using synthetic data\n\
+-T2 time GeodecicExact::Direct using synthetic data\n\
+-T3 time GeodecicExact::Inverse with synthetic data\n\
+\n\
+-c requires an instrumented version of Geodesic.\n";
+  return retval;
+}
+
+Math::real angdiff(Math::real a1, Math::real a2) {
+  Math::real d = a2 - a1;
+  if (d >= 180)
+    d -= 360;
+  else if (d < -180)
+    d += 360;
+  return d;
+}
+
+Math::real azidiff(Math::real lat,
+                   Math::real lon1, Math::real lon2,
+                   Math::real azi1, Math::real azi2) {
+  Math::real
+    phi = lat * Math::degree(),
+    alpha1 = azi1 * Math::degree(),
+    alpha2 = azi2 * Math::degree(),
+    dlam = angdiff(lon1, lon2) * Math::degree();
+  Math::real res = sin(alpha2-alpha1)*cos(dlam)
+    -cos(alpha2-alpha1)*sin(dlam)*sin(phi)
+    // -sin(alpha1)*cos(alpha2)*(1-cos(dlam))*cos(phi)*cos(phi)
+    ;
+  return res;
+}
+
+Math::real dist(Math::real a, Math::real f,
+                Math::real lat0, Math::real lon0,
+                Math::real lat1, Math::real lon1) {
+  //  typedef GeographicLib::Math::real real;
+  //  real s12;
+  //  GeographicLib::Geodesic::
+  //    WGS84.Inverse(real(lat0), real(lon0), real(lat1), real(lon1), s12);
+  //  return Math::real(s12);
+  a *= Math::degree();
+  if (abs(lat0 + lat1) > Math::real(179.998)) {
+    // Near pole, transform into polar coordinates
+    Math::real
+      r0 = 90 - abs(lat0),
+      r1 = 90 - abs(lat1),
+      lam0 = lon0 * Math::degree(),
+      lam1 = lon1 * Math::degree();
+    return (a / (1 - f)) *
+      hypot(r0 * cos(lam0) - r1 * cos(lam1), r0 * sin(lam0) - r1 * sin(lam1));
+  } else {
+    // Otherwise use cylindrical formula
+    Math::real
+      phi = lat0 * Math::degree(),
+      cphi = abs(lat0) <= 45 ? cos(phi)
+      : sin((90 - abs(lat0)) * Math::degree()),
+      e2 = f * (2 - f),
+      sinphi = sin(phi),
+      n = 1/sqrt(1 - e2 * sinphi * sinphi),
+      // See Wikipedia article on latitude
+      degreeLon = a * cphi * n,
+      degreeLat = a * (1 - e2) * n * n * n,
+      dlon = angdiff(lon1, lon0),
+      dlat = lat1 - lat0;
+    dlat *= degreeLat;
+    dlon *= degreeLon;
+    return hypot(dlat, dlon);
+  }
+}
+
+// err[0] error in position of point 2 for the direct problem.
+// err[1] error in azimuth at point 2 for the direct problem.
+// err[2] error in m12 for the direct problem & inverse (except near conjugacy)
+// err[3] error in s12 for the inverse problem.
+// err[4] error in the azimuths for the inverse problem scaled by m12.
+// err[5] consistency of the azimuths for the inverse problem.
+// err[6] area error direct & inverse (except near conjugacy)
+template<class test>
+void GeodError(const test& tgeod,
+               Math::real lat1, Math::real lon1, Math::real azi1,
+               Math::real lat2, Math::real lon2, Math::real azi2,
+               Math::real s12, Math::real /*a12*/,
+               Math::real m12, Math::real S12,
+               vector<Math::real>& err) {
+  Math::real tlat1, tlon1, tazi1, tlat2, tlon2, tazi2, ts12, tm12a, tm12b,
+    tM12, tM21, tS12a, tS12b /*, ta12*/;
+  Math::real rlat1, rlon1, razi1, rlat2, rlon2, razi2, rm12;
+  tgeod.Direct(lat1, lon1, azi1,  s12,
+               tlat2, tlon2, tazi2, tm12a,
+               tM12, tM21, tS12a);
+  tS12a -= tgeod.EllipsoidArea() * (tazi2-azi2)/720;
+  tgeod.Direct(lat2, lon2, azi2, -s12,
+               tlat1, tlon1, tazi1, tm12b,
+               tM12, tM21, tS12b);
+  tS12b -= tgeod.EllipsoidArea() * (tazi1-azi1)/720;
+  err[0] = max(dist(tgeod.EquatorialRadius(), tgeod.Flattening(),
+                    lat2, lon2, tlat2, tlon2),
+               dist(tgeod.EquatorialRadius(), tgeod.Flattening(),
+                    lat1, lon1, tlat1, tlon1));
+  err[1] = max(abs(azidiff(lat2, lon2, tlon2, azi2, tazi2)),
+               abs(azidiff(lat1, lon1, tlon1, azi1, tazi1))) *
+    tgeod.EquatorialRadius();
+  err[2] = max(abs(tm12a - m12), abs(tm12b + m12));
+  if (!isnan(S12))
+    err[6] = max(abs(tS12a - S12), abs(tS12b + S12)) / tgeod.EquatorialRadius();
+
+  /* ta12 = */ tgeod.Inverse(lat1, lon1, lat2, lon2,
+                             ts12, tazi1, tazi2, tm12a,
+                             tM12, tM21, tS12a);
+  tS12a -= tgeod.EllipsoidArea() * ((tazi2-azi2)-(tazi1-azi1))/720;
+  err[3] = abs(ts12 - s12);
+  err[4] = max(abs(angdiff(azi1, tazi1)), abs(angdiff(azi2, tazi2))) *
+    Math::degree() * abs(m12);
+  if (lat1 + lat2 == 0)
+    err[4] = min(err[4],
+                 max(abs(angdiff(azi1, tazi2)), abs(angdiff(azi2, tazi1))) *
+                 Math::degree() * abs(m12));
+  // m12 and S12 are very sensitive with the inverse problem near conjugacy
+  if (!(s12 > tgeod.EquatorialRadius() && m12 < 10e3)) {
+    err[2] = max(err[2], abs(tm12a - m12));
+    if (!isnan(S12))
+      err[6] = max(err[6], abs(tS12a - S12) / tgeod.EquatorialRadius());
+  }
+  if (s12 > tgeod.EquatorialRadius()) {
+    tgeod.Direct(lat1, lon1, tazi1,   ts12/2, rlat2, rlon2, razi2, rm12);
+    tgeod.Direct(lat2, lon2, tazi2, - ts12/2, rlat1, rlon1, razi1, rm12);
+    err[5] = dist(tgeod.EquatorialRadius(), tgeod.Flattening(),
+                  rlat1, rlon1, rlat2, rlon2);
+  } else {
+    tgeod.Direct(lat1, lon1, tazi1,
+                 ts12 + tgeod.EquatorialRadius(),
+                 rlat2, rlon2, razi2, rm12);
+    tgeod.Direct(lat2, lon2, tazi2, tgeod.EquatorialRadius(),
+                 rlat1, rlon1, razi1, rm12);
+    err[5] = dist(tgeod.EquatorialRadius(), tgeod.Flattening(),
+                  rlat1, rlon1, rlat2, rlon2);
+    tgeod.Direct(lat1, lon1, tazi1, - tgeod.EquatorialRadius(),
+                 rlat2, rlon2, razi2, rm12);
+    tgeod.Direct(lat2, lon2, tazi2,
+                 - ts12 - tgeod.EquatorialRadius(),
+                 rlat1, rlon1, razi1, rm12);
+    err[5] = max(err[5], dist(tgeod.EquatorialRadius(), tgeod.Flattening(),
+                              rlat1, rlon1, rlat2, rlon2));
+  }
+}
+
+int main(int argc, char* argv[]) {
+  Utility::set_digits();
+  Math::real a = Constants::WGS84_a();
+  Math::real f = Constants::WGS84_f();
+  bool timing = false;
+  int timecase = 0; // 0 = line, 1 = line ang, 2 = direct, 3 = inverse
+  bool accuracytest = true;
+  bool coverage = false;
+  bool exact = false;
+  if (argc == 2) {
+    string arg = argv[1];
+    if (arg == "-a") {
+      accuracytest = true;
+      coverage = false;
+      timing = false;
+      exact = false;
+    } else if (arg == "-E") {
+      accuracytest = true;
+      coverage = false;
+      timing = false;
+      exact = true;
+    } else if (arg == "-F") {
+      accuracytest = true;
+      coverage = false;
+      timing = false;
+      exact = true;
+      string s;
+      getline(cin, s);
+      istringstream str(s);
+      str >> a >> f;
+    } else if (arg == "-c") {
+      accuracytest = false;
+      coverage = true;
+      timing = false;
+      exact = false;
+    } else if (arg == "-t0") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 0;
+      exact = false;
+    } else if (arg == "-t1") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 1;
+      exact = false;
+    } else if (arg == "-t2") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 2;
+      exact = false;
+    } else if (arg == "-t3") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 3;
+      exact = false;
+    } else if (arg == "-T0") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 0;
+      exact = true;
+    } else if (arg == "-T1") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 1;
+      exact = true;
+    } else if (arg == "-T2") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 2;
+      exact = true;
+    } else if (arg == "-T3") {
+      accuracytest = false;
+      coverage = false;
+      timing = true;
+      timecase = 3;
+      exact = true;
+    } else
+      return usage(arg == "-h" ? 0 : 1);
+  } else if (argc > 2)
+    return usage(1);
+
+  if (timing) {
+    if (!exact) {
+      const Geodesic& geod = Geodesic::WGS84();
+      unsigned cnt = 0;
+      Math::real s = 0;
+      Math::real dl;
+      switch (timecase) {
+      case 0:
+        // Time Line
+        dl = 2e7/1000;
+        for (int i = 0; i <= 90; ++i) {
+          Math::real lat1 = i;
+          for (int j = 0; j <= 180; ++j) {
+            Math::real azi1 = j;
+            const GeodesicLine l(geod, lat1, 0.0, azi1);
+            for (int k = 0; k <= 1000; ++k) {
+              Math::real s12 = dl * k;
+              Math::real lat2, lon2;
+              l.Position(s12, lat2, lon2);
+              ++cnt;
+              s += lat2;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      case 1:
+        // Time Line ang
+        dl = Math::real(180)/1000;
+        for (int i = 0; i <= 90; ++i) {
+          Math::real lat1 = i;
+          for (int j = 0; j <= 180; ++j) {
+            Math::real azi1 = j;
+            GeodesicLine l(geod, lat1, 0.0, azi1);
+            for (int k = 0; k <= 1000; ++k) {
+              Math::real s12 = dl * k;
+              Math::real lat2, lon2;
+              l.ArcPosition(s12, lat2, lon2);
+              ++cnt;
+              s += lat2;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      case 2:
+        // Time Direct
+        dl = 2e7/200;
+        for (int i = 0; i <= 90; ++i) {
+          Math::real lat1 = i;
+          for (int j = 0; j <= 180; ++j) {
+            Math::real azi1 = j;
+            for (int k = 0; k <= 200; ++k) {
+              Math::real s12 = dl * k;
+              Math::real lat2, lon2;
+              geod.Direct(lat1, 0.0, azi1, s12, lat2, lon2);
+              ++cnt;
+              s += lat2;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      case 3:
+        // Time Inverse
+        for (int i = 1; i <= 179; i += 2) {
+          Math::real lat1 = i * Math::real(0.5);
+          for (int j = -179; j <= 179; j += 2) {
+            Math::real lat2 = j * Math::real(0.5);
+            for (int k = 1; k <= 359; k += 2) {
+              Math::real lon2 = k * Math::real(0.5);
+              Math::real s12;
+              geod.Inverse(lat1, 0.0, lat2, lon2, s12);
+              ++cnt;
+              s += s12;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      }
+    } else {
+      const GeodesicExact& geod = GeodesicExact::WGS84();
+      unsigned cnt = 0;
+      Math::real s = 0;
+      Math::real dl;
+      switch (timecase) {
+      case 0:
+        // Time Line
+        dl = 2e7/1000;
+        for (int i = 0; i <= 90; ++i) {
+          Math::real lat1 = i;
+          for (int j = 0; j <= 180; ++j) {
+            Math::real azi1 = j;
+            const GeodesicLineExact l(geod, lat1, 0.0, azi1);
+            for (int k = 0; k <= 1000; ++k) {
+              Math::real s12 = dl * k;
+              Math::real lat2, lon2;
+              l.Position(s12, lat2, lon2);
+              ++cnt;
+              s += lat2;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      case 1:
+        // Time Line ang
+        dl = Math::real(180)/1000;
+        for (int i = 0; i <= 90; ++i) {
+          Math::real lat1 = i;
+          for (int j = 0; j <= 180; ++j) {
+            Math::real azi1 = j;
+            GeodesicLineExact l(geod, lat1, 0.0, azi1);
+            for (int k = 0; k <= 1000; ++k) {
+              Math::real s12 = dl * k;
+              Math::real lat2, lon2;
+              l.ArcPosition(s12, lat2, lon2);
+              ++cnt;
+              s += lat2;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      case 2:
+        // Time Direct
+        dl = 2e7/200;
+        for (int i = 0; i <= 90; ++i) {
+          Math::real lat1 = i;
+          for (int j = 0; j <= 180; ++j) {
+            Math::real azi1 = j;
+            for (int k = 0; k <= 200; ++k) {
+              Math::real s12 = dl * k;
+              Math::real lat2, lon2;
+              geod.Direct(lat1, 0.0, azi1, s12, lat2, lon2);
+              ++cnt;
+              s += lat2;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      case 3:
+        // Time Inverse
+        for (int i = 1; i <= 179; i += 2) {
+          Math::real lat1 = i * Math::real(0.5);
+          for (int j = -179; j <= 179; j += 2) {
+            Math::real lat2 = j * Math::real(0.5);
+            for (int k = 1; k <= 359; k += 2) {
+              Math::real lon2 = k * Math::real(0.5);
+              Math::real s12;
+              geod.Inverse(lat1, 0.0, lat2, lon2, s12);
+              ++cnt;
+              s += s12;
+            }
+          }
+        }
+        cout << cnt << " " << s << "\n";
+        break;
+      }
+    }
+  }
+  else if (accuracytest || coverage) {
+    const Geodesic geod(a, f);
+    const GeodesicExact geode(a, f);
+
+    const unsigned NUMERR = 7;
+
+    cout << fixed << setprecision(2);
+    vector<Math::real> erra(NUMERR);
+    vector<Math::real> err(NUMERR, 0.0);
+    vector<unsigned> errind(NUMERR);
+    unsigned cnt = 0;
+    string s;
+    while (getline(cin, s)) {
+      istringstream str(s);
+      Math::real lat1l, lon1l, azi1l, lat2l, lon2l, azi2l,
+        s12l, a12l, m12l, S12l;
+      if (!(str >> lat1l >> lon1l >> azi1l
+                >> lat2l >> lon2l >> azi2l
+                >> s12l >> a12l >> m12l))
+        break;
+      if (!(str >> S12l))
+        S12l = Math::NaN();
+      if (coverage) {
+#if defined(GEOD_DIAG) && GEOD_DIAG
+        Math::real
+          lat1 = lat1l, lon1 = lon1l,
+          lat2 = lat2l, lon2 = lon2l,
+          azi1, azi2, s12, m12;
+        geod.Inverse(lat1, lon1, lat2, lon2, s12, azi1, azi2, m12);
+        cout << geod.coverage << " " << geod.niter << "\n";
+#endif
+      } else {
+        exact ?
+          GeodError< GeodesicExact >
+          (geode, lat1l, lon1l, azi1l,
+           lat2l, lon2l, azi2l,
+           s12l, a12l, m12l, S12l,
+           erra) :
+          GeodError< Geodesic >
+          (geod, lat1l, lon1l, azi1l,
+           lat2l, lon2l, azi2l,
+           s12l, a12l, m12l, S12l,
+           erra);
+        for (unsigned i = 0; i < NUMERR; ++i) {
+          if (isfinite(err[i]) && !(erra[i] <= err[i])) {
+            err[i] = erra[i];
+            errind[i] = cnt;
+          }
+        }
+        ++cnt;
+      }
+    }
+    if (accuracytest) {
+      Math::real mult = Math::real(Math::extra_digits() == 0 ? 1e9l :
+                                   Math::extra_digits() <= 3 ? 1e12l : 1e15l);
+      for (unsigned i = 0; i < NUMERR; ++i)
+        cout << i << " " << mult * err[i]
+             << " " << errind[i] << "\n";
+    }
+  }
+}

libs/geographiclib/develop/Geodesic30.hpp

@@ -0,0 +1,838 @@
+/**
+ * \file Geodesic30.hpp
+ * \brief Header for GeographicLib::Geodesic30 class
+ *
+ * Copyright (c) Charles Karney (2009-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ **********************************************************************/
+
+#if !defined(GEOGRAPHICLIB_GEODESICEXACT_HPP)
+#define GEOGRAPHICLIB_GEODESICEXACT_HPP 1
+
+#include <GeographicLib/Constants.hpp>
+
+#if !defined(GEOGRAPHICLIB_GEODESICEXACT_ORDER)
+/**
+ * The order of the expansions used by Geodesic30.
+ **********************************************************************/
+#  define GEOGRAPHICLIB_GEODESICEXACT_ORDER 30
+#endif
+
+namespace GeographicLib {
+
+  template<typename real> class GeodesicLine30;
+
+  /**
+   * \brief %Geodesic calculations
+   *
+   * The shortest path between two points on an ellipsoid at (\e lat1, \e lon1)
+   * and (\e lat2, \e lon2) is called the geodesic.  Its length is \e s12 and
+   * the geodesic from point 1 to point 2 has azimuths \e azi1 and \e azi2 at
+   * the two end points.  (The azimuth is the heading measured clockwise from
+   * north.  \e azi2 is the "forward" azimuth, i.e., the heading that takes you
+   * beyond point 2 not back to point 1.)
+   *
+   * Given \e lat1, \e lon1, \e azi1, and \e s12, we can determine \e lat2, \e
+   * lon2, and \e azi2.  This is the \e direct geodesic problem and its
+   * solution is given by the function Geodesic30::Direct.  (If \e s12 is
+   * sufficiently large that the geodesic wraps more than halfway around the
+   * earth, there will be another geodesic between the points with a smaller \e
+   * s12.)
+   *
+   * Given \e lat1, \e lon1, \e lat2, and \e lon2, we can determine \e azi1, \e
+   * azi2, and \e s12.  This is the \e inverse geodesic problem, whose solution
+   * is given by Geodesic30::Inverse.  Usually, the solution to the inverse
+   * problem is unique.  In cases where there are multiple solutions (all with
+   * the same \e s12, of course), all the solutions can be easily generated
+   * once a particular solution is provided.
+   *
+   * The standard way of specifying the direct problem is the specify the
+   * distance \e s12 to the second point.  However it is sometimes useful
+   * instead to specify the arc length \e a12 (in degrees) on the auxiliary
+   * sphere.  This is a mathematical construct used in solving the geodesic
+   * problems.  The solution of the direct problem in this form is provided by
+   * Geodesic30::ArcDirect.  An arc length in excess of 180&deg; indicates that
+   * the geodesic is not a shortest path.  In addition, the arc length between
+   * an equatorial crossing and the next extremum of latitude for a geodesic is
+   * 90&deg;.
+   *
+   * This class can also calculate several other quantities related to
+   * geodesics.  These are:
+   * - <i>reduced length</i>.  If we fix the first point and increase \e azi1
+   *   by \e dazi1 (radians), the second point is displaced \e m12 \e dazi1 in
+   *   the direction \e azi2 + 90&deg;.  The quantity \e m12 is called
+   *   the "reduced length" and is symmetric under interchange of the two
+   *   points.  On a curved surface the reduced length obeys a symmetry
+   *   relation, \e m12 + \e m21 = 0.  On a flat surface, we have \e m12 = \e
+   *   s12.  The ratio <i>s12</i>/\e m12 gives the azimuthal scale for an
+   *   azimuthal equidistant projection.
+   * - <i>geodesic scale</i>.  Consider a reference geodesic and a second
+   *   geodesic parallel to this one at point 1 and separated by a small
+   *   distance \e dt.  The separation of the two geodesics at point 2 is \e
+   *   M12 \e dt where \e M12 is called the "geodesic scale".  \e M21 is
+   *   defined similarly (with the geodesics being parallel at point 2).  On a
+   *   flat surface, we have \e M12 = \e M21 = 1.  The quantity 1/\e M12 gives
+   *   the scale of the Cassini-Soldner projection.
+   * - <i>area</i>.  Consider the quadrilateral bounded by the following lines:
+   *   the geodesic from point 1 to point 2, the meridian from point 2 to the
+   *   equator, the equator from \e lon2 to \e lon1, the meridian from the
+   *   equator to point 1.  The area of this quadrilateral is represented by \e
+   *   S12 with a clockwise traversal of the perimeter counting as a positive
+   *   area and it can be used to compute the area of any simple geodesic
+   *   polygon.
+   *
+   * Overloaded versions of Geodesic30::Direct, Geodesic30::ArcDirect,
+   * and Geodesic30::Inverse allow these quantities to be returned.  In
+   * addition there are general functions Geodesic30::GenDirect, and
+   * Geodesic30::GenInverse which allow an arbitrary set of results to be
+   * computed.  The quantities \e m12, \e M12, \e M21 which all specify the
+   * behavior of nearby geodesics obey addition rules.  Let points 1, 2, and 3
+   * all lie on a single geodesic, then
+   * - \e m13 = \e m12 \e M23 + \e m23 \e M21
+   * - \e M13 = \e M12 \e M23 &minus; (1 &minus; \e M12 \e M21) \e m23 / \e m12
+   * - \e M31 = \e M32 \e M21 &minus; (1 &minus; \e M23 \e M32) \e m12 / \e m23
+   *
+   * Additional functionality is provided by the GeodesicLine30 class, which
+   * allows a sequence of points along a geodesic to be computed.
+   *
+   * The calculations are accurate to better than 15 nm (15 nanometers).  See
+   * Sec. 9 of
+   * <a href="https://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a>
+   * for details.
+   *
+   * The algorithms are described in
+   * - C. F. F. Karney,
+   *   <a href="https://doi.org/10.1007/s00190-012-0578-z">
+   *   Algorithms for geodesics</a>,
+   *   J. Geodesy <b>87</b>, 43--55 (2013);
+   *   DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
+   *   10.1007/s00190-012-0578-z</a>;
+   *   addenda: <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
+   *   geod-addenda.html</a>.
+   * .
+   * For more information on geodesics see \ref geodesic.
+   **********************************************************************/
+
+  template<typename real>
+  class Geodesic30 {
+  private:
+    friend class GeodesicLine30<real>;
+    static const int nA1_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nC1_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nC1p_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nA2_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nC2_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nA3_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nA3x_ = nA3_;
+    static const int nC3_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nC3x_ = (nC3_ * (nC3_ - 1)) / 2;
+    static const int nC4_ = GEOGRAPHICLIB_GEODESICEXACT_ORDER;
+    static const int nC4x_ = (nC4_ * (nC4_ + 1)) / 2;
+    static const unsigned maxit_ = 50;
+
+    static const real tiny_;
+    static const real tol0_;
+    static const real tol1_;
+    static const real tol2_;
+    static const real 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,
+      OUT_ALL  = 0x7F80U,
+    };
+
+    static real SinCosSeries(bool sinp,
+                             real sinx, real cosx, const real c[], int n)
+     ;
+    static inline real AngRound(real x) {
+      // The makes the smallest gap in x = 1/16 - nextafter(1/16, 0) = 1/2^57
+      // for reals = 0.7 pm on the earth if x is an angle in degrees.  (This
+      // is about 1000 times more resolution than we get with angles around 90
+      // degrees.)  We use this to avoid having to deal with near singular
+      // cases when x is non-zero but tiny (e.g., 1.0e-200).
+      const real z = real(0.0625); // 1/16
+      volatile real y = std::abs(x);
+      // The compiler mustn't "simplify" z - (z - y) to y
+      y = y < z ? z - (z - y) : y;
+      return x < 0 ? -y : y;
+    }
+    static inline void SinCosNorm(real& sinx, real& cosx) {
+      using std::hypot;
+      real r = hypot(sinx, cosx);
+      sinx /= r;
+      cosx /= r;
+    }
+    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 ssig2, real csig2,
+                 real cbet1, real cbet2,
+                 real& s12s, real& m12a, real& m0,
+                 bool scalep, real& M12, real& M21,
+                 real C1a[], real C2a[]) const;
+    real InverseStart(real sbet1, real cbet1, real sbet2, real cbet2,
+                      real lam12,
+                      real& salp1, real& calp1,
+                      real& salp2, real& calp2,
+                      real C1a[], real C2a[]) const;
+    real Lambda12(real sbet1, real cbet1, real sbet2, real cbet2,
+                  real salp1, real calp1,
+                  real& salp2, real& calp2, real& sig12,
+                  real& ssig1, real& csig1, real& ssig2, real& csig2,
+                  real& eps, real& domg12, bool diffp, real& dlam12,
+                  real C1a[], real C2a[], real C3a[])
+      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:
+
+    /**
+     * Bit masks for what calculations to do.  These masks do double duty.
+     * They signify to the GeodesicLine30::GeodesicLine30 constructor and
+     * to Geodesic30::Line what capabilities should be included in the
+     * GeodesicLine30 object.  They also specify which results to return in
+     * the general routines Geodesic30::GenDirect and
+     * Geodesic30::GenInverse routines.  GeodesicLine30::mask is a
+     * duplication of this enum.
+     **********************************************************************/
+    enum mask {
+      /**
+       * No capabilities, no output.
+       * @hideinitializer
+       **********************************************************************/
+      NONE          = 0U,
+      /**
+       * Calculate latitude \e lat2.  (It's not necessary to include this as a
+       * capability to GeodesicLine30 because this is included by default.)
+       * @hideinitializer
+       **********************************************************************/
+      LATITUDE      = 1U<<7  | CAP_NONE,
+      /**
+       * Calculate longitude \e lon2.
+       * @hideinitializer
+       **********************************************************************/
+      LONGITUDE     = 1U<<8  | CAP_C3,
+      /**
+       * Calculate azimuths \e azi1 and \e azi2.  (It's not necessary to
+       * include this as a capability to GeodesicLine30 because this is
+       * included by default.)
+       * @hideinitializer
+       **********************************************************************/
+      AZIMUTH       = 1U<<9  | CAP_NONE,
+      /**
+       * Calculate distance \e s12.
+       * @hideinitializer
+       **********************************************************************/
+      DISTANCE      = 1U<<10 | CAP_C1,
+      /**
+       * Allow distance \e s12 to be used as input in the direct geodesic
+       * problem.
+       * @hideinitializer
+       **********************************************************************/
+      DISTANCE_IN   = 1U<<11 | CAP_C1 | CAP_C1p,
+      /**
+       * Calculate reduced length \e m12.
+       * @hideinitializer
+       **********************************************************************/
+      REDUCEDLENGTH = 1U<<12 | CAP_C1 | CAP_C2,
+      /**
+       * Calculate geodesic scales \e M12 and \e M21.
+       * @hideinitializer
+       **********************************************************************/
+      GEODESICSCALE = 1U<<13 | CAP_C1 | CAP_C2,
+      /**
+       * Calculate area \e S12.
+       * @hideinitializer
+       **********************************************************************/
+      AREA          = 1U<<14 | CAP_C4,
+      /**
+       * All capabilities, calculate everything.
+       * @hideinitializer
+       **********************************************************************/
+      ALL           = OUT_ALL| CAP_ALL,
+    };
+
+    /** \name Constructor
+     **********************************************************************/
+    ///@{
+    /**
+     * Constructor for an ellipsoid with
+     *
+     * @param[in] a equatorial radius (meters).
+     * @param[in] f flattening of ellipsoid.  Setting \e f = 0 gives a sphere.
+     *   Negative \e f gives a prolate ellipsoid.  If \e f > 1, set flattening
+     *   to 1/\e f.
+     * @exception GeographicErr if \e a or (1 &minus; \e f ) \e a is not
+     *   positive.
+     **********************************************************************/
+    Geodesic30(real a, real f);
+    ///@}
+
+    /** \name Direct geodesic problem specified in terms of distance.
+     **********************************************************************/
+    ///@{
+    /**
+     * Perform the direct geodesic calculation where the length of the geodesic
+     * is specified in terms of distance.
+     *
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] azi1 azimuth at point 1 (degrees).
+     * @param[in] s12 distance between point 1 and point 2 (meters); it can be
+     *   signed.
+     * @param[out] lat2 latitude of point 2 (degrees).
+     * @param[out] lon2 longitude of point 2 (degrees).
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] m12 reduced length of geodesic (meters).
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless).
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless).
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
+     * @return \e a12 arc length of between point 1 and point 2 (degrees).
+     *
+     * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
+     * azi1 should be in the range [&minus;540&deg;, 540&deg;).  The values of
+     * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
+     * 180&deg;).
+     *
+     * If either point is at a pole, the azimuth is defined by keeping the
+     * longitude fixed and writing \e lat = 90&deg; &minus; &epsilon; or
+     * &minus;90&deg; + &epsilon; and taking the limit &epsilon; &rarr; 0 from
+     * above.  An arc length greater that 180&deg; signifies a geodesic which
+     * is not a shortest path.  (For a prolate ellipsoid, an additional
+     * condition is necessary for a shortest path: the longitudinal extent must
+     * not exceed of 180&deg;.)
+     *
+     * The following functions are overloaded versions of Geodesic30::Direct
+     * which omit some of the output parameters.  Note, however, that the arc
+     * length is always computed and returned as the function value.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Direct.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Direct.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Direct.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Direct.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Direct.
+     **********************************************************************/
+    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);
+    }
+    ///@}
+
+    /** \name Direct geodesic problem specified in terms of arc length.
+     **********************************************************************/
+    ///@{
+    /**
+     * Perform the direct geodesic calculation where the length of the geodesic
+     * is specified in terms of arc length.
+     *
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] azi1 azimuth at point 1 (degrees).
+     * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
+     *   be signed.
+     * @param[out] lat2 latitude of point 2 (degrees).
+     * @param[out] lon2 longitude of point 2 (degrees).
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] s12 distance between point 1 and point 2 (meters).
+     * @param[out] m12 reduced length of geodesic (meters).
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless).
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless).
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
+     *
+     * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
+     * azi1 should be in the range [&minus;540&deg;, 540&deg;).  The values of
+     * \e lon2 and \e azi2 returned are in the range [&minus;180&deg;,
+     * 180&deg;).
+     *
+     * If either point is at a pole, the azimuth is defined by keeping the
+     * longitude fixed and writing \e lat = 90&deg; &minus; &epsilon; or
+     * &minus;90&deg; + &epsilon; and taking the limit &epsilon; &rarr; 0 from
+     * above.  An arc length greater that 180&deg; signifies a geodesic which
+     * is not a shortest path.  (For a prolate ellipsoid, an additional
+     * condition is necessary for a shortest path: the longitudinal extent must
+     * not exceed of 180&deg;.)
+     *
+     * The following functions are overloaded versions of Geodesic30::Direct
+     * which omit some of the output parameters.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::ArcDirect.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::ArcDirect.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::ArcDirect.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::ArcDirect.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::ArcDirect.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::ArcDirect.
+     **********************************************************************/
+    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);
+    }
+    ///@}
+
+    /** \name General version of the direct geodesic solution.
+     **********************************************************************/
+    ///@{
+
+    /**
+     * The general direct geodesic calculation.  Geodesic30::Direct and
+     * Geodesic30::ArcDirect are defined in terms of this function.
+     *
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] azi1 azimuth at point 1 (degrees).
+     * @param[in] arcmode boolean flag determining the meaning of the second
+     *   parameter.
+     * @param[in] s12_a12 if \e arcmode is false, this is the distance between
+     *   point 1 and point 2 (meters); otherwise it is the arc length between
+     *   point 1 and point 2 (degrees); it can be signed.
+     * @param[in] outmask a bitor'ed combination of Geodesic30::mask values
+     *   specifying which of the following parameters should be set.
+     * @param[out] lat2 latitude of point 2 (degrees).
+     * @param[out] lon2 longitude of point 2 (degrees).
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] s12 distance between point 1 and point 2 (meters).
+     * @param[out] m12 reduced length of geodesic (meters).
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless).
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless).
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
+     * @return \e a12 arc length of between point 1 and point 2 (degrees).
+     *
+     * The Geodesic30::mask values possible for \e outmask are
+     * - \e outmask |= Geodesic30::LATITUDE for the latitude \e lat2.
+     * - \e outmask |= Geodesic30::LONGITUDE for the latitude \e lon2.
+     * - \e outmask |= Geodesic30::AZIMUTH for the latitude \e azi2.
+     * - \e outmask |= Geodesic30::DISTANCE for the distance \e s12.
+     * - \e outmask |= Geodesic30::REDUCEDLENGTH for the reduced length \e
+     *   m12.
+     * - \e outmask |= Geodesic30::GEODESICSCALE for the geodesic scales \e
+     *   M12 and \e M21.
+     * - \e outmask |= Geodesic30::AREA for the area \e S12.
+     * .
+     * The function value \e a12 is always computed and returned and this
+     * equals \e s12_a12 is \e arcmode is true.  If \e outmask includes
+     * Geodesic30::DISTANCE and \e arcmode is false, then \e s12 = \e
+     * s12_a12.  It is not necessary to include Geodesic30::DISTANCE_IN in
+     * \e outmask; this is automatically included is \e arcmode is false.
+     **********************************************************************/
+    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;
+    ///@}
+
+    /** \name Inverse geodesic problem.
+     **********************************************************************/
+    ///@{
+    /**
+     * Perform the inverse geodesic calculation.
+     *
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] lat2 latitude of point 2 (degrees).
+     * @param[in] lon2 longitude of point 2 (degrees).
+     * @param[out] s12 distance between point 1 and point 2 (meters).
+     * @param[out] azi1 azimuth at point 1 (degrees).
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] m12 reduced length of geodesic (meters).
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless).
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless).
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
+     * @return \e a12 arc length of between point 1 and point 2 (degrees).
+     *
+     * \e lat1 and \e lat2 should be in the range [&minus;90&deg;, 90&deg;]; \e
+     * lon1 and \e lon2 should be in the range [&minus;540&deg;, 540&deg;).
+     * The values of \e azi1 and \e azi2 returned are in the range
+     * [&minus;180&deg;, 180&deg;).
+     *
+     * If either point is at a pole, the azimuth is defined by keeping the
+     * longitude fixed and writing \e lat = 90&deg; &minus; &epsilon; or
+     * &minus;90&deg; + &epsilon; and taking the limit &epsilon; &rarr; 0 from
+     * above.  If the routine fails to converge, then all the requested outputs
+     * are set to Math::NaN().  (Test for such results with Math::isnan.)  This
+     * is not expected to happen with ellipsoidal models of the earth; please
+     * report all cases where this occurs.
+     *
+     * The following functions are overloaded versions of Geodesic30::Inverse
+     * which omit some of the output parameters.  Note, however, that the arc
+     * length is always computed and returned as the function value.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Inverse.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Inverse.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Inverse.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Inverse.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Inverse.
+     **********************************************************************/
+    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);
+    }
+
+    /**
+     * See the documentation for Geodesic30::Inverse.
+     **********************************************************************/
+    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);
+    }
+    ///@}
+
+    /** \name General version of inverse geodesic solution.
+     **********************************************************************/
+    ///@{
+    /**
+     * The general inverse geodesic calculation.  Geodesic30::Inverse is
+     * defined in terms of this function.
+     *
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] lat2 latitude of point 2 (degrees).
+     * @param[in] lon2 longitude of point 2 (degrees).
+     * @param[in] outmask a bitor'ed combination of Geodesic30::mask values
+     *   specifying which of the following parameters should be set.
+     * @param[out] s12 distance between point 1 and point 2 (meters).
+     * @param[out] azi1 azimuth at point 1 (degrees).
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] m12 reduced length of geodesic (meters).
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless).
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless).
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>).
+     * @return \e a12 arc length of between point 1 and point 2 (degrees).
+     *
+     * The Geodesic30::mask values possible for \e outmask are
+     * - \e outmask |= Geodesic30::DISTANCE for the distance \e s12.
+     * - \e outmask |= Geodesic30::AZIMUTH for the latitude \e azi2.
+     * - \e outmask |= Geodesic30::REDUCEDLENGTH for the reduced length \e
+     *   m12.
+     * - \e outmask |= Geodesic30::GEODESICSCALE for the geodesic scales \e
+     *   M12 and \e M21.
+     * - \e outmask |= Geodesic30::AREA for the area \e S12.
+     * .
+     * The arc length is always computed and returned as the function value.
+     **********************************************************************/
+    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;
+    ///@}
+
+    /** \name Interface to GeodesicLine30.
+     **********************************************************************/
+    ///@{
+
+    /**
+     * Set up to compute several points on a single geodesic.
+     *
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] azi1 azimuth at point 1 (degrees).
+     * @param[in] caps bitor'ed combination of Geodesic30::mask values
+     *   specifying the capabilities the GeodesicLine30 object should
+     *   possess, i.e., which quantities can be returned in calls to
+     *   GeodesicLine::Position.
+     *
+     * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
+     * azi1 should be in the range [&minus;540&deg;, 540&deg;).
+     *
+     * The Geodesic30::mask values are
+     * - \e caps |= Geodesic30::LATITUDE for the latitude \e lat2; this is
+     *   added automatically
+     * - \e caps |= Geodesic30::LONGITUDE for the latitude \e lon2
+     * - \e caps |= Geodesic30::AZIMUTH for the latitude \e azi2; this is
+     *   added automatically
+     * - \e caps |= Geodesic30::DISTANCE for the distance \e s12
+     * - \e caps |= Geodesic30::REDUCEDLENGTH for the reduced length \e m12
+     * - \e caps |= Geodesic30::GEODESICSCALE for the geodesic scales \e M12
+     *   and \e M21
+     * - \e caps |= Geodesic30::AREA for the area \e S12
+     * - \e caps |= Geodesic30::DISTANCE_IN permits the length of the
+     *   geodesic to be given in terms of \e s12; without this capability the
+     *   length can only be specified in terms of arc length.
+     * .
+     * The default value of \e caps is Geodesic30::ALL which turns on all
+     * the capabilities.
+     *
+     * If the point is at a pole, the azimuth is defined by keeping the \e lon1
+     * fixed and writing \e lat1 = &plusmn;(90&deg; &minus; &epsilon;) and
+     * taking the limit &epsilon; &rarr; 0+.
+     **********************************************************************/
+    GeodesicLine30<real>
+      Line(real lat1, real lon1, real azi1, unsigned caps = ALL)
+      const;
+
+    ///@}
+
+    /** \name Inspector functions.
+     **********************************************************************/
+    ///@{
+
+    /**
+     * @return \e a the equatorial radius of the ellipsoid (meters).  This is
+     *   the value used in the constructor.
+     **********************************************************************/
+    real EquatorialRadius() const { return _a; }
+
+    /**
+     * @return \e f the  flattening of the ellipsoid.  This is the
+     *   value used in the constructor.
+     **********************************************************************/
+    real Flattening() const { return _f; }
+
+    /// \cond SKIP
+    /**
+     * <b>DEPRECATED</b>
+     * @return \e r the inverse flattening of the ellipsoid.
+     **********************************************************************/
+    real InverseFlattening() const { return 1/_f; }
+    /// \endcond
+
+    /**
+     * @return total area of ellipsoid in meters<sup>2</sup>.  The area of a
+     *   polygon encircling a pole can be found by adding
+     *   Geodesic30::EllipsoidArea()/2 to the sum of \e S12 for each side of
+     *   the polygon.
+     **********************************************************************/
+    real EllipsoidArea() const
+    { return 4 * Math::pi<real>() * _c2; }
+    ///@}
+
+    /**
+     * A global instantiation of Geodesic30 with the parameters for the WGS84
+     * ellipsoid.
+     **********************************************************************/
+    static const Geodesic30 WGS84;
+
+  };
+
+} // namespace GeographicLib
+
+#endif  // GEOGRAPHICLIB_GEODESICEXACT_HPP

libs/geographiclib/develop/GeodesicLine30.cpp

@@ -0,0 +1,270 @@
+/**
+ * \file GeodesicLine30.cpp
+ * \brief Implementation for GeographicLib::GeodesicLine30 class
+ *
+ * Copyright (c) Charles Karney (2009-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ *
+ * This is a reformulation of the geodesic problem.  The notation is as
+ * follows:
+ * - at a general point (no suffix or 1 or 2 as suffix)
+ *   - phi = latitude
+ *   - beta = latitude on auxiliary sphere
+ *   - omega = longitude on auxiliary sphere
+ *   - lambda = longitude
+ *   - alpha = azimuth of great circle
+ *   - sigma = arc length along great circle
+ *   - s = distance
+ *   - tau = scaled distance (= sigma at multiples of pi/2)
+ * - at northwards equator crossing
+ *   - beta = phi = 0
+ *   - omega = lambda = 0
+ *   - alpha = alpha0
+ *   - sigma = s = 0
+ * - a 12 suffix means a difference, e.g., s12 = s2 - s1.
+ * - s and c prefixes mean sin and cos
+ **********************************************************************/
+
+#include "GeodesicLine30.hpp"
+
+namespace GeographicLib {
+
+  using namespace std;
+
+  template<typename real>
+  GeodesicLine30<real>::GeodesicLine30(const Geodesic30<real>& g,
+                                       real lat1, real lon1, real azi1,
+                                       unsigned caps)
+    : _a(g._a)
+    , _f(g._f)
+    , _b(g._b)
+    , _c2(g._c2)
+    , _f1(g._f1)
+      // Always allow latitude and azimuth
+    , _caps(caps | LATITUDE | AZIMUTH)
+  {
+    azi1 = Math::AngNormalize(azi1);
+    // Guard against underflow in salp0
+    azi1 = Geodesic30<real>::AngRound(azi1);
+    lon1 = Math::AngNormalize(lon1);
+    _lat1 = lat1;
+    _lon1 = lon1;
+    _azi1 = azi1;
+    // alp1 is in [0, pi]
+    real alp1 = azi1 * Math::degree<real>();
+    // Enforce sin(pi) == 0 and cos(pi/2) == 0.  Better to face the ensuing
+    // problems directly than to skirt them.
+    _salp1 =     azi1  == -180 ? 0 : sin(alp1);
+    _calp1 = abs(azi1) ==   90 ? 0 : cos(alp1);
+    real cbet1, sbet1, phi;
+    phi = lat1 * Math::degree<real>();
+    // Ensure cbet1 = +epsilon at poles
+    sbet1 = _f1 * sin(phi);
+    cbet1 = abs(lat1) == 90 ? Geodesic30<real>::tiny_ : cos(phi);
+    Geodesic30<real>::SinCosNorm(sbet1, cbet1);
+
+    // Evaluate alp0 from sin(alp1) * cos(bet1) = sin(alp0),
+    _salp0 = _salp1 * cbet1; // alp0 in [0, pi/2 - |bet1|]
+    // Alt: calp0 = hypot(sbet1, calp1 * cbet1).  The following
+    // is slightly better (consider the case salp1 = 0).
+    _calp0 = hypot(_calp1, _salp1 * sbet1);
+    // Evaluate sig with tan(bet1) = tan(sig1) * cos(alp1).
+    // sig = 0 is nearest northward crossing of equator.
+    // With bet1 = 0, alp1 = pi/2, we have sig1 = 0 (equatorial line).
+    // With bet1 =  pi/2, alp1 = -pi, sig1 =  pi/2
+    // With bet1 = -pi/2, alp1 =  0 , sig1 = -pi/2
+    // Evaluate omg1 with tan(omg1) = sin(alp0) * tan(sig1).
+    // With alp0 in (0, pi/2], quadrants for sig and omg coincide.
+    // No atan2(0,0) ambiguity at poles since cbet1 = +epsilon.
+    // With alp0 = 0, omg1 = 0 for alp1 = 0, omg1 = pi for alp1 = pi.
+    _ssig1 = sbet1; _somg1 = _salp0 * sbet1;
+    _csig1 = _comg1 = sbet1 != 0 || _calp1 != 0 ? cbet1 * _calp1 : 1;
+    Geodesic30<real>::SinCosNorm(_ssig1, _csig1); // sig1 in (-pi, pi]
+    Geodesic30<real>::SinCosNorm(_somg1, _comg1);
+
+    _k2 = Math::sq(_calp0) * g._ep2;
+    real eps = _k2 / (2 * (1 + sqrt(1 + _k2)) + _k2);
+
+    if (_caps & CAP_C1) {
+      _A1m1 = Geodesic30<real>::A1m1f(eps);
+      Geodesic30<real>::C1f(eps, _C1a);
+      _B11 = Geodesic30<real>::SinCosSeries(true, _ssig1, _csig1, _C1a, nC1_);
+      real s = sin(_B11), c = cos(_B11);
+      // tau1 = sig1 + B11
+      _stau1 = _ssig1 * c + _csig1 * s;
+      _ctau1 = _csig1 * c - _ssig1 * s;
+      // Not necessary because C1pa reverts C1a
+      //    _B11 = -SinCosSeries(true, _stau1, _ctau1, _C1pa, nC1p_);
+    }
+
+    if (_caps & CAP_C1p)
+      Geodesic30<real>::C1pf(eps, _C1pa);
+
+    if (_caps & CAP_C2) {
+      _A2m1 = Geodesic30<real>::A2m1f(eps);
+      Geodesic30<real>::C2f(eps, _C2a);
+      _B21 = Geodesic30<real>::SinCosSeries(true, _ssig1, _csig1, _C2a, nC2_);
+    }
+
+    if (_caps & CAP_C3) {
+      g.C3f(eps, _C3a);
+      _A3c = -_f * _salp0 * g.A3f(eps);
+      _B31 = Geodesic30<real>::SinCosSeries(true, _ssig1, _csig1,
+                                            _C3a, nC3_-1);
+    }
+
+    if (_caps & CAP_C4) {
+      g.C4f(_k2, _C4a);
+      // Multiplier = a^2 * e^2 * cos(alpha0) * sin(alpha0)
+      _A4 = Math::sq(_a) * _calp0 * _salp0 * g._e2;
+      _B41 = Geodesic30<real>::SinCosSeries(false, _ssig1, _csig1, _C4a, nC4_);
+    }
+  }
+
+  template<typename real>
+  real GeodesicLine30<real>::GenPosition(bool arcmode, real s12_a12,
+                                            unsigned outmask,
+                                            real& lat2, real& lon2, real& azi2,
+                                            real& s12, real& m12,
+                                            real& M12, real& M21,
+                                            real& S12)
+  const {
+    outmask &= _caps & OUT_ALL;
+    if (!( Init() && (arcmode || (_caps & DISTANCE_IN & OUT_ALL)) ))
+      // Uninitialized or impossible distance calculation requested
+      return Math::NaN<real>();
+
+    // Avoid warning about uninitialized B12.
+    real sig12, ssig12, csig12, B12 = 0, AB1 = 0;
+    if (arcmode) {
+      // Interpret s12_a12 as spherical arc length
+      sig12 = s12_a12 * Math::degree<real>();
+      real s12a = abs(s12_a12);
+      s12a -= 180 * floor(s12a / 180);
+      ssig12 = s12a ==  0 ? 0 : sin(sig12);
+      csig12 = s12a == 90 ? 0 : cos(sig12);
+    } else {
+      // Interpret s12_a12 as distance
+      real
+        tau12 = s12_a12 / (_b * (1 + _A1m1)),
+        s = sin(tau12),
+        c = cos(tau12);
+      // tau2 = tau1 + tau12
+      B12 = - Geodesic30<real>::SinCosSeries(true, _stau1 * c + _ctau1 * s,
+                                     _ctau1 * c - _stau1 * s,
+                                     _C1pa, nC1p_);
+      sig12 = tau12 - (B12 - _B11);
+      ssig12 = sin(sig12);
+      csig12 = cos(sig12);
+    }
+
+    real omg12, lam12, lon12;
+    real ssig2, csig2, sbet2, cbet2, somg2, comg2, salp2, calp2;
+    // sig2 = sig1 + sig12
+    ssig2 = _ssig1 * csig12 + _csig1 * ssig12;
+    csig2 = _csig1 * csig12 - _ssig1 * ssig12;
+    if (outmask & (DISTANCE | REDUCEDLENGTH | GEODESICSCALE)) {
+      if (arcmode)
+        B12 = Geodesic30<real>::SinCosSeries(true, ssig2, csig2, _C1a, nC1_);
+      AB1 = (1 + _A1m1) * (B12 - _B11);
+    }
+    // sin(bet2) = cos(alp0) * sin(sig2)
+    sbet2 = _calp0 * ssig2;
+    // Alt: cbet2 = hypot(csig2, salp0 * ssig2);
+    cbet2 = hypot(_salp0, _calp0 * csig2);
+    if (cbet2 == 0)
+      // I.e., salp0 = 0, csig2 = 0.  Break the degeneracy in this case
+      cbet2 = csig2 = Geodesic30<real>::tiny_;
+    // tan(omg2) = sin(alp0) * tan(sig2)
+    somg2 = _salp0 * ssig2; comg2 = csig2;  // No need to normalize
+    // tan(alp0) = cos(sig2)*tan(alp2)
+    salp2 = _salp0; calp2 = _calp0 * csig2; // No need to normalize
+    // omg12 = omg2 - omg1
+    omg12 = atan2(somg2 * _comg1 - comg2 * _somg1,
+                  comg2 * _comg1 + somg2 * _somg1);
+
+    if (outmask & DISTANCE)
+      s12 = arcmode ? _b * ((1 + _A1m1) * sig12 + AB1) : s12_a12;
+
+    if (outmask & LONGITUDE) {
+      lam12 = omg12 + _A3c *
+        ( sig12 +
+          (Geodesic30<real>::SinCosSeries(true, ssig2, csig2, _C3a, nC3_-1)
+           - _B31));
+      lon12 = lam12 / Math::degree<real>();
+      lon12 = Math::AngNormalize(lon12);
+      lon2 = Math::AngNormalize(_lon1 + lon12);
+    }
+
+    if (outmask & LATITUDE)
+      lat2 = atan2(sbet2, _f1 * cbet2) / Math::degree<real>();
+
+    if (outmask & AZIMUTH)
+      // minus signs give range [-180, 180). 0- converts -0 to +0.
+      azi2 = 0 - atan2(-salp2, calp2) / Math::degree<real>();
+
+    if (outmask & (REDUCEDLENGTH | GEODESICSCALE)) {
+      real
+        ssig1sq = Math::sq(_ssig1),
+        ssig2sq = Math::sq( ssig2),
+        w1 = sqrt(1 + _k2 * ssig1sq),
+        w2 = sqrt(1 + _k2 * ssig2sq),
+        B22 = Geodesic30<real>::SinCosSeries(true, ssig2, csig2, _C2a, nC2_),
+        AB2 = (1 + _A2m1) * (B22 - _B21),
+        J12 = (_A1m1 - _A2m1) * sig12 + (AB1 - AB2);
+      if (outmask & REDUCEDLENGTH)
+        // Add parens around (_csig1 * ssig2) and (_ssig1 * csig2) to ensure
+        // accurate cancellation in the case of coincident points.
+        m12 = _b * ((w2 * (_csig1 * ssig2) - w1 * (_ssig1 * csig2))
+                  - _csig1 * csig2 * J12);
+      if (outmask & GEODESICSCALE) {
+        M12 = csig12 + (_k2 * (ssig2sq - ssig1sq) *  ssig2 / (w1 + w2)
+                        - csig2 * J12) * _ssig1 / w1;
+        M21 = csig12 - (_k2 * (ssig2sq - ssig1sq) * _ssig1 / (w1 + w2)
+                        - _csig1 * J12) * ssig2 / w2;
+      }
+    }
+
+    if (outmask & AREA) {
+      real
+        B42 = Geodesic30<real>::SinCosSeries(false, ssig2, csig2, _C4a, nC4_);
+      real salp12, calp12;
+      if (_calp0 == 0 || _salp0 == 0) {
+        // alp12 = alp2 - alp1, used in atan2 so no need to normalize
+        salp12 = salp2 * _calp1 - calp2 * _salp1;
+        calp12 = calp2 * _calp1 + salp2 * _salp1;
+        // The right thing appears to happen if alp1 = +/-180 and alp2 = 0, viz
+        // salp12 = -0 and alp12 = -180.  However this depends on the sign
+        // being attached to 0 correctly.  The following ensures the correct
+        // behavior.
+        if (salp12 == 0 && calp12 < 0) {
+          salp12 = Geodesic30<real>::tiny_ * _calp1;
+          calp12 = -1;
+        }
+      } else {
+        // tan(alp) = tan(alp0) * sec(sig)
+        // tan(alp2-alp1) = (tan(alp2) -tan(alp1)) / (tan(alp2)*tan(alp1)+1)
+        // = calp0 * salp0 * (csig1-csig2) / (salp0^2 + calp0^2 * csig1*csig2)
+        // If csig12 > 0, write
+        //   csig1 - csig2 = ssig12 * (csig1 * ssig12 / (1 + csig12) + ssig1)
+        // else
+        //   csig1 - csig2 = csig1 * (1 - csig12) + ssig12 * ssig1
+        // No need to normalize
+        salp12 = _calp0 * _salp0 *
+          (csig12 <= 0 ? _csig1 * (1 - csig12) + ssig12 * _ssig1 :
+           ssig12 * (_csig1 * ssig12 / (1 + csig12) + _ssig1));
+        calp12 = Math::sq(_salp0) + Math::sq(_calp0) * _csig1 * csig2;
+      }
+      S12 = _c2 * atan2(salp12, calp12) + _A4 * (B42 - _B41);
+    }
+
+    return arcmode ? s12_a12 : sig12 / Math::degree<real>();
+  }
+
+  template class GeodesicLine30<double>;
+#if GEOGRAPHICLIB_HAVE_LONG_DOUBLE
+  template class GeodesicLine30<long double>;
+#endif
+
+} // namespace GeographicLib

libs/geographiclib/develop/GeodesicLine30.hpp

@@ -0,0 +1,593 @@
+/**
+ * \file GeodesicLine30.hpp
+ * \brief Header for GeographicLib::GeodesicLine30 class
+ *
+ * Copyright (c) Charles Karney (2009-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ **********************************************************************/
+
+#if !defined(GEOGRAPHICLIB_GEODESICLINEEXACT_HPP)
+#define GEOGRAPHICLIB_GEODESICLINEEXACT_HPP 1
+
+#include <GeographicLib/Constants.hpp>
+#include "Geodesic30.hpp"
+
+namespace GeographicLib {
+
+  /**
+   * \brief A geodesic line
+   *
+   * GeodesicLine30 facilitates the determination of a series of points on a
+   * single geodesic.  The starting point (\e lat1, \e lon1) and the azimuth \e
+   * azi1 are specified in the constructor.  GeodesicLine30.Position returns
+   * the location of point 2 a distance \e s12 along the geodesic.
+   * Alternatively GeodesicLine30.ArcPosition gives the position of point 2
+   * an arc length \e a12 along the geodesic.
+   *
+   * The default copy constructor and assignment operators work with this
+   * class.  Similarly, a vector can be used to hold GeodesicLine30 objects.
+   *
+   * The calculations are accurate to better than 15 nm (15 nanometers).  See
+   * Sec. 9 of
+   * <a href="https://arxiv.org/abs/1102.1215v1">arXiv:1102.1215v1</a> for
+   * details.
+   *
+   * The algorithms are described in
+   * - C. F. F. Karney,
+   *   <a href="https://doi.org/10.1007/s00190-012-0578-z">
+   *   Algorithms for geodesics</a>,
+   *   J. Geodesy <b>87</b>, 43--55 (2013);
+   *   DOI: <a href="https://doi.org/10.1007/s00190-012-0578-z">
+   *   10.1007/s00190-012-0578-z</a>;
+   *   <a href="https://geographiclib.sourceforge.io/geod-addenda.html">
+   *   addenda</a>.
+   * .
+   * For more information on geodesics see \ref geodesic.
+   **********************************************************************/
+
+  template<typename real>
+  class GeodesicLine30 {
+  private:
+    friend class Geodesic30<real>;
+    static const int nC1_ = Geodesic30<real>::nC1_;
+    static const int nC1p_ = Geodesic30<real>::nC1p_;
+    static const int nC2_ = Geodesic30<real>::nC2_;
+    static const int nC3_ = Geodesic30<real>::nC3_;
+    static const int nC4_ = Geodesic30<real>::nC4_;
+
+    real _lat1, _lon1, _azi1;
+    real _a, _f, _b, _c2, _f1, _salp0, _calp0, _k2,
+      _salp1, _calp1, _ssig1, _csig1, _stau1, _ctau1, _somg1, _comg1,
+      _A1m1, _A2m1, _A3c, _B11, _B21, _B31, _A4, _B41;
+    // index zero elements of _C1a, _C1pa, _C2a, _C3a are unused
+    real _C1a[nC1_ + 1], _C1pa[nC1p_ + 1], _C2a[nC2_ + 1], _C3a[nC3_],
+      _C4a[nC4_];    // all the elements of _C4a are used
+    unsigned _caps;
+
+    enum captype {
+      CAP_NONE = Geodesic30<real>::CAP_NONE,
+      CAP_C1   = Geodesic30<real>::CAP_C1,
+      CAP_C1p  = Geodesic30<real>::CAP_C1p,
+      CAP_C2   = Geodesic30<real>::CAP_C2,
+      CAP_C3   = Geodesic30<real>::CAP_C3,
+      CAP_C4   = Geodesic30<real>::CAP_C4,
+      CAP_ALL  = Geodesic30<real>::CAP_ALL,
+      OUT_ALL  = Geodesic30<real>::OUT_ALL,
+    };
+  public:
+
+    /**
+     * Bit masks for what calculations to do.  They signify to the
+     * GeodesicLine30::GeodesicLine30 constructor and to
+     * Geodesic30::Line what capabilities should be included in the
+     * GeodesicLine30 object.  This is merely a duplication of
+     * Geodesic30::mask.
+     **********************************************************************/
+    enum mask {
+      /**
+       * No capabilities, no output.
+       * @hideinitializer
+       **********************************************************************/
+      NONE          = Geodesic30<real>::NONE,
+      /**
+       * Calculate latitude \e lat2.  (It's not necessary to include this as a
+       * capability to GeodesicLine30 because this is included by default.)
+       * @hideinitializer
+       **********************************************************************/
+      LATITUDE      = Geodesic30<real>::LATITUDE,
+      /**
+       * Calculate longitude \e lon2.
+       * @hideinitializer
+       **********************************************************************/
+      LONGITUDE     = Geodesic30<real>::LONGITUDE,
+      /**
+       * Calculate azimuths \e azi1 and \e azi2.  (It's not necessary to
+       * include this as a capability to GeodesicLine30 because this is
+       * included by default.)
+       * @hideinitializer
+       **********************************************************************/
+      AZIMUTH       = Geodesic30<real>::AZIMUTH,
+      /**
+       * Calculate distance \e s12.
+       * @hideinitializer
+       **********************************************************************/
+      DISTANCE      = Geodesic30<real>::DISTANCE,
+      /**
+       * Allow distance \e s12 to be used as input in the direct geodesic
+       * problem.
+       * @hideinitializer
+       **********************************************************************/
+      DISTANCE_IN   = Geodesic30<real>::DISTANCE_IN,
+      /**
+       * Calculate reduced length \e m12.
+       * @hideinitializer
+       **********************************************************************/
+      REDUCEDLENGTH = Geodesic30<real>::REDUCEDLENGTH,
+      /**
+       * Calculate geodesic scales \e M12 and \e M21.
+       * @hideinitializer
+       **********************************************************************/
+      GEODESICSCALE = Geodesic30<real>::GEODESICSCALE,
+      /**
+       * Calculate area \e S12.
+       * @hideinitializer
+       **********************************************************************/
+      AREA          = Geodesic30<real>::AREA,
+      /**
+       * All capabilities, calculate everything.
+       * @hideinitializer
+       **********************************************************************/
+      ALL           = Geodesic30<real>::ALL,
+    };
+
+    /** \name Constructors
+     **********************************************************************/
+    ///@{
+
+    /**
+     * Constructor for a geodesic line staring at latitude \e lat1, longitude
+     * \e lon1, and azimuth \e azi1 (all in degrees).
+     *
+     * @param[in] g A Geodesic30 object used to compute the necessary
+     *   information about the GeodesicLine30.
+     * @param[in] lat1 latitude of point 1 (degrees).
+     * @param[in] lon1 longitude of point 1 (degrees).
+     * @param[in] azi1 azimuth at point 1 (degrees).
+     * @param[in] caps bitor'ed combination of GeodesicLine30::mask values
+     *   specifying the capabilities the GeodesicLine30 object should
+     *   possess, i.e., which quantities can be returned in calls to
+     *   GeodesicLine::Position.
+     *
+     * \e lat1 should be in the range [&minus;90&deg;, 90&deg;]; \e lon1 and \e
+     * azi1 should be in the range [&minus;540&deg;, 540&deg;).
+     *
+     * The GeodesicLine30::mask values are
+     * - \e caps |= GeodesicLine30::LATITUDE for the latitude \e lat2; this
+     *   is added automatically
+     * - \e caps |= GeodesicLine30::LONGITUDE for the latitude \e lon2
+     * - \e caps |= GeodesicLine30::AZIMUTH for the latitude \e azi2; this is
+     *   added automatically
+     * - \e caps |= GeodesicLine30::DISTANCE for the distance \e s12
+     * - \e caps |= GeodesicLine30::REDUCEDLENGTH for the reduced length \e
+         m12
+     * - \e caps |= GeodesicLine30::GEODESICSCALE for the geodesic scales \e
+     *   M12 and \e M21
+     * - \e caps |= GeodesicLine30::AREA for the area \e S12
+     * - \e caps |= GeodesicLine30::DISTANCE_IN permits the length of the
+     *   geodesic to be given in terms of \e s12; without this capability the
+     *   length can only be specified in terms of arc length.
+     * .
+     * The default value of \e caps is GeodesicLine30::ALL which turns on
+     * all the capabilities.
+     *
+     * If the point is at a pole, the azimuth is defined by keeping the \e lon1
+     * fixed and writing \e lat1 = &plusmn;(90&deg; &minus; &epsilon;) and
+     * taking the limit &epsilon; &rarr; 0+.
+     **********************************************************************/
+    GeodesicLine30(const Geodesic30<real>& g, real lat1, real lon1, real azi1,
+                 unsigned caps = ALL)
+     ;
+
+    /**
+     * A default constructor.  If GeodesicLine30::Position is called on the
+     * resulting object, it returns immediately (without doing any
+     * calculations).  The object can be set with a call to
+     * Geodesic30::Line.  Use Init() to test whether object is still in this
+     * uninitialized state.
+     **********************************************************************/
+    GeodesicLine30() : _caps(0U) {}
+    ///@}
+
+    /** \name Position in terms of distance
+     **********************************************************************/
+    ///@{
+
+    /**
+     * Compute the position of point 2 which is a distance \e s12 (meters)
+     * from point 1.
+     *
+     * @param[in] s12 distance between point 1 and point 2 (meters); it can be
+     *   signed.
+     * @param[out] lat2 latitude of point 2 (degrees).
+     * @param[out] lon2 longitude of point 2 (degrees); requires that the
+     *   GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::LONGITUDE.
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] m12 reduced length of geodesic (meters); requires that the
+     *   GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::REDUCEDLENGTH.
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless); requires that the GeodesicLine30 object was
+     *   constructed with \e caps |= GeodesicLine30::GEODESICSCALE.
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless); requires that the GeodesicLine30 object was
+     *   constructed with \e caps |= GeodesicLine30::GEODESICSCALE.
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>); requires
+     *   that the GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::AREA.
+     * @return \e a12 arc length of between point 1 and point 2 (degrees).
+     *
+     * The values of \e lon2 and \e azi2 returned are in the range
+     * [&minus;180&deg;, 180&deg;).
+     *
+     * The GeodesicLine30 object \e must have been constructed with \e caps
+     * |= GeodesicLine30::DISTANCE_IN; otherwise Math::NaN() is returned and
+     * no parameters are set.  Requesting a value which the GeodesicLine30
+     * object is not capable of computing is not an error; the corresponding
+     * argument will not be altered.
+     *
+     * The following functions are overloaded versions of
+     * GeodesicLine30::Position which omit some of the output parameters.
+     * Note, however, that the arc length is always computed and returned as
+     * the function value.
+     **********************************************************************/
+    real Position(real s12,
+                        real& lat2, real& lon2, real& azi2,
+                        real& m12, real& M12, real& M21,
+                        real& S12) const {
+      real t;
+      return GenPosition(false, s12,
+                         LATITUDE | LONGITUDE | AZIMUTH |
+                         REDUCEDLENGTH | GEODESICSCALE | AREA,
+                         lat2, lon2, azi2, t, m12, M12, M21, S12);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::Position.
+     **********************************************************************/
+    real Position(real s12, real& lat2, real& lon2) const {
+      real t;
+      return GenPosition(false, s12,
+                         LATITUDE | LONGITUDE,
+                         lat2, lon2, t, t, t, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::Position.
+     **********************************************************************/
+    real Position(real s12, real& lat2, real& lon2,
+                        real& azi2) const {
+      real t;
+      return GenPosition(false, s12,
+                         LATITUDE | LONGITUDE | AZIMUTH,
+                         lat2, lon2, azi2, t, t, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::Position.
+     **********************************************************************/
+    real Position(real s12, real& lat2, real& lon2,
+                        real& azi2, real& m12) const {
+      real t;
+      return GenPosition(false, s12,
+                         LATITUDE | LONGITUDE |
+                         AZIMUTH | REDUCEDLENGTH,
+                         lat2, lon2, azi2, t, m12, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::Position.
+     **********************************************************************/
+    real Position(real s12, real& lat2, real& lon2,
+                        real& azi2, real& M12, real& M21)
+      const {
+      real t;
+      return GenPosition(false, s12,
+                         LATITUDE | LONGITUDE |
+                         AZIMUTH | GEODESICSCALE,
+                         lat2, lon2, azi2, t, t, M12, M21, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::Position.
+     **********************************************************************/
+    real Position(real s12,
+                        real& lat2, real& lon2, real& azi2,
+                        real& m12, real& M12, real& M21)
+      const {
+      real t;
+      return GenPosition(false, s12,
+                         LATITUDE | LONGITUDE | AZIMUTH |
+                         REDUCEDLENGTH | GEODESICSCALE,
+                         lat2, lon2, azi2, t, m12, M12, M21, t);
+    }
+
+    ///@}
+
+    /** \name Position in terms of arc length
+     **********************************************************************/
+    ///@{
+
+    /**
+     * Compute the position of point 2 which is an arc length \e a12 (degrees)
+     * from point 1.
+     *
+     * @param[in] a12 arc length between point 1 and point 2 (degrees); it can
+     *   be signed.
+     * @param[out] lat2 latitude of point 2 (degrees).
+     * @param[out] lon2 longitude of point 2 (degrees); requires that the
+     *   GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::LONGITUDE.
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] s12 distance between point 1 and point 2 (meters); requires
+     *   that the GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::DISTANCE.
+     * @param[out] m12 reduced length of geodesic (meters); requires that the
+     *   GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::REDUCEDLENGTH.
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless); requires that the GeodesicLine30 object was
+     *   constructed with \e caps |= GeodesicLine30::GEODESICSCALE.
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless); requires that the GeodesicLine30 object was
+     *   constructed with \e caps |= GeodesicLine30::GEODESICSCALE.
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>); requires
+     *   that the GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::AREA.
+     *
+     * The values of \e lon2 and \e azi2 returned are in the range
+     * [&minus;180&deg;, 180&deg;).
+     *
+     * Requesting a value which the GeodesicLine30 object is not capable of
+     * computing is not an error; the corresponding argument will not be
+     * altered.
+     *
+     * The following functions are overloaded versions of
+     * GeodesicLine30::ArcPosition which omit some of the output parameters.
+     **********************************************************************/
+    void ArcPosition(real a12, real& lat2, real& lon2, real& azi2,
+                     real& s12, real& m12, real& M12, real& M21,
+                     real& S12) const {
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE |
+                  REDUCEDLENGTH | GEODESICSCALE | AREA,
+                  lat2, lon2, azi2, s12, m12, M12, M21, S12);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::ArcPosition.
+     **********************************************************************/
+    void ArcPosition(real a12, real& lat2, real& lon2)
+      const {
+      real t;
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE,
+                  lat2, lon2, t, t, t, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::ArcPosition.
+     **********************************************************************/
+    void ArcPosition(real a12,
+                     real& lat2, real& lon2, real& azi2)
+      const {
+      real t;
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE | AZIMUTH,
+                  lat2, lon2, azi2, t, t, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::ArcPosition.
+     **********************************************************************/
+    void ArcPosition(real a12, real& lat2, real& lon2, real& azi2,
+                     real& s12) const {
+      real t;
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE | AZIMUTH | DISTANCE,
+                  lat2, lon2, azi2, s12, t, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::ArcPosition.
+     **********************************************************************/
+    void ArcPosition(real a12, real& lat2, real& lon2, real& azi2,
+                     real& s12, real& m12) const {
+      real t;
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE | AZIMUTH |
+                  DISTANCE | REDUCEDLENGTH,
+                  lat2, lon2, azi2, s12, m12, t, t, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::ArcPosition.
+     **********************************************************************/
+    void ArcPosition(real a12, real& lat2, real& lon2, real& azi2,
+                     real& s12, real& M12, real& M21)
+      const {
+      real t;
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE | AZIMUTH |
+                  DISTANCE | GEODESICSCALE,
+                  lat2, lon2, azi2, s12, t, M12, M21, t);
+    }
+
+    /**
+     * See the documentation for GeodesicLine30::ArcPosition.
+     **********************************************************************/
+    void ArcPosition(real a12, real& lat2, real& lon2, real& azi2,
+                     real& s12, real& m12, real& M12, real& M21)
+      const {
+      real t;
+      GenPosition(true, a12,
+                  LATITUDE | LONGITUDE | AZIMUTH |
+                  DISTANCE | REDUCEDLENGTH | GEODESICSCALE,
+                  lat2, lon2, azi2, s12, m12, M12, M21, t);
+    }
+    ///@}
+
+    /** \name The general position function.
+     **********************************************************************/
+    ///@{
+
+    /**
+     * The general position function.  GeodesicLine30::Position and
+     * GeodesicLine30::ArcPosition are defined in terms of this function.
+     *
+     * @param[in] arcmode boolean flag determining the meaning of the second
+     *   parameter; if arcmode is false, then the GeodesicLine30 object must
+     *   have been constructed with \e caps |= GeodesicLine30::DISTANCE_IN.
+     * @param[in] s12_a12 if \e arcmode is false, this is the distance between
+     *   point 1 and point 2 (meters); otherwise it is the arc length between
+     *   point 1 and point 2 (degrees); it can be signed.
+     * @param[in] outmask a bitor'ed combination of GeodesicLine30::mask
+     *   values specifying which of the following parameters should be set.
+     * @param[out] lat2 latitude of point 2 (degrees).
+     * @param[out] lon2 longitude of point 2 (degrees); requires that the
+     *   GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::LONGITUDE.
+     * @param[out] azi2 (forward) azimuth at point 2 (degrees).
+     * @param[out] s12 distance between point 1 and point 2 (meters); requires
+     *   that the GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::DISTANCE.
+     * @param[out] m12 reduced length of geodesic (meters); requires that the
+     *   GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::REDUCEDLENGTH.
+     * @param[out] M12 geodesic scale of point 2 relative to point 1
+     *   (dimensionless); requires that the GeodesicLine30 object was
+     *   constructed with \e caps |= GeodesicLine30::GEODESICSCALE.
+     * @param[out] M21 geodesic scale of point 1 relative to point 2
+     *   (dimensionless); requires that the GeodesicLine30 object was
+     *   constructed with \e caps |= GeodesicLine30::GEODESICSCALE.
+     * @param[out] S12 area under the geodesic (meters<sup>2</sup>); requires
+     *   that the GeodesicLine30 object was constructed with \e caps |=
+     *   GeodesicLine30::AREA.
+     * @return \e a12 arc length of between point 1 and point 2 (degrees).
+     *
+     * The GeodesicLine30::mask values possible for \e outmask are
+     * - \e outmask |= GeodesicLine30::LATITUDE for the latitude \e lat2.
+     * - \e outmask |= GeodesicLine30::LONGITUDE for the latitude \e lon2.
+     * - \e outmask |= GeodesicLine30::AZIMUTH for the latitude \e azi2.
+     * - \e outmask |= GeodesicLine30::DISTANCE for the distance \e s12.
+     * - \e outmask |= GeodesicLine30::REDUCEDLENGTH for the reduced length
+     *   \e m12.
+     * - \e outmask |= GeodesicLine30::GEODESICSCALE for the geodesic scales
+     *   \e M12 and \e M21.
+     * - \e outmask |= GeodesicLine30::AREA for the area \e S12.
+     * .
+     * Requesting a value which the GeodesicLine30 object is not capable of
+     * computing is not an error; the corresponding argument will not be
+     * altered.  Note, however, that the arc length is always computed and
+     * returned as the function value.
+     **********************************************************************/
+    real GenPosition(bool arcmode, real s12_a12, unsigned outmask,
+                           real& lat2, real& lon2, real& azi2,
+                           real& s12, real& m12, real& M12, real& M21,
+                           real& S12) const;
+
+    ///@}
+
+    /** \name Inspector functions
+     **********************************************************************/
+    ///@{
+
+    /**
+     * @return true if the object has been initialized.
+     **********************************************************************/
+    bool Init() const { return _caps != 0U; }
+
+    /**
+     * @return \e lat1 the latitude of point 1 (degrees).
+     **********************************************************************/
+    real Latitude() const
+    { return Init() ? _lat1 : Math::NaN<real>(); }
+
+    /**
+     * @return \e lon1 the longitude of point 1 (degrees).
+     **********************************************************************/
+    real Longitude() const
+    { return Init() ? _lon1 : Math::NaN<real>(); }
+
+    /**
+     * @return \e azi1 the azimuth (degrees) of the geodesic line at point 1.
+     **********************************************************************/
+    real Azimuth() const
+    { return Init() ? _azi1 : Math::NaN<real>(); }
+
+    /**
+     * @return \e azi0 the azimuth (degrees) of the geodesic line as it crosses
+     * the equator in a northward direction.
+     **********************************************************************/
+    real EquatorialAzimuth() const {
+      using std::atan2;
+      return Init() ?
+        atan2(_salp0, _calp0) / Math::degree<real>() : Math::NaN<real>();
+    }
+
+    /**
+     * @return \e a1 the arc length (degrees) between the northward equatorial
+     * crossing and point 1.
+     **********************************************************************/
+    real EquatorialArc() const {
+      using std::atan2;
+      return Init() ?
+        atan2(_ssig1, _csig1) / Math::degree<real>() : Math::NaN<real>();
+    }
+
+    /**
+     * @return \e a the equatorial radius of the ellipsoid (meters).  This is
+     *   the value inherited from the Geodesic30 object used in the
+     *   constructor.
+     **********************************************************************/
+    real EquatorialRadius() const
+    { return Init() ? _a : Math::NaN<real>(); }
+
+    /**
+     * @return \e f the flattening of the ellipsoid.  This is the value
+     *   inherited from the Geodesic30 object used in the constructor.
+     **********************************************************************/
+    real Flattening() const
+    { return Init() ? _f : Math::NaN<real>(); }
+
+    /// \cond SKIP
+    /**
+     * <b>DEPRECATED</b>
+     * @return \e r the inverse flattening of the ellipsoid.
+     **********************************************************************/
+    real InverseFlattening() const
+    { return Init() ? 1/_f : Math::NaN<real>(); }
+    /// \endcond
+
+    /**
+     * @return \e caps the computational capabilities that this object was
+     *   constructed with.  LATITUDE and AZIMUTH are always included.
+     **********************************************************************/
+    unsigned Capabilities() const { return _caps; }
+
+    /**
+     * @param[in] testcaps a set of bitor'ed GeodesicLine30::mask values.
+     * @return true if the GeodesicLine30 object has all these capabilities.
+     **********************************************************************/
+    bool Capabilities(unsigned testcaps) const {
+      testcaps &= OUT_ALL;
+      return (_caps & testcaps) == testcaps;
+    }
+    ///@}
+
+  };
+
+} // namespace GeographicLib
+
+#endif  // GEOGRAPHICLIB_GEODESICLINEEXACT_HPP

libs/geographiclib/develop/HarmTest.cpp

@@ -0,0 +1,240 @@
+#include <fstream>
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <ctime>
+#include <GeographicLib/MagneticModel.hpp>
+#include <GeographicLib/SphericalHarmonic.hpp>
+#include <GeographicLib/CircularEngine.hpp>
+#include <GeographicLib/MagneticCircle.hpp>
+#include <GeographicLib/NormalGravity.hpp>
+#include <GeographicLib/GravityModel.hpp>
+#include <GeographicLib/GravityCircle.hpp>
+#include <GeographicLib/Utility.hpp>
+#include <GeographicLib/Geoid.hpp>
+
+#if defined(_MSC_VER)
+// Squelch warnings about constant conditional and enum-float expressions,
+// potentially uninitialized local variables, and unreachable code
+#  pragma warning (disable: 4127 5055 4701 4702)
+#endif
+
+using namespace std;
+using namespace GeographicLib;
+
+int main() {
+  typedef GeographicLib::Math::real real;
+  try {
+    {
+      cout << fixed;
+      // GravityModel egm("egm2008");
+      // GravityModel egm("egm84-mod","/home/ckarney/geographiclib/gravity");
+      GravityModel egm("egm2008");
+      real lat, lon;
+      while (cin >> lat >> lon) {
+        real g = egm.GeoidHeight(lat, lon);
+        cout << setprecision(6)
+             << setw(12) << lat
+             << setw(12) << lon
+             << setprecision(12) << setw(20)
+             << g << "\n";
+        /*
+          real h;
+          cin >> h;
+        real Dg01, xi, eta;
+        egm.Anomaly(lat, lon, h, Dg01, xi, eta);
+        Dg01 *= 1e5;
+        xi *= Math::ds;
+        eta *= Math::ds;
+        cout << setprecision(12) << g << " " << Dg01 << " "
+             << xi << " " << eta << "\n";
+        */
+      }
+    return 0;
+    }
+    {
+      cout << fixed << setprecision(6);
+      GravityModel egm("egm2008");
+      real lat, lon;
+      while (cin >> lat >> lon) {
+        real
+          h = 0,
+          gamma = egm.ReferenceEllipsoid().SurfaceGravity(lat),
+          U0 = egm.ReferenceEllipsoid().SurfacePotential();
+        real deltax, deltay, deltaz;
+        real T = egm.Disturbance(lat, lon, h, deltax, deltay, deltaz);
+        h = T/gamma;
+        real h0 = h;
+        //        cout << "I " << U0 << " " << h << "\n";
+        for (int i = 0; i < 10; ++i) {
+          real gx, gy, gz;
+          real W = egm.Gravity(lat, lon, h, gx, gy, gz);
+          h += (W-U0)/gamma;
+          if (i == 9)
+          cout << i << " " << W << " " << h <<  " " << h - h0 << "\n";
+        }
+      }
+      return 0;
+    }
+    if (false) {
+    {
+      cout << fixed;
+      Geoid egm("egm2008-1");
+      real
+        lat0 = real(89.234412), lon0 = real(-179.234257),
+        dl = 180/real(1123), h = 0;
+      for (int j = 0; j < 2000; ++j) {
+        real lat = lat0 - j*dl/2;
+        for (int i = 0; i < 4000; ++i) {
+          real lon = lon0 + i * dl/2;
+        // cout << setprecision(3) << lon << " "
+        // << setprecision(12) << egm.GeoidHeight(lat,lon) << "\n";
+          h+=egm(lat,lon);
+        }
+      }
+      cout << setprecision(3) << 0/*lon*/ << " "
+           << setprecision(12) << h/*c.GeoidHeight(lon)*/ << "\n";
+      return 0;
+    }
+
+    {
+      cout << fixed;
+      GravityModel egm("egm2008","/home/ckarney/geographiclib/gravity");
+      real lat = 33, lon0 = 44, dlon = real(0.001);
+      GravityCircle c = egm.Circle(lat, 0.0);
+      for (int i = 0; i < 100000; ++i) {
+        real lon = lon0 + i * dlon;
+        // cout << setprecision(3) << lon << " "
+        // << setprecision(12) << egm.GeoidHeight(lat,lon) << "\n";
+        cout << setprecision(3) << lon << " "
+             << setprecision(12) << c.GeoidHeight(lon) << "\n";
+      }
+      return 0;
+    }
+    }
+    if (true) {
+      cout << setprecision(8);
+      // MagneticModel mag("/scratch/WMM2010NewLinux/WMM2010ISO.COF");
+      MagneticModel mag1("wmm2010");
+      MagneticModel mag2("emm2010");
+      MagneticModel mag3("igrf11");
+      real lat, lon, h, t, bx, by, bz, bxt, byt, bzt;
+      while (cin >> t >> lat >> lon >> h) {
+        mag1(t, lat, lon, h, bx, by, bz, bxt, byt, bzt);
+        cout << by << " " << bx << " " << -bz << " "
+             << byt << " " << bxt << " " << -bzt << "\n";
+        MagneticCircle circ1(mag1.Circle(t, lat, h));
+        circ1(lon, bx, by, bz, bxt, byt, bzt);
+        cout << by << " " << bx << " " << -bz << " "
+             << byt << " " << bxt << " " << -bzt << "\n";
+        /*
+        mag2(t, lat, lon, h, bx, by, bz, bxt, byt, bzt);
+        cout << by << " " << bx << " " << -bz << " "
+             << byt << " " << bxt << " " << -bzt << "\n";
+        MagneticCircle circ2(mag2.Circle(t, lat, h));
+        circ2(lon, bx, by, bz, bxt, byt, bzt);
+        cout << by << " " << bx << " " << -bz << " "
+             << byt << " " << bxt << " " << -bzt << "\n";
+        mag3(t, lat, lon, h, bx, by, bz, bxt, byt, bzt);
+        cout << by << " " << bx << " " << -bz << " "
+             << byt << " " << bxt << " " << -bzt << "\n";
+        MagneticCircle circ3(mag3.Circle(t, lat, h));
+        circ3(lon, bx, by, bz, bxt, byt, bzt);
+        cout << by << " " << bx << " " << -bz << " "
+             << byt << " " << bxt << " " << -bzt << "\n";
+        */
+      }
+      return 0;
+    }
+    int type = 2;
+    int N;
+    string name;
+    switch (type) {
+    case 0:
+      N = 360;
+      name = "data_EGM96_360.dat";
+      break;
+    case 1:
+      N = 2190;
+      name = "data_EGM2008_2190.dat";
+      break;
+    case 2:
+    default:
+      N = 5;
+      name = "harmtest.dat";
+      break;
+    }
+    int k = ((N + 1) * (N + 2)) / 2;
+    vector<real> C(k);
+    vector<real> S(k);
+    name = "/scratch/egm2008/harm/" + name;
+    {
+      ifstream f(name.c_str(), ios::binary);
+      if (!f.good())
+        throw GeographicErr("Cannot open coefficient file");
+      f.read(reinterpret_cast<char*>(&C[0]), k * sizeof(real));
+      f.read(reinterpret_cast<char*>(&S[0]), k * sizeof(real));
+    }
+    //    for (int i = 0; i < k; ++i)
+    //      cout << i << " " << C[i] << " " << S[i] << "\n";
+    real lat, lon;
+    cout << setprecision(17);
+    real a = real(0.9L), r = real(1.2L);
+    SphericalHarmonic harm(C, S, N, a, SphericalHarmonic::FULL);
+    vector<real> Z;
+    while (cin >> lat >> lon) {
+      real
+        phi = Math::degree() * lat,
+        lam = Math::degree() * lon,
+        x = r * (abs(lat) == 90 ? 0 : cos(phi)) * cos(lam),
+        y = r * (abs(lat) == 90 ? 0 : cos(phi)) * sin(lam),
+        z = r * sin(phi);
+      real
+        d = real(1e-7L),
+        dx1 = (harm(x+d, y, z) - harm(x-d, y, z))/(2*d),
+        dy1 = (harm(x, y+d, z) - harm(x, y-d, z))/(2*d),
+        dz1 = (harm(x, y, z+d) - harm(x, y, z-d))/(2*d),
+        dx2, dy2, dz2;
+      real
+        v1 = harm(x, y, z);
+      real
+        v2 = harm(x, y, z, dx2, dy2, dz2);
+      cout << v1 << " " << v2 << "\n";
+      cout << dx1 << " " << dx2 << "\n"
+           << dy1 << " " << dy2 << "\n"
+           << dz1 << " " << dz2 << "\n";
+      CircularEngine circ1 = harm.Circle(r * cos(phi), z, false);
+      CircularEngine circ2 = harm.Circle(r * cos(phi), z, true);
+      real v3, dx3, dy3, dz3;
+      v3 = circ2(lon, dx3, dy3, dz3);
+      cout << v3 << " " << dx3 << " " << dy3 << " " << dz3 << "\n";
+    }
+    cout << "start timing" << endl;
+    real phi, lam, sum, z, p, dx, dy, dz;
+    lat = 33; phi = lat * Math::degree();
+    z = r * sin(phi);
+    p = r * cos(phi);
+    sum = 0;
+    for (int i = 0; i < 100; ++i) {
+      lam = (44 + real(0.01)*i) * Math::degree();
+      sum += harm(p * cos(lam), p * sin(lam), z, dx, dy, dz);
+    }
+    cout << "sum a " << sum << endl;
+    CircularEngine circ(harm.Circle(p, z, true));
+    sum = 0;
+    for (int i = 0; i < 100000; ++i) {
+      lon = (44 + real(0.01)*i);
+      sum += circ(lon, dx, dy, dz);
+    }
+    cout << "sum b " << sum << endl;
+  }
+  catch (const exception& e) {
+    cerr << "Caught exception: " << e.what() << "\n";
+    return 1;
+  }
+  catch (...) {
+    cerr << "Caught unknown exception\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/LevelEllipsoid.cpp

@@ -0,0 +1,403 @@
+#include <iostream>
+#include <iomanip>
+#include <vector>
+#include <array>
+#include <limits>
+#include <algorithm>
+#include <boost/numeric/odeint.hpp>
+
+#include <GeographicLib/Math.hpp>
+#include <GeographicLib/Utility.hpp>
+#include <GeographicLib/NormalGravity.hpp>
+
+using namespace GeographicLib;
+namespace ode = boost::numeric::odeint;
+typedef Math::real real;
+typedef std::array<real, 2> point;
+
+// Initialize with the end points of a line
+class LineDistance {
+  point _a, _b;
+  real _l, _nx, _ny;
+public:
+  LineDistance(const point& a, const point& b);
+  real Distance(const point& p) const;
+  // negative on left of line, pos on right
+  real Displacement(const point& p) const;
+};
+
+// Use Ramer-Douglas-Peucker to simplify a polyline
+class LineSimplifier {
+  // Simplify range [a,b]
+  static void InternalSimplify(std::vector<point>& p, unsigned a, unsigned b,
+                               real thresh);
+public:
+  static void Simplify(std::vector<point>& p, real thresh);
+};
+
+inline LineDistance::LineDistance(const point& a, const point& b)
+  : _a(a), _b(b) {
+  using std::hypot;
+  real dx = _b[0] - _a[0], dy = _b[1] - _a[1];
+  _l = hypot(dx, dy);
+  if (_l > std::numeric_limits<real>::epsilon()) {
+  _nx = dx / _l; _ny = dy / _l;
+  } else {
+    _l = 0; _nx = 1; _ny = 0;
+  }
+}
+
+inline real LineDistance::Distance(const point& p) const {
+  using std::abs; using std::hypot;
+  real x = p[0] - _a[0], y = p[1] - _a[1];
+  if (_l != 0) {
+    real X = x * _nx + y * _ny, Y = abs(x * _ny - y * _nx);
+    X = X < 0 ? -X : (X > _l ? X - _l : 0);
+    return X != 0 ? hypot(X, Y) : Y;
+  } else
+    return hypot(x, y);
+}
+
+inline real LineDistance::Displacement(const point& p) const {
+  real x = p[0] - _a[0], y = p[1] - _a[1];
+  return x * _ny - y * _nx;
+}
+
+inline void LineSimplifier::Simplify(std::vector<point>& p, real thresh) {
+  using std::isnan;
+  unsigned n = p.size();
+  InternalSimplify(p, 0, n-1, thresh);
+  unsigned i = 0, j = 0;
+  while (i < n) {               // Squeeze out nans
+    if (!isnan(p[i][0])) { p[j] = p[i]; ++j; }
+    ++i;
+  }
+  p.resize(j);
+}
+
+void LineSimplifier::InternalSimplify(std::vector<point>& p,
+                                      unsigned a, unsigned b,
+                                      real thresh) {
+  if (b - a + 1 <= 2) return;
+  // Assume b > a+1, thresh > 0
+  real maxdist = -1;
+  unsigned maxind = (a + b) / 2;         // initial value isn't used
+  LineDistance dist(p[a], p[b]);
+  for (unsigned i = a + 1; i <= b - 1; ++i) {
+    real d = dist.Distance(p[i]);
+    if (d > maxdist) { maxdist = d; maxind = i; }
+  }
+  if (maxdist > thresh) {
+    InternalSimplify(p, a, maxind, thresh);
+    InternalSimplify(p, maxind, b, thresh);
+  } else {
+    for (unsigned i = a+1; i <= b-1; ++i)
+      p[i][0] = Math::NaN();
+  }
+}
+
+class PointTest {
+  real _xmin, _xmax, _ymin, _ymax;
+public:
+  PointTest(real xmin, real xmax, real ymin, real ymax)
+    : _xmin(xmin)
+    , _xmax(xmax)
+    , _ymin(ymin)
+    , _ymax(ymax) {}
+  bool operator() (const point &p) const {
+    return p[0] >= _xmin && p[0] <= _xmax && p[1] >= _ymin && p[1] <= _ymax;
+  }
+};
+
+class GravityInt {
+  const NormalGravity& _grav;
+  unsigned _dir;
+public:
+  GravityInt(const NormalGravity& grav, unsigned dir)
+    : _grav(grav)
+    , _dir(dir & 3u)
+  {}
+  void operator() (const point& p, point &d, const real /*t*/) const {
+    real Y = 0, gX, gY, gZ;
+    _grav.U(p[0], Y, p[1], gX, gY, gZ);
+    Math::norm(gX, gZ);
+    switch (_dir) {
+    case  0: d[0] =  gX; d[1] =  gZ; break; // parallel to g
+    case  1: d[0] = -gZ; d[1] =  gX; break; // counterclockwise 90d from g
+    case  2: d[0] = -gX; d[1] = -gZ; break; // antiparallel to g
+    case  3:
+    default: d[0] =  gZ; d[1] = -gX; break; // clockwise 90d from g
+    }
+  }
+};
+
+class GravityFollow {
+public:
+  static void follow(const GravityInt& sys, const point& p0, real t0, real ds,
+                     std::vector<point>& points, const PointTest& box) {
+    ode::result_of::make_dense_output
+      < ode::runge_kutta_dopri5< point, real > >::type integrator =
+      ode::make_dense_output(real(1.0e-8), real(0*1.0e-8),
+                             ode::runge_kutta_dopri5< point, real >() );
+    integrator.initialize(p0, t0, real(1e-2));
+    integrator.do_step(sys);
+    int n = 1, i = 1;
+    point out;
+    while (true) {
+      real tout = i * ds;
+      while (integrator.current_time() < tout) {
+        integrator.do_step(sys); ++n;
+      }
+      integrator.calc_state(tout, out);
+      points.push_back(out);
+      if (!box(out)) break;
+      ++i;
+    }
+  }
+};
+
+// Evaluate gravity potential or gradient along 1 axis
+class GravityEval {
+  const NormalGravity& _grav;
+  bool _zaxis, _grad;
+  real _disp;
+public:
+  GravityEval(const NormalGravity& grav, bool zaxis, bool grad, real disp)
+    : _grav(grav)
+    , _zaxis(zaxis)
+    , _grad(grad)
+    , _disp(disp) {};
+  real operator() (real v) const {
+    real
+      X = _zaxis ? _disp : v,
+      Y = 0,
+      Z = _zaxis ? v : _disp,
+      gX, gY, gZ,
+      U = _grav.U(X, Y, Z, gX, gY, gZ);
+    return _grad ? gX : U;
+  }
+};
+
+class Interpolator {
+private:
+  int _maxit;
+  real _eps;
+public:
+  Interpolator()
+    : _maxit(20), _eps(sqrt(std::numeric_limits<real>::epsilon())) {}
+  real operator() (const GravityEval& gravfun,
+                   real x0, real x1, real val) {
+    using std::abs;
+    real y0 = gravfun(x0) - val, y1;
+    for (int i = 0, trip = 0; i < _maxit; ++i) {
+      y1 = gravfun(x1) - val;
+      if (y1 == y0) break;
+      real x2 = x1 - y1 * (x1 - x0) / (y1 - y0);
+      x0 = x1; x1 = x2; y0 = y1;
+      if (abs(x0 - x1) <= _eps) ++trip;
+      if (trip > 2) break;
+    }
+    return x1;
+  }
+};
+
+void dump(const std::vector<point>& points) {
+  for (auto p = points.begin(); p != points.end(); ++p) {
+    std::cout << (*p)[0] << "," << (*p)[1] << ";\n";
+  }
+}
+
+int main() {
+  Utility::set_digits();
+  using std::sqrt;
+  real a = 1, GM = 1, omega = 3/Math::real(10), f = 2/Math::real(10);
+  real gX, gY, gZ, b = (1 - f) * a, E = a * sqrt(f * (2 - f));
+  real thresh = real(0.5e-3);
+  NormalGravity grav(a, GM, omega, f, true);
+  real xmax = 3, ymax = 2;
+  PointTest box(0, xmax, 0, ymax);
+  Interpolator intpol;
+  real
+    X0 = a,
+    U0 = grav.U(X0, 0, 0, gX, gY, gZ),
+    Uc = grav.U(0, 0, 0, gX, gY, gZ);
+  real Xnull, Unull;
+  {
+    GravityEval eval(grav, false, true, 0);
+    Xnull = intpol(eval, 1, 2, 0);
+    Unull = grav.U(Xnull, 0, 0, gX, gY, gZ);
+  }
+  int ndiv = 20;
+  real del = (U0 - Unull) / 20;
+  std::cout << std::setprecision(6);
+  std::cout << "a=" << a << "; b=" << b << "; Rnull=" << Xnull
+            << "; xmax=" << xmax << "; ymax=" << ymax << ";\n";
+  std::cout << "q={}; p={}; qa={}; pa={};\n";
+  std::vector<point> points;
+  point p0;
+  int k = 0;
+  for (int i = 0; i <= 90; i += 10) {
+    points.clear();
+    real s, c;
+    Math::sincosd(real(i), s, c);
+    p0[0] = a * c; p0[1] = b * s;
+    points.push_back(p0);
+    if (i == 0) {
+      p0[0] = xmax; p0[1] = 0;
+      points.push_back(p0);
+    } else if (i == 90) {
+      p0[0] = 0; p0[1] = ymax;
+      points.push_back(p0);
+    } else {
+      GravityInt sys(grav, 2u);
+      GravityFollow::follow(sys, p0, real(0), 1/real(100), points, box);
+    }
+    LineSimplifier::Simplify(points, thresh);
+    std::cout << "q{" << ++k << "}=[\n";
+    dump(points);
+    std::cout << "];\n";
+  }
+  {
+    points.clear();
+    p0[0] = Xnull; p0[1] = 0;
+    points.push_back(p0);
+    p0[0] = Xnull; p0[1] = real(0.001);
+    points.push_back(p0);
+    GravityInt sys(grav, 2u);
+    GravityFollow::follow(sys, p0, real(0), 1/real(100), points, box);
+    LineSimplifier::Simplify(points, thresh);
+    std::cout << "q{" << ++k << "}=[\n";
+    dump(points);
+    std::cout << "];\n";
+  }
+  for (int i = 1; i <= 5; ++i) {
+    points.clear();
+    p0[0] = xmax;
+    p0[1] = real(i == 1 ? 0.45 :
+                 (i == 2 ? 0.9 :
+                  (i == 3 ? 1.25 :
+                   (i == 4 ? 1.6 : 1.9))));
+    points.push_back(p0);
+    if (!box(p0)) break;
+    GravityInt sys(grav, 2u);
+    GravityFollow::follow(sys, p0, real(0), 1/real(100), points, box);
+    LineSimplifier::Simplify(points, thresh);
+    std::cout << "q{" << ++k << "}=[\n";
+    dump(points);
+    std::cout << "];\n";
+  }
+  k = 0;
+  for (int i = 0; i <= 90; i += 10) {
+    points.clear();
+    real s, c;
+    Math::sincosd(real(i), s, c);
+    p0[0] = a * c; p0[1] = b * s;
+    points.push_back(p0);
+    if (i == 0) {
+      p0[0] = E; p0[1] = 0;
+      points.push_back(p0);
+    } else if (i == 90) {
+      p0[0] = 0; p0[1] = 0;
+      points.push_back(p0);
+    } else {
+      GravityInt sys(grav, 0u);
+      GravityFollow::follow(sys, p0, real(0), 1/real(100), points, box);
+    }
+    LineSimplifier::Simplify(points, thresh);
+    std::cout << "qa{" << ++k << "}=[\n";
+    dump(points);
+    std::cout << "];\n";
+  }
+  k = 0;
+  {
+    GravityEval eval(grav, false, false, 0);
+    for (int i = 0; i < ndiv; ++i) {
+      points.clear();
+      real U = U0 - del * i;
+      p0[0] = intpol(eval, X0, Xnull, U); p0[1] = 0;
+      points.push_back(p0);
+      GravityInt sysa(grav, 3u);
+      GravityFollow::follow(sysa, p0, real(0), 1/real(100), points, box);
+      LineSimplifier::Simplify(points, thresh);
+      std::cout << "p{" << ++k << "}=[\n";
+      dump(points);
+      std::cout << "];\n";
+    }
+    for (int i = ndiv - 1;; --i) {
+      points.clear();
+      real U = U0 - del * i;
+      p0[0] = intpol(eval, Xnull + 3, Xnull, U); p0[1] = 0;
+      if (!box(p0)) break;
+      points.push_back(p0);
+      GravityInt sysa(grav, 1u);
+      GravityFollow::follow(sysa, p0, real(0), 1/real(100), points, box);
+      LineSimplifier::Simplify(points, thresh);
+      std::cout << "p{" << ++k << "}=[\n";
+      dump(points);
+      std::cout << "];\n";
+    }
+  }
+  {
+    GravityEval eval(grav, true, false, 0);
+    for (int i = ndiv + 1;; ++i) {
+      points.clear();
+      real U = U0 - del * i;
+      p0[1] = intpol(eval, X0, Xnull, U); p0[0] = 0;
+      if (!box(p0)) break;
+      points.push_back(p0);
+      GravityInt sysa(grav, 1u);
+      GravityFollow::follow(sysa, p0, real(0), 1/real(100), points, box);
+      LineSimplifier::Simplify(points, thresh);
+      std::cout << "p{" << ++k << "}=[\n";
+      dump(points);
+      std::cout << "];\n";
+    }
+  }
+  {
+    real dy = real(0.02);
+    GravityEval eval(grav, false, false, dy);
+    points.clear();
+    p0[0] = Xnull; p0[1] = 0;
+    points.push_back(p0);
+    p0[0] = intpol(eval, Xnull - real(0.1), Xnull, Unull); p0[1] = dy;
+    points.push_back(p0);
+    GravityInt sysa(grav, 3u);
+    GravityFollow::follow(sysa, p0, real(0), 1/real(100), points, box);
+    LineSimplifier::Simplify(points, thresh);
+    std::cout << "p{" << ++k << "}=[\n";
+    dump(points);
+    std::cout << "];\n";
+  }
+  {
+    real dy = real(0.02);
+    GravityEval eval(grav, false, false, dy);
+    points.clear();
+    p0[0] = Xnull; p0[1] = 0;
+    points.push_back(p0);
+    p0[0] = intpol(eval, Xnull + real(0.1), Xnull, Unull); p0[1] = dy;
+    points.push_back(p0);
+    GravityInt sysa(grav, 1u);
+    GravityFollow::follow(sysa, p0, real(0), 1/real(100), points, box);
+    LineSimplifier::Simplify(points, thresh);
+    std::cout << "p{" << ++k << "}=[\n";
+    dump(points);
+    std::cout << "];\n";
+  }
+  k = 0;
+  {
+    GravityEval eval(grav, false, false, 0);
+    for (int i = 1;; ++i) {
+      points.clear();
+      real U = U0 + 5 * del * i;
+      if (U > Uc) break;
+      p0[0] = intpol(eval, 0, X0, U); p0[1] = 0;
+      points.push_back(p0);
+      GravityInt sysa(grav, 3u);
+      GravityFollow::follow(sysa, p0, real(0), 1/real(100), points, box);
+      LineSimplifier::Simplify(points, thresh);
+      std::cout << "pa{" << ++k << "}=[\n";
+      dump(points);
+      std::cout << "];\n";
+    }
+  }
+}

libs/geographiclib/develop/LevelEllipsoid.m

@@ -0,0 +1,57 @@
+% Assume output from LevelEllipsoid has been read in
+doinside=1;
+red=[179,27,27]/255;
+white=[1,1,1];
+black=[0,0,0];
+blue=[0,19,56]/100;
+green=[9,45,27]/100;
+gray=[0.9,0.9,0.9];
+thick=1;
+pltsize=[5.9 4.12];
+set(gcf,'Units','pixels');
+set(gcf,'Position',50+150*[0 0 pltsize]);
+set(gcf,'Units','normalized');
+hold off;
+if ~doinside
+  fill([p{1}(:,1);0],[p{1}(:,2);0], gray, 'EdgeColor', 'none');
+  hold on
+end
+nq=size(q,2);
+for i=1:nq;
+  plot(q{i}(:,1),q{i}(:,2), 'Color', green);
+  if i == 1, hold on; end
+end
+np=size(p,2);
+for i=1:np
+  color = blue;
+  if i == 1
+    plot(p{i}(:,1),p{i}(:,2), 'Color', red, 'LineWidth', thick);
+  else
+    plot(p{i}(:,1),p{i}(:,2), 'Color', blue);
+  end
+end
+if doinside
+  nq=size(qa,2);
+  for i=1:nq;
+    plot(qa{i}(:,1),qa{i}(:,2), 'Color', green);
+  end
+  np=size(pa,2);
+  for i=1:np
+    plot(pa{i}(:,1),pa{i}(:,2), 'Color', blue);
+  end
+  plot([0, 0.6], [0, 0], 'Color', black, 'LineWidth', thick);
+end
+hold off
+xlabel('R'); ylabel('Z');
+axis equal;
+axis([0,xmax,0,ymax]);
+set(gcf,'PaperOrientation', 'landscape');
+set(gcf,'PaperSize',pltsize);
+set(gcf,'PaperPosition',[0 0 pltsize]);
+set(gca, 'XTick',[0:3]);
+set(gca, 'YTick',[0:2]);
+set(gca, 'LooseInset',[0.07 0.09 0.03 0.02]);
+ylabelh=get(gca,'ylabel');
+set(ylabelh,'rotation',0);
+set(ylabelh,'Position', get(ylabelh, 'Position') + [-0.1 0.2 0]);
+print('-dsvg', ['normal-gravity-potential-', num2str(doinside), '.svg']);

libs/geographiclib/develop/M12zero.cpp

@@ -0,0 +1,42 @@
+#include <iostream>
+#include <GeographicLib/Utility.hpp>
+#include <GeographicLib/GeodesicExact.hpp>
+#include <GeographicLib/GeodesicLineExact.hpp>
+
+int main(int argc, const char* const argv[]) {
+  try {
+    using namespace GeographicLib;
+    typedef Math::real real;
+    Utility::set_digits();
+    real a = 2 / Math::pi();
+    if (argc != 2) {
+      std::cout << "Usage: M12zero f\n";
+      return 1;
+    }
+    real f = Utility::fract<real>(std::string(argv[1]));
+
+    GeodesicExact g(a, f);
+    GeodesicLineExact l(g, 90, 0, 0);
+    real s12 = a * Math::pi()/2;
+    real lat2, lon2, azi2, m12, M12, M21, ds12;
+    for (int i = 0; i < 20; ++i) {
+      using std::abs;
+      l.Position(s12, lat2, lon2, azi2, m12, M12, M21);
+      ds12 = m12 * M12 / (1 - M12 * M21);
+      s12 = abs(s12 + ds12);
+    }
+    real q;
+    g.Inverse(0,0,90,0,q);
+    int p = 16 + Math::extra_digits();
+    p = 20;
+    std::cout << std::fixed << std::setprecision(p) << s12 << " " << q << "\n";
+  }
+  catch (const std::exception& e) {
+    std::cerr << "Caught exception: " << e.what() << "\n";
+    return 1;
+  }
+  catch (...) {
+    std::cerr << "Caught unknown exception\n";
+    return 1;
+  }
+}

libs/geographiclib/develop/NaNTester.cpp

@@ -0,0 +1,126 @@
+/**
+ * \file ProjTest.cpp
+ * \brief Command line utility for testing transverse Mercator projections
+ *
+ * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ **********************************************************************/
+
+#include "GeographicLib/EllipticFunction.hpp"
+#include "GeographicLib/TransverseMercator.hpp"
+#include "GeographicLib/TransverseMercatorExact.hpp"
+#include "GeographicLib/PolarStereographic.hpp"
+#include <iostream>
+#include <cassert>
+
+using namespace GeographicLib;
+using namespace std;
+
+int main() {
+  {
+    Math::real x, y, gamma, k;
+    x = y = gamma = k = 0;
+    PolarStereographic::UPS().Forward(true, Math::NaN<Math::real>(), 30.0,
+                                    x, y, gamma, k);
+    cout << x << " " << y << " " << k << "\n";
+    x = y = gamma = k = 0;
+    PolarStereographic::UPS().Forward(false, -80.0, Math::NaN<Math::real>(),
+                                    x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << "\n";
+  }
+  {
+    Math::real lat, lon, gamma, k;
+    lat = lon = gamma = k = 0;
+    PolarStereographic::UPS().Reverse(true, Math::NaN<Math::real>(), 0.0,
+                                    lat, lon, gamma, k);
+    cout << lat << " " << lon << " " << gamma << " " << k << "\n";
+    lat = lon = gamma = k = 0;
+    PolarStereographic::UPS().Reverse(false,0.0, Math::NaN<Math::real>(),
+                                    lat, lon, gamma, k);
+    cout << lat << " " << lon << " " << gamma << " " << k << "\n";
+  }
+  {
+    Math::real x, y, gamma, k;
+    x = y = gamma = k = 0;
+    TransverseMercator::UTM().Forward(0.0, Math::NaN<Math::real>(), 0.0,
+                                    x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << " " << k << "\n";
+    x = y = gamma = k = 0;
+    TransverseMercator::UTM().Forward(0.0, 0.0, Math::NaN<Math::real>(),
+                                    x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << " " << k << "\n";
+    x = y = gamma = k = 0;
+    TransverseMercator::UTM().Forward(Math::NaN<Math::real>(), 0.0, 0.0,
+                                    x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << " " << k << "\n";
+  }
+  {
+    Math::real lat, lon, gamma, k;
+    lat = lon = gamma = k = 0;
+    TransverseMercator::UTM().Reverse(0.0, Math::NaN<Math::real>(), 0.0,
+                                    lat, lon, gamma, k);
+    cout << lat << " " << lon << " " << gamma << " " << k << "\n";
+    lat = lon = gamma = k = 0;
+    TransverseMercator::UTM().Reverse(0.0, 0.0, Math::NaN<Math::real>(),
+                                    lat, lon, gamma, k);
+    cout << lat << " " << lon << " " << gamma << " " << k << "\n";
+    lat = lon = gamma = k = 0;
+    TransverseMercator::UTM().Reverse(Math::NaN<Math::real>(), 0.0, 0.0,
+                                    lat, lon, gamma, k);
+    cout << lon << "\n";
+  }
+  {
+    Math::real x, y, gamma, k;
+    x = y = gamma = k = 0;
+    TransverseMercatorExact::UTM().Forward(0.0, Math::NaN<Math::real>(), 0.0,
+                                         x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << " " << k << "\n";
+    x = y = gamma = k = 0;
+    TransverseMercatorExact::UTM().Forward(0.0, 0.0, Math::NaN<Math::real>(),
+                                         x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << " " << k << "\n";
+    x = y = gamma = k = 0;
+    TransverseMercatorExact::UTM().Forward(Math::NaN<Math::real>(), 0.0, 0.0,
+                                         x, y, gamma, k);
+    cout << x << " " << y << " " << gamma << " " << k << "\n";
+  }
+  {
+    Math::real lat, lon, gamma, k;
+    lat = lon = gamma = k = 0;
+    TransverseMercatorExact::UTM().Reverse(0.0, Math::NaN<Math::real>(), 0.0,
+                                         lat, lon, gamma, k);
+    cout << lat << " " << lon << " " << gamma << " " << k << "\n";
+    lat = lon = gamma = k = 0;
+    TransverseMercatorExact::UTM().Reverse(0.0, 0.0, Math::NaN<Math::real>(),
+                                         lat, lon, gamma, k);
+    cout << lat << " " << lon << " " << gamma << " " << k << "\n";
+    lat = lon = gamma = k = 0;
+    TransverseMercatorExact::UTM().Reverse(Math::NaN<Math::real>(), 0.0, 0.0,
+                                         lat, lon, gamma, k);
+    cout << lon << "\n";
+  }
+  {
+    EllipticFunction e(Math::NaN<Math::real>());
+    Math::real sn, cn, dn;
+    cout << " k2 " << e.k2()
+         << " kp2 " << e.kp2()
+         << " K " << e.K()
+         << " E " << e.E()
+         << " KE " << e.KE() << "\n";
+    sn = cn = dn = 0;
+    e.sncndn(Math::real(0.1), sn, cn, dn);
+    cout << " sncndn " << sn << " " << cn << " " << dn
+         << " E(phi) " << e.E(Math::real(0.1))
+         << " E(sncndn) " << e.E(sn, cn, dn) << "\n";
+  }
+  {
+    EllipticFunction e(Math::real(0.1));
+    Math::real sn, cn, dn;
+    sn = cn = dn = 0;
+    e.sncndn(Math::NaN<Math::real>(), sn, cn, dn);
+    cout << " sncndn " << sn << " " << cn << " " << dn
+         << " E(phi) " << e.E(Math::NaN<Math::real>())
+         << " E(sncndn) " << e.E(sn, cn, dn) << "\n";
+  }
+}

libs/geographiclib/develop/NormalTest.cpp

@@ -0,0 +1,134 @@
+#include <iostream>
+#include <iomanip>
+#include <limits>
+#include <vector>
+#include <random>
+#include <algorithm>
+
+#include <GeographicLib/NormalGravity.hpp>
+#include <GeographicLib/Ellipsoid.hpp>
+#include <GeographicLib/Utility.hpp>
+
+using namespace std;
+using namespace GeographicLib;
+
+void checkdiff(const NormalGravity& g,
+               Math::real X, Math::real Y, Math::real Z,
+               Math::real eps, Math::real tol) {
+  Math::real gX, gY, gZ, t,
+    V   = g.V0(X, Y, Z, gX, gY, gZ),
+    VXp = g.V0(X+eps, Y, Z, t, t, t), VXm = g.V0(X-eps, Y, Z, t, t, t),
+    VYp = g.V0(X, Y+eps, Z, t, t, t), VYm = g.V0(X, Y-eps, Z, t, t, t),
+    VZp = g.V0(X, Y, Z+eps, t, t, t), VZm = g.V0(X, Y, Z-eps, t, t, t),
+    ggX = (VXp - VXm) / (2*eps),
+    ggY = (VYp - VYm) / (2*eps),
+    ggZ = (VZp - VZm) / (2*eps),
+    lap = (VXp + VXm + VYp + VYm + VZp + VZm - 6 * V) / Math::sq(eps);
+  if (!(abs(ggX - gX) < tol &&
+        abs(ggY - gY) < tol &&
+        abs(ggZ - gZ) < tol &&
+        abs(lap) < tol))
+    cout << "Failure with f = " << g.Flattening() << " X,Y,Z = "
+         << X << " " << Y << " " << Z << "\n"
+         << "     " << ggX - gX << " " << ggY - gY << " " << ggZ - gZ
+         << " " << lap << "\n";
+}
+
+// Copied from NormalGravity (w/o the series expansion)
+Math::real atanzz(Math::real x, bool alt) {
+  // This routine obeys the identity
+  //   atanzz(x, alt) = atanzz(-x/(1+x), !alt)
+  //
+  // Require x >= -1.  Best to call with alt, s.t. x >= 0; this results in
+  // a call to atan, instead of asin, or to asinh, instead of atanh.
+  Math::real z = sqrt(abs(x));
+  return x == 0 ? 1 :
+    (alt
+     ? (!(x < 0) ? asinh(z) : asin(z)) / sqrt(abs(x) / (1 + x))
+     : (!(x < 0) ? atan(z) : atanh(z)) / z);
+}
+
+// Copied from NormalGravity (w/o the series expansion)
+Math::real Qf(Math::real x, bool alt) {
+  // Compute
+  //   Q(z) = (((1 + 3/z^2) * atan(z) - 3/z)/2) / z^3
+  //        = q(z)/z^3 with q(z) defined by H+M, Eq 2-57 with z = E/u
+  //   z = sqrt(x)
+  Math::real y = alt ? -x / (1 + x) : x;
+  return ((1 + 3/y) * atanzz(x, alt) - 3/y) / (2 * y);
+}
+
+int main() {
+  Utility::set_digits();
+  Math::real
+    eps = sqrt(sqrt(numeric_limits<Math::real>::epsilon())),
+    tol = eps;
+  cout << setprecision(13);
+  cout << "eps = " << eps << "\n";
+  unsigned s = random_device()(); // Set seed from random_device
+  mt19937 r(s);                   // Initialize URNG
+  if (1)  {
+    Math::real GM = 1, a = 1, omega = 1;
+    {
+      Math::real f = 1/Math::real(5);
+      NormalGravity g(a, GM, omega, f, true);
+      Math::real gX, gY, gZ, U, X = Math::real(5)/10, Y = 0,
+        R = hypot(X, Y),
+        b = a*(1-f),  E2 = abs(a*a - b*b), E = sqrt(E2);
+      U = g.V0(X, Y, 0, gX, gY, gZ);
+      cout << U << " " << gX << " " << gY << " " << gZ << "\n";
+      U = g.V0(X, Y, eps, gX, gY, gZ);
+      cout << U << " " << gX << " " << gY << " " << gZ << "\n";
+      cout << "gz = "
+           << -(GM/(E*sqrt(E2-R*R)) +
+                Math::sq(a*b*omega)*b /
+                (E2*E2*E*Qf(Math::sq(E/b), false)) *
+                (2*E2/3 - R*R)/sqrt(E2-R*R)) << "\n";
+    }
+    {
+      Math::real f = -1/Math::real(5);
+      NormalGravity g(a, GM, omega, f, true);
+      Math::real gX, gY, gZ, U, Z = Math::real(4)/10, Y = 0,
+        b = a*(1-f),  E2 = abs(a*a - b*b), E = sqrt(E2);
+      U = g.V0(0, Y, Z, gX, gY, gZ);
+      cout << U << " " << gX << " " << gY << " " << gZ << "\n";
+      U = g.V0(eps, Y, Z, gX, gY, gZ);
+      cout << U << " " << gX << " " << gY << " " << gZ << "\n";
+      cout << "-R*gR = " << eps*hypot(gX, gY) << " "
+           << 2*(GM/(2*E) + Math::sq(a*b*omega)*b /
+                 (E2*E2*E*Qf(Math::sq(E/a), true)) *
+                 (3*Z*Z - E2)/12) << "\n";
+    }
+    {
+      Math::real f = 0;
+      NormalGravity g(a, GM, omega, f, true);
+      Math::real gX, gY, gZ, U,
+        X = Math::real(1)/10, Z = Math::real(4)/10, Y = 0,
+        b = a*(1-f), R = hypot(X, Z), beta = atan2(Z, X);
+      U = g.V0(X, Y, Z, gX, gY, gZ);
+      cout << U << " " << gX << " " << gY << " " << gZ << "\n";
+      cout << "U = " << GM/R + Math::sq(a*b*omega)*b/(R*R*R) *
+        (3*Math::sq(sin(beta)) - 1)/6 << "\n";
+      cout << "gR = " << gX*cos(beta) + gZ*sin(beta) << " "
+           << -GM/(R*R) - (Math::sq(a*b*omega)*b/(R*R*R*R) *
+                           (3*Math::sq(sin(beta)) - 1)/2)
+           << " gbeta = " << -gX*sin(beta) + gZ*cos(beta) << " "
+           << Math::sq(a*b*omega)*b/(R*R*R*R) * sin(beta)*cos(beta) << "\n";
+    }
+    return 0;
+  }
+  uniform_real_distribution<double> dis;
+  vector<Math::real> f = {-99, -0.1, 0, 0.1, 0.99};
+  Math::real GM = 1, omega = 1;
+  int num = 1000;
+  for (size_t i = 0; i < f.size(); ++i) {
+    Math::real a = 1 / sqrt(1 - f[i]);
+    cout << a << " " << a*(1-f[i]) << "\n";
+    NormalGravity g(a, GM, omega, f[i], true);
+    for (int j = 0; j < num; ++j) {
+      Math::real
+        X = Math::real(dis(r)), Y = Math::real(dis(r)), Z = Math::real(dis(r));
+      checkdiff(g, X, Y, Z, eps, tol);
+    }
+  }
+}

libs/geographiclib/develop/ProjTest.cpp

@@ -0,0 +1,343 @@
+/**
+ * \file ProjTest.cpp
+ * \brief Command line utility for testing transverse Mercator projections
+ *
+ * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ **********************************************************************/
+
+#include "GeographicLib/LambertConformalConic.hpp"
+#include "GeographicLib/PolarStereographic.hpp"
+#include "GeographicLib/TransverseMercator.hpp"
+#include "GeographicLib/TransverseMercatorExact.hpp"
+#include "GeographicLib/Constants.hpp"
+#include <string>
+#include <limits>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <iomanip>
+#include <stdexcept>
+
+class DataFile {
+private:
+  typedef GeographicLib::Math::real real;
+  std::istream* _istr;
+  std::ifstream _fstr;
+  //  String headers
+  //    COORDINATES
+  //    PROJECTION
+  //    DATUM
+  std::string _coords, _proj, _datum;
+  //  Numeric headers
+  //    CENTRAL MERIDIAN
+  //    FALSE EASTING
+  //    FALSE NORTHING
+  //    LATITUDE OF TRUE SCALE
+  //    LONGITUDE DOWN FROM POLE
+  //    ORIGIN LATITUDE
+  //    SCALE FACTOR
+  //    STANDARD PARALLEL ONE
+  //    STANDARD PARALLEL TWO
+  real _lon0, _fe, _fn, _latts, _lonfp, _lat0, _k0, _lat1, _lat2;
+  DataFile(const DataFile&);
+  DataFile& operator=(const DataFile&);
+public:
+  DataFile(const std::string& file) {
+    if (!(file.empty() || file == "-")) {
+      _fstr.open(file.c_str());
+      if (!_fstr.good())
+        throw std::out_of_range("Cannot open " + file);
+    }
+
+    _istr = (file.empty() || file == "-") ? &std::cin : &_fstr;
+    _coords = _proj = _datum = "NONE";
+    _lon0 = _fe = _fn = _latts = _lonfp = _lat0 = _k0 = _lat1 = _lat2
+      = GeographicLib::Math::NaN<real>();
+    std::string s;
+    while (getline(*_istr, s)) {
+      if (s.empty() || s[0] == '#')
+        continue;
+      if (s.substr(0, 13) == "END OF HEADER")
+        break;
+      std::string::size_type p = s.find(':');
+      if (p == std::string::npos)
+        continue;
+      std::string key(s, 0, p);
+      p = s.find_first_not_of(" \t\v\r\n\f", p + 1);
+      if (p == std::string::npos)
+        continue;
+      s = s.substr(p);
+      p = s.find_last_not_of(" \t\v\r\n\f");
+      s = s.substr(0, p + 1);
+      std::istringstream is(s);
+      if (key == "COORDINATES") {
+        _coords = s;
+      } else if (key == "PROJECTION") {
+        _proj = s;
+      } else if (key == "DATUM") {
+        _datum = s;
+      } else if (key == "CENTRAL MERIDIAN") {
+        is >> _lon0;
+      } else if (key == "FALSE EASTING") {
+        is >> _fe;
+      } else if (key == "FALSE NORTHING") {
+        is >> _fn;
+      } else if (key == "LATITUDE OF TRUE SCALE") {
+        is >> _latts;
+      } else if (key == "LONGITUDE DOWN FROM POLE") {
+        is >> _lonfp;
+      } else if (key == "ORIGIN LATITUDE") {
+        is >> _lat0;
+      } else if (key == "SCALE FACTOR") {
+        is >> _k0;
+      } else if (key == "STANDARD PARALLEL ONE") {
+        is >> _lat1;
+      } else if (key == "STANDARD PARALLEL TWO") {
+        is >> _lat2;
+      }
+    }
+  }
+  bool Next(real &x, real &y) {
+    char c;
+    // Avoid a warning about void* changed to bool
+    return (*_istr >> x >> c >> y) ? true : false;
+  }
+  bool Next(real &x, real &y, real &z) {
+    char c, d;
+    // Avoid a warning about void* changed to bool
+    return (*_istr >> x >> c >> y >> d >> z) ? true : false;
+  }
+  const std::string& coords() const { return _coords; }
+  const std::string& proj() const { return _proj; }
+  const std::string& datum() const { return _datum; }
+  real lon0() const { return _lon0; }
+  real fe() const { return _fe; }
+  real fn() const { return _fn; }
+  real latts() const { return _latts; }
+  real lonfp() const { return _lonfp; }
+  real lat0() const { return _lat0; }
+  real k0() const { return _k0; }
+  real lat1() const { return _lat1; }
+  real lat2() const { return _lat2; }
+};
+
+GeographicLib::Math::real
+dist(GeographicLib::Math::real a, GeographicLib::Math::real f,
+     GeographicLib::Math::real lat0, GeographicLib::Math::real lon0,
+     GeographicLib::Math::real lat1, GeographicLib::Math::real lon1) {
+  using namespace GeographicLib;
+  using std::cos; using std::sin; using std::sqrt; using std::hypot;
+  typedef Math::real real;
+  real
+    phi = lat0 * Math::degree(),
+    e2 = f * (2 - f),
+    sinphi = sin(phi),
+    n = 1/sqrt(1 - e2 * sinphi * sinphi),
+      // See Wikipedia article on latitude
+    hlon = cos(phi) * n,
+    hlat = (1 - e2) * n * n * n,
+    dlon = lon1 - lon0;
+  if (dlon >= 180) dlon -= 360;
+  else if (dlon < -180) dlon += 360;
+  return a * Math::degree() *
+    hypot((lat1 - lat0) * hlat, dlon * hlon);
+}
+
+int usage(int retval) {
+  ( retval ? std::cerr : std::cout ) <<
+"ProjTest latlonfile projfile\n\
+\n\
+Checks projections against NGS GoldData files\n";
+  return retval;
+}
+
+int main(int argc, char* argv[]) {
+  using namespace GeographicLib;
+  typedef Math::real real;
+
+  if (argc != 3)
+    return usage(1);
+  try {
+    real eps = 2*std::numeric_limits<real>::epsilon();
+
+    enum { undef = 0, tm = 1, ps = 2, lcc = 3 };
+    int m = 0;
+    std::string geof(argv[++m]);
+    DataFile geo(geof);
+    std::string projf(argv[++m]);
+    DataFile proj(projf);
+    if (geo.coords() != "Geodetic")
+      throw std::out_of_range("Unsupported coordinates " + geo.coords());
+    real a, f;
+    if (geo.datum() == "WGE") {
+      a = Constants::WGS84_a();
+      f = Constants::WGS84_f();
+    } else if (geo.datum() == "Test_sphere") {
+      a = 20000000/Math::pi();
+      f = 0;
+    } else if (geo.datum() == "Test_SRMmax") {
+      a = 6400000;
+      f = 1/real(150);
+    } else
+      throw std::out_of_range("Unsupported datum " + geo.datum());
+    if (proj.datum() != geo.datum())
+      throw std::out_of_range("Datum mismatch " +
+                              geo.datum() + " " + proj.datum());
+    real fe, fn, lat0, lat1, lat2, lon0, k1;
+    int type = undef;
+    if (proj.proj() == "Lambert Conformal Conic (1 parallel)") {
+      fe = proj.fe();
+      fn = proj.fn();
+      k1 = proj.k0();
+      lon0 = proj.lon0();
+      lat1 = lat2 = lat0 = proj.lat0();
+      type = lcc;
+    } else if (proj.proj() == "Lambert Conformal Conic (2 parallel)") {
+      fe = proj.fe();
+      fn = proj.fn();
+      k1 = 1;
+      lon0 = proj.lon0();
+      lat0 = proj.lat0();
+      lat1 = proj.lat1();
+      lat2 = proj.lat2();
+      type = lcc;
+    } else if (proj.proj() == "Mercator") {
+      using std::isfinite;
+      fe = proj.fe();
+      fn = proj.fn();
+      k1 = proj.k0();
+      if (!isfinite(k1))
+        k1 = 1;
+      lon0 = proj.lon0();
+      lat0 = 0;
+      lat1 = proj.latts();
+      if (!isfinite(lat1))
+        lat1 = 0;
+      lat2 = -lat1;
+      type = lcc;
+    } else if (proj.proj() == "Polar Stereographic") {
+      using std::isfinite;
+      fe = proj.fe();
+      fn = proj.fn();
+      lon0 = proj.lonfp();
+      lat0 = 90;
+      k1 = proj.k0();
+      if (!isfinite(k1))
+        k1 = 1;
+      real latts = proj.latts();
+      if (isfinite(latts)) {
+        if (latts < 0)
+          lat0 = -lat0;
+        LambertConformalConic tx(a, f, lat0, 1);
+        real x, y, gam, k;
+        tx.Forward(0, latts, 10, x, y, gam, k);
+        k1 = 1/k;
+      }
+      lat1 = lat2 = lat0;
+      type = lon0 == 0 ? ps : lcc;
+    } else if (proj.proj() == "Transverse Mercator") {
+      fe = proj.fe();
+      fn = proj.fn();
+      lon0 = proj.lon0();
+      k1 = proj.k0();
+      lat0 = lat1 = lat2 = proj.lat0();
+      type = tm;
+    } else
+      throw std::out_of_range("Unsupported projection " + proj.proj());
+    LambertConformalConic tx(a, f, lat1, lat2, k1);
+    PolarStereographic txa(a, f, k1);
+    TransverseMercator txb(a, f, k1);
+    TransverseMercatorExact txc(a, f == 0 ? real(0.1)/eps : f, k1);
+    std::cout << a << " 1/" << 1/f << " " << k1 << " " << type << "\n";
+    real x0, y0, gam, k;
+    if (type == lcc)
+      tx.Forward(lon0, lat0, lon0, x0, y0, gam, k);
+    else if (type == tm)
+      txb.Forward(lon0, lat0, lon0, x0, y0, gam, k);
+    else
+      x0 = y0 = 0;
+    real lata, lona, xa, ya;
+    unsigned count = 0;
+    real maxerrx = 0, maxerry = 0, maxerr = 0, maxerrk = 0, maxerrr = 0;
+    std::cout << std::fixed << std::setprecision(7);
+    while (geo.Next(lata, lona) && proj.Next(xa, ya)) {
+      using std::abs; using std::hypot;
+      ++count;
+      // Suppress bogus uninitialized warnings for lat and lon
+      real lat = 0, lon = 0, x, y, xx, yy;
+      x = xa - fe + x0;
+      y = ya - fn + y0;
+      xx = x/k1;
+      yy = y/k1;
+      switch (type) {
+      case lcc:
+        tx.Reverse(lon0, x, y, lat, lon, gam, k);
+        break;
+      case ps:
+        txa.Reverse(lat1 > 0, x, y, lat, lon, gam, k);
+        break;
+      case tm:
+        txb.Reverse(lon0, x, y, lat, lon, gam, k);
+        break;
+      default: break;           // To suppress warning
+      }
+      real errr = dist(a, f, lata, lona, lat, lon);
+      maxerrr = std::max(errr, maxerrr);
+      if (!(errr < 1e-6))
+        std::cout << "REV: "
+                  << lata << " "
+                  << lona << " "
+                  << xx << " "
+                  << yy << " "
+                  << errr << "\n";
+
+      switch (type) {
+      case lcc:
+        tx.Forward(lon0, lata, lona, x, y, gam, k);
+        break;
+      case ps:
+        txa.Forward(lat1 > 0, lata, lona, x, y, gam, k);
+        break;
+      case tm:
+        txb.Forward(lon0, lata, lona, x, y, gam, k);
+        break;
+      default: break;           // To suppress warning
+      }
+      x -= x0;
+      y -= y0;
+      x += fe;
+      y += fn;
+      real
+        errx = abs(x - xa),
+        erry = abs(y - ya),
+        err = hypot(errx, erry),
+        errk = err/std::max(real(1),k);
+      std::ostringstream sx, sxa, sy, sya;
+      sx << std::fixed << std::setprecision(6) << x;
+      sxa << std::fixed << std::setprecision(6) << xa;
+      sy << std::fixed << std::setprecision(6) << y;
+      sya << std::fixed << std::setprecision(6) << ya;
+      // if (sx.str() != sxa.str() || sy.str() != sya.str())
+      if (!(errk < 1e-6))
+        std::cout << "FOR: "
+                  << lata << " "
+                  << lona << " "
+                  << xx << " "
+                  << yy << " "
+                  << errk << "\n";
+      maxerrx = std::max(errx, maxerrx);
+      maxerry = std::max(erry, maxerry);
+      maxerr = std::max(err, maxerr);
+      maxerrk = std::max(errk, maxerrk);
+    }
+    std::cout << count << " records; maxerr "
+              << maxerrk << " " << maxerrr << "\n";
+  }
+  catch (const std::exception& e) {
+    std::cout << "ERROR: " << e.what() << "\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/TMTest.cpp

@@ -0,0 +1,227 @@
+/**
+ * \file TMTest.cpp
+ * \brief Command line utility for testing transverse Mercator projections
+ *
+ * Copyright (c) Charles Karney (2008-2022) <charles@karney.com> and licensed
+ * under the MIT/X11 License.  For more information, see
+ * https://geographiclib.sourceforge.io/
+ **********************************************************************/
+
+#include "GeographicLib/TransverseMercator.hpp"
+#include "GeographicLib/TransverseMercatorExact.hpp"
+#include "GeographicLib/Geodesic.hpp"
+#include "GeographicLib/Constants.hpp"
+#include <vector>
+#include <algorithm>
+#include <string>
+#include <limits>
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <iomanip>
+#include <stdexcept>
+
+using namespace GeographicLib;
+
+GeographicLib::Math::real
+dist(GeographicLib::Math::real a, GeographicLib::Math::real f,
+     GeographicLib::Math::extended lat0, GeographicLib::Math::extended lon0,
+     GeographicLib::Math::real lat1, GeographicLib::Math::real lon1) {
+  using namespace GeographicLib;
+  using std::cos; using std::sin; using std::sqrt; using std::hypot;
+  typedef Math::real real;
+  real
+    phi = real(lat0) * Math::degree(),
+    e2 = f * (2 - f),
+    sinphi = sin(phi),
+    n = 1/sqrt(1 - e2 * sinphi * sinphi),
+    // See Wikipedia article on latitude
+    hlon = cos(phi) * n,
+    hlat = (1 - e2) * n * n * n;
+  Math::extended dlon = Math::extended(lon1) - lon0;
+  if (dlon >= 180) dlon -= 360;
+  else if (dlon < -180) dlon += 360;
+  return a * Math::degree() *
+    hypot(real(Math::extended(lat1) - lat0) * hlat, real(dlon) * hlon);
+}
+
+int usage(int retval) {
+  ( retval ? std::cerr : std::cout ) <<
+"TMTest [-s] [-d]\n\
+\n\
+Read in TMcoords.dat on standard input and test TransverseMercatorExact\n\
+or (if -s is given) TransverseMercator.  If -d dump the error for each\n\
+point; otherwise summarise errors.  If -tf, perform a timing test of the\n\
+forward projection.  If -tr,  perform a timing test of the reverse\n\
+projection.\n";
+  return retval;
+}
+
+int main(int argc, char* argv[]) {
+  using namespace GeographicLib;
+  typedef Math::real real;
+  bool series = false;
+  bool dump = false;
+  bool timefor = false, timerev = false;
+  for (int m = 1; m < argc; ++m) {
+    std::string arg(argv[m]);
+    if (arg == "-s") {
+      series = true;
+    } else if (arg == "-d") {
+      dump = true;
+      timefor = false;
+      timerev = false;
+    } else if (arg == "-tf") {
+      dump = false;
+      timefor = true;
+      timerev = false;
+    } else if (arg == "-tr") {
+      dump = false;
+      timefor = false;
+      timerev = true;
+    } else
+      return usage(arg != "-h");
+  }
+
+  if (timefor || timerev) {
+    real s = 0;
+    int count = 0;
+    real dlat = real(0.015), dlon = real(0.015), dx = 2e3, dy = 2e3;
+    if (series) {
+      const TransverseMercator& tm = TransverseMercator::UTM();
+      if (timefor) {
+        real x, y, gam, k;
+        for (real lat = -80.0; lat <= 84.0; lat += dlat)
+          for (real lon = -3.0; lon <= 3.0; lon += dlon) {
+            tm.Forward(0.0, lat, lon, x, y, gam, k);
+            s += k;
+            ++count;
+          }
+      } else {
+        real lat, lon, gam, k;
+        for (real x = -400e3; x <= 400e3; x += dx)
+          for (real y = -9000e3; y <= 9500e3; y += dy) {
+            tm.Reverse(0.0, x, y, lat, lon, gam, k);
+            s += k;
+            ++count;
+          }
+      }
+    } else {
+      const TransverseMercatorExact tm(Constants::WGS84_a<real>(),
+                                       Constants::WGS84_f<real>(),
+                                       Constants::UTM_k0<real>(),
+                                       true);
+      if (timefor) {
+        real x, y, gam, k;
+        for (real lat = -80.0; lat <= 84.0; lat += dlat)
+          for (real lon = -3.0; lon <= 3.0; lon += dlon) {
+            tm.Forward(0.0, lat, lon, x, y, gam, k);
+            s += k;
+            ++count;
+          }
+      } else {
+        real lat, lon, gam, k;
+        for (real x = -400e3; x <= 400e3; x += dx)
+          for (real y = -9000e3; y <= 9500e3; y += dy) {
+            tm.Reverse(0.0, x, y, lat, lon, gam, k);
+            s += k;
+            ++count;
+          }
+      }
+    }
+    std::cout << count << " " << s << "\n";
+    return 0;
+  }
+
+  try {
+    real minlat = series ? 0 : (dump ? -100 : -15);
+    const unsigned nbins = 101;
+    std::vector<real> d(nbins);
+    std::vector<real> errv(nbins, 0);
+    std::vector<real> errvg(nbins, 0);
+    std::vector<real> errvk(nbins, 0);
+    real esterr = real(sizeof(real) == sizeof(double) ?
+                       (series ? 3e-9 : 8e-9) :
+                       (series ? 4e-12 : 4e-12));
+    for (unsigned i = 0; i < nbins; ++i)
+      d[i] = real(100e3 * i);
+    d[0] = 10e3;
+    d[nbins - 1] = 10001966;
+    const TransverseMercator& tm = TransverseMercator::UTM();
+    const TransverseMercatorExact tme(Constants::WGS84_a<real>(),
+                                      Constants::WGS84_f<real>(),
+                                      Constants::UTM_k0<real>(),
+                                      true);
+    real
+      a = series ? tm.EquatorialRadius() : tme.EquatorialRadius(),
+      f = series ? tm.Flattening() : tme.Flattening();
+    const Geodesic geod(a, f);
+    Math::extended lat0l, lon0l, x0l, y0l, gam0l, k0l;
+    while (std::cin >> lat0l >> lon0l >> x0l >> y0l >> gam0l >> k0l) {
+      using std::abs; using std::sin; using std::hypot;
+      real
+        lat0 = real(lat0l),
+        lon0 = real(lon0l),
+        x0 = real(x0l),
+        y0 = real(y0l),
+        gam0 = real(gam0l),
+        k0 = real(k0l);
+      if (lat0 < minlat)
+        continue;
+      real azi1, azi2, s12, m12;
+      real errf, errgf, errkf, errr, errgr, errkr;
+      geod.Inverse(std::max(lat0,real(0)), lon0, std::max(lat0,real(0)), -lon0,
+                   s12, azi1, azi2, m12);
+      s12 /= 2;
+      real lat, lon, x, y, gam, k;
+      if (series) {
+        tm.Forward(0, lat0, lon0, x, y, gam, k);
+      } else
+        tme.Forward(0, lat0, lon0, x, y, gam, k);
+      errf = real(hypot(Math::extended(x) - x0l,
+                        Math::extended(y) - y0l)) / k0;
+      errgf = real(abs(Math::extended(gam) - gam0));
+      errkf = real(abs(Math::extended(k) - k0));
+      if (series) {
+        tm.Reverse(0, x0, y0, lat, lon, gam, k);
+      } else
+        tme.Reverse(0, x0, y0, lat, lon, gam, k);
+      errr = dist(a, f, lat0l, lon0l, lat, lon);
+      errgr = real(abs(Math::extended(gam) - gam0));
+      errkr = real(abs(Math::extended(k) - k0));
+
+      real
+        err = std::max(errf, errr),
+        errg = std::max(errgf, errgr)
+        - esterr/(a * sin((90 - lat0) * Math::degree())
+                  * Math::degree()),
+        errk = std::max(errkf, errkr) / k0;
+      if (dump)
+        std::cout << std::fixed << std::setprecision(12)
+                  << lat0 << " " << lon0 << " "
+                  << std::scientific << std::setprecision(4)
+                  << errf << " " << errr << " "
+                  << errgf << " " << errgr << " "
+                  << errkf << " " << errkr << "\n";
+      else
+        for (unsigned i = 0; i < nbins; ++i) {
+          if (s12 <= d[i]) {
+            errv[i] = std::max(err, errv[i]);
+            errvg[i] = std::max(errg, errvg[i]);
+            errvk[i] = std::max(errk, errvk[i]);
+          }
+        }
+    }
+    if (!dump)
+      for (unsigned i = 0; i < nbins; ++i)
+        std::cout << int(d[i]/1000) << " "
+                  << errv[i] << " "
+                  << errvg[i] << " "
+                  << errvk[i] << "\n";
+  }
+  catch (const std::exception& e) {
+    std::cout << "ERROR: " << e.what() << "\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/dms-sample.txt

@@ -0,0 +1,29 @@
+1d2:3 0
+1:2:3 0
+1d2'3" 0
+1D2'3" 0
+1°2′3″ 0
+1º2‵3˝ 0
+1⁰2´3“ 0
+1˚2`3” 0
+1∘2‘3‟ 0
+1*2’3´´ 0
+1:2‛3`` 0
+1∘2‛3: 0
++1:2:3 0
++1d2'3" 0
++1D2'3" 0
++1°2′3″ 0
++1º2‵3˝ 0
++1⁰2´3“ 0
++1˚2`3” 0
++1∘2‘3‟ 0
++1*2’3´´ 0
++1:2‛3`` 0
+-1:2:3 0
+‐1d2'3" 0
+‑1D2'3" 0
+–1°2′3″ 0
+—1º2‵3˝ 0
+−1⁰2´3“ 0
+—1˚2`3” 0

libs/geographiclib/develop/intersect.cpp

@@ -0,0 +1,78 @@
+#include <iostream>
+#include <iomanip>
+#include <GeographicLib/Utility.hpp>
+#include <GeographicLib/Gnomonic.hpp>
+#include <GeographicLib/Geodesic.hpp>
+
+class vector3 {
+public:
+  GeographicLib::Math::real _x, _y, _z;
+  vector3(GeographicLib::Math::real x,
+          GeographicLib::Math::real y,
+          GeographicLib::Math::real z = 1) throw()
+    : _x(x)
+    , _y(y)
+    , _z(z) {}
+  vector3 cross(const vector3& b) const throw() {
+    return vector3(_y * b._z - _z * b._y,
+                   _z * b._x - _x * b._z,
+                   _x * b._y - _y * b._x);
+  }
+  void norm() throw() {
+    _x /= _z;
+    _y /= _z;
+    _z = 1;
+  }
+};
+
+int main() {
+  GeographicLib::Utility::set_digits();
+  GeographicLib::Math::real lata1, lona1, lata2, lona2;
+  GeographicLib::Math::real latb1, lonb1, latb2, lonb2;
+  std::cin >> lata1 >> lona1 >> lata2 >> lona2
+           >> latb1 >> lonb1 >> latb2 >> lonb2;
+  const GeographicLib::Geodesic& geod = GeographicLib::Geodesic::WGS84();
+  const GeographicLib::Gnomonic gn(geod);
+  GeographicLib::Math::real
+    lat0 = (lata1 + lata2 + latb1 + latb2)/4,
+    // Possibly need to deal with longitudes wrapping around
+    lon0 = (lona1 + lona2 + lonb1 + lonb2)/4;
+  std::cout << std::setprecision(16);
+  std::cout << "Initial guess " << lat0 << " " << lon0 << "\n";
+  GeographicLib::Math::real s120 = 1, s121;
+  for (int i = 0; i < 20; ++i) {
+    GeographicLib::Math::real xa1, ya1, xa2, ya2;
+    GeographicLib::Math::real xb1, yb1, xb2, yb2;
+    gn.Forward(lat0, lon0, lata1, lona1, xa1, ya1);
+    gn.Forward(lat0, lon0, lata2, lona2, xa2, ya2);
+    gn.Forward(lat0, lon0, latb1, lonb1, xb1, yb1);
+    gn.Forward(lat0, lon0, latb2, lonb2, xb2, yb2);
+    // See Hartley and Zisserman, Multiple View Geometry, Sec. 2.2.1
+    vector3 va1(xa1, ya1); vector3 va2(xa2, ya2);
+    vector3 vb1(xb1, yb1); vector3 vb2(xb2, yb2);
+    // la is homogeneous representation of line A1,A2
+    // lb is homogeneous representation of line B1,B2
+    vector3 la = va1.cross(va2);
+    vector3 lb = vb1.cross(vb2);
+    // p0 is homogeneous representation of intersection of la and lb
+    vector3 p0 = la.cross(lb);
+    p0.norm();
+    GeographicLib::Math::real lat1, lon1;
+    gn.Reverse(lat0, lon0, p0._x, p0._y, lat1, lon1);
+    geod.Inverse(lat1, lon1, lat0, lon1, s121);
+    std::cout << "Increment " << s121 << " ratio " << s121/s120 << "\n";
+    lat0 = lat1;
+    lon0 = lon1;
+    s120 = s121;
+  }
+  std::cout << "Final result " << lat0 << " " << lon0 << "\n";
+  GeographicLib::Math::real azi1, azi2;
+  geod.Inverse(lata1, lona1, lat0, lon0, azi1, azi2);
+  std::cout << "Azimuths on line A " << azi2 << " ";
+  geod.Inverse(lat0, lon0, lata2, lona2, azi1, azi2);
+  std::cout << azi1 << "\n";
+  geod.Inverse(latb1, lonb1, lat0, lon0, azi1, azi2);
+  std::cout << "Azimuths on line B " << azi2 << " ";
+  geod.Inverse(lat0, lon0, latb2, lonb2, azi1, azi2);
+  std::cout << azi1 << "\n";
+}

libs/geographiclib/develop/listallfiles.sh

@@ -0,0 +1,10 @@
+#! /bin/sh
+cd $HOME
+for d in geographiclib \
+             git/geographiclib-{c,fortran,java,js,octave,python,doc}; do
+    (
+        cd $d
+        git ls-tree -r HEAD --name-only |
+            sed -e "s%^%$HOME/$d/%"
+    )
+done

libs/geographiclib/develop/matlab-toolbox-config.sh

@@ -0,0 +1,81 @@
+#! /bin/sh
+(
+    cat <<EOF
+<deployment-project plugin="plugin.toolbox" plugin-version="1.0">
+  <configuration build-checksum="" file="$ROOT/geographiclib.prj" location="$ROOT" name="geographiclib" target="target.toolbox" target-name="Package Toolbox">
+    <param.appname>geographiclib</param.appname>
+    <param.authnamewatermark>Charles Karney</param.authnamewatermark>
+    <param.email>charles@karney.com</param.email>
+    <param.company />
+    <param.summary>MATLAB implementations of a subset of the C++ library, GeographicLib</param.summary>
+    <param.description>GeographicLib toolbox
+Version $VERSION $DATE
+
+EOF
+    cat $ROOT/geographiclib-blurb.txt
+    cat <<EOF
+</param.description>
+    <param.screenshot>\${PROJECT_ROOT}/geodesic.png</param.screenshot>
+    <param.version>$VERSION</param.version>
+    <param.output>\${PROJECT_ROOT}/geographiclib.mltbx</param.output>
+    <param.products.name>
+      <item>MATLAB</item>
+    </param.products.name>
+    <param.products.id>
+      <item>1</item>
+    </param.products.id>
+    <param.products.version>
+      <item>7.9</item>
+    </param.products.version>
+    <param.platforms />
+    <param.guid />
+    <param.exclude.filters />
+    <param.examples>&lt;?xml version="1.0" encoding="utf-8"?&gt;
+&lt;examples/&gt;</param.examples>
+    <param.demosxml />
+    <param.apps />
+    <param.docs />
+    <unset>
+      <param.company />
+      <param.output />
+      <param.platforms />
+      <param.exclude.filters />
+      <param.demosxml />
+      <param.apps />
+      <param.docs />
+    </unset>
+    <fileset.rootdir>
+      <file>\${PROJECT_ROOT}/geographiclib</file>
+    </fileset.rootdir>
+    <fileset.rootfiles>
+      <file>${PROJECT_ROOT}/geographiclib</file>
+    </fileset.rootfiles>
+    <fileset.depfun.included />
+    <fileset.depfun.excluded />
+    <fileset.package />
+    <build-deliverables>
+      <file location="$ROOT" name="geographiclib.mltbx" optional="false">$ROOT/geographiclib.mltbx</file>
+    </build-deliverables>
+    <workflow />
+    <matlab />
+    <platform>
+      <unix>true</unix>
+      <mac>false</mac>
+      <windows>false</windows>
+      <win2k>false</win2k>
+      <winxp>false</winxp>
+      <vista>false</vista>
+      <linux>true</linux>
+      <solaris>false</solaris>
+      <osver />
+      <os32>false</os32>
+      <os64>true</os64>
+      <arch>glnxa64</arch>
+      <matlab>true</matlab>
+    </platform>
+  </configuration>
+</deployment-project>
+EOF
+) > $ROOT/geographiclib.prj
+
+exit

libs/geographiclib/develop/reformat.cpp

@@ -0,0 +1,472 @@
+#include <iostream>
+#include <sstream>
+#include <fstream>
+#include <vector>
+#include <stdexcept>
+#include <algorithm>
+#include <cmath>
+
+int cosind(int N, int /*M*/, int n, int m) {
+  return m * N - m * (m - 1) / 2 + n;
+}
+
+int sinind(int N, int M, int n, int m) {
+  return cosind(N, M, n, m) - (N + 1);
+}
+
+void storecos(std::vector<double>& C, int N, int M,
+              double v, int n, int m) {
+  if (v == 0)
+    return;
+  if (n < 0 || n > N || m < 0 || m > M)
+    throw std::runtime_error("Invalid coefficient");
+  C[cosind(N, M, n, m)] = v;
+}
+
+void storesin(std::vector<double>& S, int N, int M,
+              double v, int n, int m) {
+  if (v == 0)
+    return;
+  if (n < 0 || n > N || m <= 0 || m > M)
+    throw std::runtime_error("Invalid coefficient");
+  S[sinind(N, M, n, m)] = v;
+}
+
+void storecoeff(std::ostream& str,
+                const std::vector<double>& C, const std::vector<double>& S,
+                int M, int N) {
+  int K = (M + 1) * (2*N - M + 2) / 2;
+  int ind[2] = {N, M};
+  str.write(reinterpret_cast<const char*>(ind), 2 * sizeof(int));
+  str.write(reinterpret_cast<const char*>(&C[0]), K * sizeof(double));
+  str.write(reinterpret_cast<const char*>(&S[0]), (K-N-1) * sizeof(double));
+}
+
+int main(int argc, char* argv[]) {
+  if (argc != 3) {
+    std::cerr << "Usage: reformat model outfile\n";
+    return 1;
+  }
+  std::string model = argv[1];
+  std::string file = argv[2];
+  try {
+    // model = wmm2010
+    //http://ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2010_Sph_Windows_Linux.zip
+    // unpack
+    // coefficients are in EMM_Sph_Windows_Linux/WMM2010.COF
+    // N=739, M=718
+    std::ofstream fout(file.c_str(), std::ios::binary);
+    if (model == "emm2010") {
+      // download
+      // http://ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2010_Sph_Windows_Linux.zip
+      // unpack
+      // coefficients are in EMM_Sph_Windows_Linux/EMM2010.COF
+      // and EMM_Sph_Windows_Linux/EMM2010SV.COF
+      std::string id = "EMM2010A";
+      fout.write(id.c_str(), 8);
+      for (int i = 0; i < 2; ++i) {
+        std::string filename = i == 0 ? "EMM_Sph_Windows_Linux/EMM2010.COF" :
+          "EMM_Sph_Windows_Linux/EMM2010SV.COF";
+        int
+          N = i == 0 ? 739 : 16,
+          M = i == 0 ? 718 : 16,
+          K = (M + 1) * (2*N - M + 2) / 2;
+        std::vector<double> C(K, 0.0);
+        std::vector<double> S(K - (N + 1), 0.0);
+        std::ifstream fin(filename.c_str());
+        std::string ss;
+        if (i == 0) std::getline(fin, ss); // Skip first line
+        while (std::getline(fin, ss)) {
+          int n, m;
+          double c, s;
+          std::istringstream is(ss);
+          if (!(is >> n >> m >> c >> s))
+            throw std::runtime_error("Short read");
+          storecos(C, N, M, c, n, m);
+          storesin(S, N, M, s, n, m);
+        }
+        storecoeff(fout, C, S, M, N);
+      }
+    } else if (model == "emm2015") {
+      // download
+      //http://www.ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2015_Sph_Linux.zip
+      // unpack
+      //
+      // * The only coefficients needed are EMM20{00-15}.COF and EMM2015SV.COF.
+      //
+      // * The other SV files can be ignored because the time dependence can be
+      //   obtained using linear interpolation for dates < 2015 (or
+      //   extrapolation for dates < 2000).
+      //
+      // * The time varying part of the field is of degree 15.  This
+      //   constitutes all the coefficients in EMM20{00-14}.COF and
+      //   EMM2015SV.COF and a subset of the coefficients in EMM2015.COF.
+      //
+      // * To this should be added a time independent short wavelength field
+      //   which is given by the degree > 15 terms in EMM2015.COF.
+      //
+      // * These time independent terms are of degree 729 and order 718.  There
+      //   are higher degree and order coefficients listed in the file, but
+      //   these are zero.
+      //
+      // * The EMM2015 coefficients as used by GeographicLib compress much
+      //   better than those for EMM2010 (660 kB instead of 3700 kB).
+      //   Presumably this is because the EMM2015 are only given with 4 decimal
+      //   digits.
+      //
+      // * The GeographicLib implementation produces the same results as listed
+      //   in EMM2015_TEST_VALUES.txt
+
+      std::string id = "EMM2015A";
+      fout.write(id.c_str(), 8);
+      for (int i = 0; i <= 17; ++i) {
+        std::string filename;
+        {
+          std::ostringstream os;
+          os << "EMM2015_linux/EMM" << (2000 + std::min(15, i))
+             << (i == 16 ? "SV" : "") << ".COF";
+          filename = os.str();
+        }
+        int
+          N = i == 17 ? 729 : 15,
+          M = i == 17 ? 718 : 15,
+          K = (M + 1) * (2*N - M + 2) / 2;
+        std::vector<double> C(K, 0.0);
+        std::vector<double> S(K - (N + 1), 0.0);
+        std::ifstream fin(filename.c_str());
+        std::string ss;
+        if (i != 16) std::getline(fin, ss); // Skip first line
+        while (std::getline(fin, ss)) {
+          int n, m;
+          double c, s;
+          std::istringstream is(ss);
+          if (!(is >> n >> m >> c >> s))
+            throw std::runtime_error("Short read " + filename + ": " + ss);
+          if (i == 15 && n > 15)
+            continue;
+          if (i == 17 && n <= 15)
+            continue;
+          storecos(C, N, M, c, n, m);
+          storesin(S, N, M, s, n, m);
+        }
+        storecoeff(fout, C, S, M, N);
+      }
+    } else if (model == "emm2017") {
+      // download
+      //https://www.ngdc.noaa.gov/geomag/EMM/data/geomag/EMM2017_Sph_Linux.zip
+      // unpack
+      //
+      // * The only coefficients needed are EMM20{00-17}.COF and EMM2017SV.COF.
+      //
+      // * The other SV files can be ignored because the time dependence can be
+      //   obtained using linear interpolation for dates < 2017 (or
+      //   extrapolation for dates < 2000).
+      //
+      // * The time varying part of the field is of degree 15.  This
+      //   constitutes all the coefficients in EMM20{00-16}.COF and
+      //   EMM2017SV.COF and a subset of the coefficients in EMM2017.COF.
+      //
+      // * To this should be added a time independent short wavelength field
+      //   which is given by the degree > 15 terms in EMM2017.COF.
+      //
+      // * These time independent terms are of degree 790 and order 790.
+      //
+      // * The GeographicLib implementation produces the same results as listed
+      //   in EMM2017TestValues.txt
+
+      std::string id = "EMM2017A";
+      fout.write(id.c_str(), 8);
+      int maxy = 17;
+      for (int i = 0; i <= 19; ++i) {
+        // i = maxy for low res components for 2000+maxy
+        // i = maxy+1 for SV at 2000+maxy
+        // i = maxy+2 for high res components for 2000+maxy
+        std::string filename;
+        {
+          std::ostringstream os;
+          os << "EMM2017_Linux/EMM" << (2000 + std::min(maxy, i))
+             << (i == maxy+1 ? "SV" : "") << ".COF";
+          filename = os.str();
+        }
+        int
+          N = i == maxy+2 ? 790 : 15,
+          M = i == maxy+2 ? 790 : 15,
+          K = (M + 1) * (2*N - M + 2) / 2;
+        std::vector<double> C(K, 0.0);
+        std::vector<double> S(K - (N + 1), 0.0);
+        std::ifstream fin(filename.c_str());
+        std::string ss;
+        if (i != maxy+1) std::getline(fin, ss); // Skip first line
+        while (std::getline(fin, ss)) {
+          int n, m;
+          double c, s;
+          std::istringstream is(ss);
+          if (!(is >> n >> m >> c >> s))
+            throw std::runtime_error("Short read " + filename + ": " + ss);
+          if (i == maxy && n > 15)
+            continue;
+          if (i == maxy+2 && n <= 15)
+            continue;
+          storecos(C, N, M, c, n, m);
+          storesin(S, N, M, s, n, m);
+        }
+        storecoeff(fout, C, S, M, N);
+      }
+    } else if (model == "wmm2010" || model == "igrf11") {
+      // Download
+      // http://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
+      // unpack
+      // wmm2010 coefficients are in geomag70_linux/WMM2010.COF
+      // igrf11 coefficients are in geomag70_linux/IGRF11.COF
+      std::string id = model == "wmm2010" ? "WMM2010A" : "IGRF11-A";
+      fout.write(id.c_str(), 8);
+      std::string filename = model == "wmm2010" ? "geomag70_linux/WMM2010.COF"
+        : "geomag70_linux-2010/IGRF11.COF";
+      std::ifstream fin(filename.c_str());
+      std::string ss;
+      bool start = true;
+      std::getline(fin, ss);
+      std::vector<double> C;
+      std::vector<double> S;
+      std::vector<double> C1;
+      std::vector<double> S1;
+      int N = 0, M = 0, N1 = 0, M1 = 0;
+      while (true) {
+        if (ss.size() == 0)
+          break;
+        std::istringstream is(ss);
+        int n, m;
+        double c, s, c1, s1;
+        if (start) {
+          std::string mm;
+          if (!(is >> mm >> mm >> N >> N1))
+            throw std::runtime_error("Short read on header");
+          M = N; M1 = N1;
+          int
+            K = (M + 1) * (2*N - M + 2) / 2,
+            K1 = (M1 + 1) * (2*N1 - M1 + 2) / 2;
+          C.resize(K);
+          S.resize(K - (N + 1));
+          C1.resize(K1);
+          S1.resize(K1 - (N1 + 1));
+          std::fill(C.begin(), C.end(), 0.0);
+          std::fill(S.begin(), S.end(), 0.0);
+          std::fill(C1.begin(), C1.end(), 0.0);
+          std::fill(S1.begin(), S1.end(), 0.0);
+          start = false;
+        } else {
+          if (!(is >> n >> m >> c >> s >> c1 >> s1))
+            throw std::runtime_error("Short read on data");
+          storecos(C, N, M, c, n, m);
+          storesin(S, N, M, s, n, m);
+          storecos(C1, N1, M1, c1, n, m);
+          storesin(S1, N1, M1, s1, n, m);
+        }
+        if (!std::getline(fin, ss))
+          ss = "";
+        if (ss.size() == 0 || (ss.size() >= 3 && ss.substr(0, 3) == "   ")) {
+          if (ss.size() > 0 && N1 != 0)
+            throw std::runtime_error("Secular coeffs not last");
+          if (ss.size() == 0 && N1 == 0)
+            throw std::runtime_error("No secular coeffs");
+          storecoeff(fout, C, S, M, N);
+          if (ss.size() == 0)
+            storecoeff(fout, C1, S1, M1, N1);
+          start = true;
+        }
+      }
+    } else if (model == "wmm2015" || model == "wmm2015v2" ||
+               model == "wmm2020" || model == "igrf12" ||
+               model == "igrf13") {
+      // Download WMM2015COF.zip
+      // http://ngdc.noaa.gov/geomag/WMM/WMM_coeff.shtml
+      // wmm2015 coefficients are in WMM2015COF/WMM.COF
+      //
+      // Download WMM2015v2COF.zip
+      // https://www.ngdc.noaa.gov/geomag/WMM/data/WMM2015/WMM2015v2COF.zip
+      // wmm2015v2 coefficients are in WMM2015v2COF/WMM.COF
+      //
+      // Download WMM2020COF.zip
+      // http://ngdc.noaa.gov/geomag/WMM/WMM_coeff.shtml
+      // wmm2020 coefficients are in WMM2020COF/WMM.COF
+      //
+      // igrf12 coefficients
+      // http://ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
+      // igrf12 coefficients are in geomag70_linux/IGRF12.COF
+      //
+      // igrf13 coefficients
+      // https://www.ngdc.noaa.gov/IAGA/vmod/geomag70_linux.tar.gz
+      // igrf13 coefficients are in geomag70_linux/IGRF13.COF
+      std::string id = model == "wmm2015" ? "WMM2015A" :
+        (model == "wmm2015v2" ?  "WMM2015B" :
+         (model == "wmm2020" ? "WMM2020A" :
+          (model == "igrf12" ? "IGRF12-A" : "IGRF13-A")));
+      fout.write(id.c_str(), 8);
+      std::string filename = model == "wmm2015" ? "WMM2015COF/WMM.COF"
+        : (model == "wmm2015v2" ? "WMM2015v2COF/WMM.COF"
+           : (model == "wmm2020" ? "WMM2020COF/WMM.COF"
+              : (model == "igrf12" ? "geomag70_linux-2015/IGRF12.COF"
+                 : "geomag70_linux-2020/IGRF13.COF")));
+      std::ifstream fin(filename.c_str());
+      std::string ss;
+      bool start = true;
+      std::getline(fin, ss);
+      std::vector<double> C;
+      std::vector<double> S;
+      std::vector<double> C1;
+      std::vector<double> S1;
+      int N = 0, M = 0, N1 = 0, M1 = 0;
+      while (true) {
+        if (ss.size() == 0)
+          break;
+        std::istringstream is(ss);
+        int n, m;
+        double c, s, c1, s1;
+        if (start) {
+          if (model == "wmm2015" || model == "wmm2015v2" ||
+              model == "wmm2020") {
+            N = 12; N1 = 12;
+          } else {
+            std::string mm;
+            if (!(is >> mm >> mm >> N >> N1))
+              throw std::runtime_error("Short read on header");
+          }
+          M = N; M1 = N1;
+          int
+            K = (M + 1) * (2*N - M + 2) / 2,
+            K1 = (M1 + 1) * (2*N1 - M1 + 2) / 2;
+          C.resize(K);
+          S.resize(K - (N + 1));
+          C1.resize(K1);
+          S1.resize(K1 - (N1 + 1));
+          std::fill(C.begin(), C.end(), 0.0);
+          std::fill(S.begin(), S.end(), 0.0);
+          std::fill(C1.begin(), C1.end(), 0.0);
+          std::fill(S1.begin(), S1.end(), 0.0);
+          start = false;
+        } else {
+          if (!(is >> n >> m >> c >> s >> c1 >> s1)) {
+            std::cerr << ss << "\n";
+            throw std::runtime_error("Short read on data");
+          }
+          storecos(C, N, M, c, n, m);
+          storesin(S, N, M, s, n, m);
+          storecos(C1, N1, M1, c1, n, m);
+          storesin(S1, N1, M1, s1, n, m);
+        }
+        if (!std::getline(fin, ss))
+          ss = "";
+        if (model == "wmm2015" || model == "wmm2015v2" ||
+            model == "wmm2020") {
+          if (ss.size() && ss[0] == '9')
+            ss = "";
+          if (ss.size() == 0) {
+            storecoeff(fout, C, S, M, N);
+            storecoeff(fout, C1, S1, M1, N1);
+          }
+        } else {
+          if (ss.size() == 0 || (ss.size() >= 3 && ss.substr(0, 3) == "   ")) {
+            if (ss.size() > 0 && N1 != 0)
+              throw std::runtime_error("Secular coeffs not last");
+            if (ss.size() == 0 && N1 == 0)
+              throw std::runtime_error("No secular coeffs");
+            storecoeff(fout, C, S, M, N);
+            if (ss.size() == 0)
+              storecoeff(fout, C1, S1, M1, N1);
+            start = true;
+          }
+        }
+      }
+    } else if (model == "egm2008" || model == "egm96" || model == "egm84"
+               || model == "wgs84") {
+      std::string id = model == "egm2008" ? "EGM2008A" :
+        (model == "egm96" ? "EGM1996A" :
+         (model == "egm84" ? "EGM1984A" : "WGS1984A"));
+      fout.write(id.c_str(), 8);
+      for (int i = 0; i < 2; ++i) {
+        std::string filename = "../gravity/";
+        filename +=
+          model == "egm2008" ?
+          (i == 0 ? "EGM2008_to2190_TideFree" : "Zeta-to-N_to2160_egm2008") :
+          (model == "egm96" ?
+           (i == 0 ? "EGM96" : "CORRCOEF") :
+           (model == "egm84" ?
+            (i == 0 ? "egm180.nor" : "zeta84") :
+            (i == 0 ? "wgs84.cof" : "zeta84")));
+        std::ifstream fin(filename.c_str());
+        if (!fin.good())
+            throw std::runtime_error("File not found");
+        int
+          N = model == "egm2008" ? (i == 0 ? 2190 : 2160) :
+          (model == "egm96" ? 360 :
+           (model == "egm84" ? (i == 0 ? 180 : -1) :
+            (i == 0 ? 20 : -1))),
+          M = model == "egm2008" ? (i == 0 ? 2159 : 2160) :
+          (model == "egm96" ? 360 :
+           (model == "egm84" ? (i == 0 ? 180 : -1) :
+            (i == 0 ? 0 : -1))),
+          K = (M + 1) * (2*N - M + 2) / 2;
+        std::vector<double> C(K, 0.0);
+        std::vector<double> S(K - (N + 1), 0.0);
+        std::string ss;
+        while (std::getline(fin, ss)) {
+          std::string::size_type p = 0;
+          while (true) {
+            p = ss.find_first_of('D', p);
+            if (p < ss.size())
+              ss[p] = 'E';
+            else
+              break;
+            ++p;
+          }
+          std::istringstream is(ss);
+          int n, m;
+          double c, s;
+          if (!(is >> n >> m >> c >> s))
+            throw std::runtime_error("Short read");
+          storecos(C, N, M, c, n, m);
+          storesin(S, N, M, s, n, m);
+        }
+        storecoeff(fout, C, S, M, N);
+      }
+    } else if (model == "dtm2006") {
+      std::string id = "DTM2006A";
+      fout.write(id.c_str(), 8);
+      int N = 2190, M = N;
+      int K = (M + 1) * (2*N - M + 2) / 2;
+      std::string filename
+        = "../gravity/Coeff_Height_and_Depth_to2190_DTM2006.0";
+      std::ifstream fin(filename.c_str());
+      std::vector<double> C(K, 0.0);
+      std::vector<double> S(K - (N + 1), 0.0);
+      std::string ss;
+      while (std::getline(fin, ss)) {
+        std::string::size_type p = 0;
+        while (true) {
+          p = ss.find_first_of('D', p);
+          if (p < ss.size())
+            ss[p] = 'E';
+          else
+            break;
+          ++p;
+        }
+        std::istringstream is(ss);
+        int n, m;
+        double c, s;
+        if (!(is >> n >> m >> c >> s))
+          throw std::runtime_error("Short read");
+        storecos(C, N, M, c, n, m);
+        storesin(S, N, M, s, n, m);
+      }
+      storecoeff(fout, C, S, M, N);
+    } else {
+      std::cerr << "UNKNOWN MODEL " << model << "\n";
+      return 1;
+    }
+  }
+  catch (const std::exception& e) {
+    std::cerr << "ERROR: " << e.what() << "\n";
+    return 1;
+  }
+  return 0;
+}

libs/geographiclib/develop/test-distribution.sh

@@ -0,0 +1,495 @@
+#! /bin/sh
+#
+# tar.gz and zip distrib files copied to $DEVELSOURCE
+# html documentation rsync'ed to $WEBDIST/htdocs/C++/$VERSION/
+#
+# Windows version ready to build in
+# $WINDOWSBUILD/GeographicLib-$VERSION/BUILD-vc10{,-x64}
+# after ./build installer is copied to
+# $DEVELSOURCE/GeographicLib-$VERSION-win{32,64}.exe
+#
+# Built version ready to install in /usr/local in
+# relc/GeographicLib-$VERSION/BUILD-system
+#
+# gita - check out from git, create package with cmake
+# gitb - check out from git, check various builds in the non-release tree
+# gitr - new release branch
+# relb - release package, build with autoconf
+# relc - release package, build with cmake
+# relx - cmake release package inventory
+# rely - autoconf release package inventory
+# instb - installed files, autoconf
+# instc - installed files, cmake
+
+set -e
+umask 0022
+
+# The following files contain version information:
+#   pom.xml
+#   CMakeLists.txt (PROJECT_VERSION_* LIBVERSION_*)
+#   NEWS
+#   configure.ac (AC_INIT, GEOGRAPHICLIB_VERSION_* LT_*)
+#   develop/test-distribution.sh
+#   doc/GeographicLib.dox.in (3 places)
+
+# maxima
+#   maxima/geodesic.mac
+
+START=`date +%s`
+DATE=`date +%F`
+VERSION=2.1.2
+SUFFIX=
+DISTVERSION=$VERSION$SUFFIX
+BRANCH=main
+TEMP=/home/scratch/geographiclib-dist
+if test `hostname` = petrel; then
+    DEVELSOURCE=$HOME/geographiclib
+    WINDEVELSOURCE=//datalake-pr-smb/vt-open/ckarney/geographiclib
+    WINDOWSBUILD=/var/tmp
+else
+    DEVELSOURCE=/u/geographiclib
+    WINDEVELSOURCE=//datalake-pr-smb/vt-open/ckarney/geographiclib
+    WINDOWSBUILD=/u/temp
+fi
+WINDOWSBUILDWIN=//datalake-pr-smb/vt-open/ckarney/temp
+GITSOURCE=file://$DEVELSOURCE
+WEBDIST=/home/ckarney/web/geographiclib-web
+mkdir -p $WEBDIST/htdocs/C++
+NUMCPUS=4
+HAVEINTEL=
+
+test -d $TEMP || mkdir $TEMP
+rm -rf $TEMP/*
+mkdir $TEMP/gita # Package creation via cmake
+mkdir $TEMP/gitb # Non release testing
+mkdir $TEMP/gitr # For release branch
+(cd $TEMP/gitr; git clone -b $BRANCH $GITSOURCE)
+(cd $TEMP/gita; git clone -b $BRANCH file://$TEMP/gitr/geographiclib)
+(cd $TEMP/gitb; git clone -b $BRANCH file://$TEMP/gitr/geographiclib)
+
+echo ==============================================================
+echo Make a source package in $TEMP/gita/geographiclib/BUILD
+
+cd $TEMP/gita/geographiclib
+cmake -S . -B BUILD
+(cd BUILD && make dist)
+# cp $TEMP/gita/geographiclib/BUILD/distrib/GeographicLib-$DISTVERSION.{zip,tar.gz} $DEVELSOURCE/data-distrib/distrib-C++/
+
+echo ==============================================================
+echo Unpack source package in $TEMP/rel bcx
+
+mkdir $TEMP/rel{b,c,x}
+tar xfpzC BUILD/distrib/GeographicLib-$DISTVERSION.tar.gz $TEMP/relb # Version for autoconf
+tar xfpzC BUILD/distrib/GeographicLib-$DISTVERSION.tar.gz $TEMP/relc # Version for cmake+mvn
+tar xfpzC BUILD/distrib/GeographicLib-$DISTVERSION.tar.gz $TEMP/relx # for listing
+
+echo ==============================================================
+echo Unpack devel cmake distribution in $TEMP/relx and list in $TEMP/files.x
+(
+    cd $TEMP/relx
+    find . -type f | sort -u > ../files.x
+)
+
+echo ==============================================================
+echo Make a release for Windows testing in $WINDOWSBUILD/GeographicLib-$VERSION
+rm -rf $WINDOWSBUILD/GeographicLib-$VERSION
+
+unzip -qq -d $WINDOWSBUILD BUILD/distrib/GeographicLib-$DISTVERSION.zip
+
+cat > $WINDOWSBUILD/GeographicLib-$VERSION/mvn-build <<'EOF'
+#! /bin/sh -exv
+unset GEOGRAPHICLIB_DATA
+# for v in 2019 2017 2015 2013 2012 2010; do
+for v in 2019 2017 2015; do
+  for a in 64 32; do
+    echo ========== maven $v-$a ==========
+    rm -rf c:/scratch/geog-mvn-$v-$a
+    mvn -Dcmake.compiler=vc$v -Dcmake.arch=$a \
+      -Dcmake.project.bin.directory=c:/scratch/geog-mvn-$v-$a install
+  done
+done
+EOF
+chmod +x $WINDOWSBUILD/GeographicLib-$VERSION/mvn-build
+cp pom.xml $WINDOWSBUILD/GeographicLib-$VERSION/
+
+# for ver in 10 11 12 14 15 16; do
+for ver in 14 15 16; do
+    for arch in win32 x64; do
+        pkg=vc$ver-$arch
+        gen="Visual Studio $ver"
+        installer=
+        # N.B. update CPACK_NSIS_INSTALL_ROOT in CMakeLists.txt and
+        # update documentation examples if VS version for binary
+        # installer changes.
+        test "$ver" = 14 && installer=y
+        (
+            echo "#! /bin/sh -exv"
+            echo echo ========== cmake $pkg ==========
+            echo b=c:/scratch/geog-$pkg
+            echo rm -rf \$b \$bc //datalake-pr-smb/vt-open/ckarney/pkg-$pkg/GeographicLib-$VERSION/\*
+            echo 'unset GEOGRAPHICLIB_DATA'
+            echo cmake -G \"$gen\" -A $arch -D BUILD_BOTH_LIBS=ON -D CMAKE_INSTALL_PREFIX=//datalake-pr-smb/vt-open/ckarney/pkg-$pkg/GeographicLib-$VERSION -D PACKAGE_DEBUG_LIBS=ON -D CONVERT_WARNINGS_TO_ERRORS=ON -D EXAMPLEDIR= -S . -B \$b
+            echo cmake --build \$b --config Debug   --target ALL_BUILD
+            echo cmake --build \$b --config Debug   --target RUN_TESTS
+            echo cmake --build \$b --config Debug   --target INSTALL
+            echo cmake --build \$b --config Release --target ALL_BUILD
+            echo cmake --build \$b --config Release --target exampleprograms
+            echo cmake --build \$b --config Release --target RUN_TESTS
+            echo cmake --build \$b --config Release --target INSTALL
+            echo cmake --build \$b --config Release --target PACKAGE
+            test "$installer" &&
+                echo cp \$b/GeographicLib-$DISTVERSION-*.exe $WINDEVELSOURCE/ ||
+                    true
+        ) > $WINDOWSBUILD/GeographicLib-$VERSION/build-$pkg
+        chmod +x $WINDOWSBUILD/GeographicLib-$VERSION/build-$pkg
+    done
+done
+cat > $WINDOWSBUILD/GeographicLib-$VERSION/test-all <<'EOF'
+#! /bin/sh
+(
+    for d in build-*; do
+        ./$d
+    done
+    ./mvn-build
+) >& build.log
+EOF
+chmod +x $WINDOWSBUILD/GeographicLib-$VERSION/test-all
+
+echo ==============================================================
+echo Set up release branch in $TEMP/gitr/geographiclib
+
+cd $TEMP/gitr/geographiclib
+git checkout release
+git config user.email charles@karney.com
+git ls-files | sort > ../files.old
+(
+    cd $TEMP/gitb/geographiclib
+    git ls-files | sort
+) > ../files.current
+cut -f3- -d/ $TEMP/files.x | sort > ../files.new
+(
+    cd ..
+    comm -23 files.old files.new | grep -v gitattributes > files.delete || true
+    comm -23 files.new files.current > files.relonly
+    comm -12 files.new files.current > files.main
+    comm -13 files.delete files.new > files.add
+    comm -13 files.old files.new > files.add
+    comm -12 files.old files.new > files.common
+)
+
+(
+    S=$TEMP/relx/GeographicLib-$VERSION
+    C=$TEMP/gitb/geographiclib
+    test -s ../files.delete && xargs git rm -f < ../files.delete
+    cat ../files.main ../files.relonly |
+        sed -e 's%[^/][^/]*$%%' -e 's%/$%%' | sort -u | grep -v '^$' |
+        while read d; do
+            test -d $d || mkdir -p $d
+        done
+    while read f; do
+        if cmp $S/$f $f >& /dev/null; then :
+        else
+            cp -p $S/$f $f
+            git add $f
+        fi
+    done < ../files.relonly
+    while read f; do
+        test -d `dirname $f` || mkdir `dirname $f`
+        git checkout $BRANCH $f
+        if cmp $S/$f $f >& /dev/null; then :
+        else
+            cp -p $S/$f $f
+            git add $f
+        fi
+    done < ../files.main
+    git checkout $BRANCH .gitattributes
+    for i in 1 2 3 4 5; do
+        find -type d -empty | xargs -r rmdir
+    done
+)
+rm -rf GeographicLib-$VERSION
+rm -f java/.gitignore
+for ((i=0; i<7; ++i)); do
+    find * -type d -empty | xargs -r rmdir
+done
+
+echo ==============================================================
+echo CMake build in $TEMP/relc/GeographicLib-$VERSION/BUILD install to $TEMP/instc
+cd $TEMP/relc/GeographicLib-$VERSION
+cmake -D BUILD_BOTH_LIBS=ON -D BUILD_DOCUMENTATION=ON -D USE_BOOST_FOR_EXAMPLES=ON -D CONVERT_WARNINGS_TO_ERRORS=ON -D CMAKE_INSTALL_PREFIX=$TEMP/instc -S . -B BUILD
+(
+    cd BUILD
+    make package_source
+    make -j$NUMCPUS all
+    make test
+    make exampleprograms
+    make install
+    # rsync -a --delete doc/html/ $WEBDIST/htdocs/C++/$VERSION/
+)
+
+echo ==============================================================
+echo List installed files fron cmake build in $TEMP/instc to $TEMP/files.c
+(
+    cd $TEMP/instc
+    find . -type f | sort -u > ../files.c
+)
+
+echo ==============================================================
+echo Make distribution from release tree with cmake
+
+cd $TEMP/relc/GeographicLib-$VERSION
+cmake -S . -B BUILD-dist
+(cd BUILD-dist && make package_source)
+
+echo ==============================================================
+echo Unpack release cmake distribution in $TEMP/relz and list in $TEMP/files.z
+mkdir $TEMP/relz
+tar xfpzC BUILD-dist/GeographicLib-$VERSION.tar.gz $TEMP/relz
+
+(
+    cd $TEMP/relz
+    find . -type f | sort -u > ../files.z
+)
+
+echo ==============================================================
+echo CMake build in $TEMP/relc/GeographicLib-$VERSION/BUILD-system install to /usr/local
+
+cmake -D BUILD_BOTH_LIBS=ON -D CONVERT_WARNINGS_TO_ERRORS=ON -S . -B BUILD-system
+(cd BUILD-system && make -j$NUMCPUS all&& make test)
+
+if test "$HAVEINTEL"; then
+    echo ==============================================================
+    echo CMake build for intel in $TEMP/relc/GeographicLib-$VERSION/BUILD-intel
+    env FC=ifort CC=icc CXX=icpc cmake -D BUILD_BOTH_LIBS=ON -D CONVERT_WARNINGS_TO_ERRORS=ON -S . -B BUILD-intel
+    (
+        cd BUILD-intel
+        make -j$NUMCPUS all
+        make test
+        make exampleprograms
+    )
+fi
+
+echo ==============================================================
+echo Check build with configure in $TEMP/relb/GeographicLib-$VERSION/BUILD-config to $TEMP/instb
+
+cd $TEMP/relb/GeographicLib-$VERSION
+mkdir BUILD-config
+cd BUILD-config
+../configure --prefix=$TEMP/instb
+make -j$NUMCPUS
+make install
+cd ..
+
+if test "$HAVEINTEL"; then
+echo ==============================================================
+echo Check build with configure + intell in $TEMP/relb/GeographicLib-$VERSION/BUILD-config-intel
+
+    mkdir BUILD-config-intel
+    cd BUILD-config-intel
+    env FC=ifort CC=icc CXX=icpc ../configure
+    make -j$NUMCPUS
+    cd ..
+fi
+
+echo ==============================================================
+echo Make source dist with autoconf $TEMP/relb/GeographicLib-$VERSION/BUILD-dist
+
+mkdir BUILD-dist
+cd BUILD-dist
+../configure
+make dist-gzip
+
+echo ==============================================================
+echo Unpack release autoconf distribution in $TEMP/rely and list in $TEMP/files.z
+
+mkdir $TEMP/rely
+tar xfpzC geographiclib-$VERSION.tar.gz $TEMP/rely
+mv $TEMP/rely/{geographiclib,GeographicLib}-$VERSION
+cd $TEMP/rely
+find . -type f | sort -u > ../files.y
+
+echo ==============================================================
+echo List installed files fron autoconf build in $TEMP/instb to $TEMP/files.b
+
+mv $TEMP/instb/share/doc/{geographiclib,GeographicLib}
+cd $TEMP/instb
+find . -type f | sort -u > ../files.b
+
+echo ==============================================================
+echo CMake build of devel tree in $TEMP/gitb/geographiclib/BUILD
+
+cd $TEMP/gitb/geographiclib
+cmake -D BUILD_BOTH_LIBS=ON -D BUILD_DOCUMENTATION=ON -D USE_BOOST_FOR_EXAMPLES=ON -D CONVERT_WARNINGS_TO_ERRORS=ON -S . -B BUILD
+(cd BUILD && make -j$NUMCPUS && make -j$NUMCPUS develprograms)
+
+cp $DEVELSOURCE/include/mpreal.h include/
+# Skip 4 for now because of various boost bugs
+for p in 1 3 5; do
+    echo ==============================================================
+    echo CMake build of devel tree at precision $p in $TEMP/gitb/geographiclib/BUILD-$p
+    mkdir BUILD-$p
+    cmake -D USE_BOOST_FOR_EXAMPLES=ON -D GEOGRAPHICLIB_PRECISION=$p -S . -B BUILD-$p
+    (
+        cd BUILD-$p
+        make -j$NUMCPUS all
+        if test $p -ne 1; then
+            make test
+        fi
+        make -j$NUMCPUS develprograms
+        make -j$NUMCPUS exampleprograms
+    )
+done
+
+echo ==============================================================
+echo Compile and run little test program
+cd $TEMP
+cat > testprogram.cpp <<EOF
+#include <iostream>
+#include <iomanip>
+
+#include <GeographicLib/Constants.hpp>
+#include <GeographicLib/DMS.hpp>
+#include <GeographicLib/LambertConformalConic.hpp>
+
+int main() {
+  using namespace GeographicLib;
+  double
+    // These are the constants for Pennsylvania South, EPSG:3364
+    // https://www.spatialreference.org/ref/epsg/3364/
+    a = Constants::WGS84_a(),   // major radius
+    f = 1/298.257222101,        // inverse flattening (GRS80)
+    lat1 = DMS::Decode(40,58),  // standard parallel 1
+    lat2 = DMS::Decode(39,56),  // standard parallel 2
+    k1 = 1,                     // scale on std parallels
+    lat0 =  DMS::Decode(39,20), // latitude of origin
+    lon0 = -DMS::Decode(77,45), // longitude of origin
+    fe = 600000,                // false easting
+    fn = 0;                     // false northing
+  LambertConformalConic PASouth(a, f, lat1, lat2, k1);
+  double x0, y0;
+  PASouth.Forward(lon0, lat0, lon0, x0, y0); // Transform origin point
+  x0 -= fe; y0 -= fn;           // Combine result with false origin
+
+  double lat = 39.95, lon = -75.17;    // Philadelphia
+  double x, y;
+  PASouth.Forward(lon0, lat, lon, x, y);
+  x -= x0; y -= y0;             // Philadelphia in PA South coordinates
+
+  std::cout << std::fixed << std::setprecision(3)
+            << x << " " << y << "\n";
+  return 0;
+}
+EOF
+for i in  b c; do
+    cp testprogram.cpp testprogram$i.cpp
+    g++ -c -g -O3 -I$TEMP/inst$i/include testprogram$i.cpp
+    g++ -g -o testprogram$i testprogram$i.o -Wl,-rpath=$TEMP/inst$i/lib \
+        -L$TEMP/inst$i/lib -lGeographicLib
+    ./testprogram$i
+done
+
+echo ==============================================================
+echo Verify library versions of cmake and autoconf builds are the same and other checks
+
+libversion=`find $TEMP/instc/lib -type f \
+-name 'libGeographicLib.so.*' -printf "%f" |
+sed 's/libGeographicLib\.so\.//'`
+test -f $TEMP/instb/lib/libGeographicLib.so.$libversion ||
+echo autoconf/cmake library so mismatch
+
+CONFIG_FILE=$TEMP/gitr/geographiclib/configure
+CONFIG_MAJOR=`grep ^GEOGRAPHICLIB_VERSION_MAJOR= $CONFIG_FILE | cut -f2 -d=`
+CONFIG_MINOR=`grep ^GEOGRAPHICLIB_VERSION_MINOR= $CONFIG_FILE | cut -f2 -d=`
+CONFIG_PATCH=`grep ^GEOGRAPHICLIB_VERSION_PATCH= $CONFIG_FILE | cut -f2 -d=`
+CONFIG_VERSIONA=`grep ^PACKAGE_VERSION= $CONFIG_FILE | cut -f2 -d= |
+cut -f2 -d\'`
+CONFIG_VERSION=$CONFIG_MAJOR.$CONFIG_MINOR
+test "$CONFIG_PATCH" = 0 || CONFIG_VERSION=$CONFIG_VERSION.$CONFIG_PATCH
+test "$CONFIG_VERSION"  = "$VERSION" || echo autoconf version number mismatch
+test "$CONFIG_VERSIONA" = "$VERSION" || echo autoconf version string mismatch
+
+cd $TEMP/relx/GeographicLib-$VERSION
+(
+    echo Files with trailing spaces:
+    find . -type f | egrep -v 'config\.guess|Makefile\.in|\.m4|\.png|\.gif' |
+        while read f; do
+            tr -d '\r' < $f | grep ' $' > /dev/null && echo $f || true
+        done
+    echo
+    echo Files with tabs:
+    find . -type f |
+        egrep -v '[Mm]akefile|\.html|\.m4|\.png|\.gif' |
+        egrep -v '\.sh|depcomp|install-sh|/config\.|configure$|compile|missing' |
+        xargs grep -l  '	' || true
+    echo
+    echo Files with multiple newlines:
+    find . -type f |
+        egrep -v \
+           '/Makefile\.in|\.1\.html|\.png|\.gif|/ltmain|/config|\.m4' |
+        while read f; do
+            tr 'X\n' 'xX' < $f | grep XXX > /dev/null && echo $f || true
+        done
+    echo
+    echo Files with no newline at end:
+    find . -type f |
+        egrep -v '\.png|\.gif' |
+        while read f; do
+            n=`tail -1 $f | wc -l`; test $n -eq 0 && echo $f || true
+        done
+    echo
+    echo Files with extra newlines at end:
+    find . -type f |
+        egrep -v '/configure|/ltmain.sh|\.png|\.gif|\.1\.html' |
+        while read f; do
+            n=`tail -1 $f | wc -w`; test $n -eq 0 && echo $f || true
+        done
+    echo
+) > $TEMP/badfiles.txt
+cat $TEMP/badfiles.txt
+cat > $TEMP/tasks.txt <<EOF
+# install built version
+sudo make -C $TEMP/relc/GeographicLib-$VERSION/BUILD-system install
+
+# commit and tag release branch
+cd $TEMP/gitr/geographiclib
+# Check .gitignore files!
+git add -A
+git commit -m "Version $VERSION ($DATE)"
+git tag -m "Version $VERSION ($DATE)" r$VERSION
+git push
+git push --tags
+
+# tag main branch
+cd $DEVELSOURCE
+git tag -m "Version $VERSION ($DATE)" v$VERSION
+git push --all
+git push --tags
+
+# Also to do
+# post release notices
+# set default download files
+# make -f makefile-admin distrib-{cgi,html}
+# update home brew
+#   dir = /usr/local/Homebrew/Library/Taps/homebrew/homebrew-core
+#   branch = geographiclib/$VERSION
+#   file = Formula/geographiclib.rb
+#   brew install --build-from-source geographiclib
+#   commit message = geographiclib $VERSION
+# update vcpkg git@github.com:microsoft/vcpkg
+#   dir = ports/geographiclib
+#   ./vcpkg install 'geographiclib[tools]'
+#   binaries in installed/x64-linux/tools/geographiclib
+#   libs in installed/x64-linux/{include,lib,debug/lib}
+#   ./vcpkg x-add-version geographiclib
+#   commit message = [geographiclib] Update to version $VERSION
+# update conda-forge
+#   url = git@github.com:conda-forge/geographiclib-cpp-feedstock
+#   conda build recipe
+# upload matlab packages
+# update binaries for cgi applications
+# trigger build on build-open
+EOF
+echo cat $TEMP/tasks.txt
+cat $TEMP/tasks.txt
+END=`date +%s`
+echo Elapsed time $((END-START)) secs