Vincenty Inverse Part 3

C++
Vincenty
Author

Sam Hutchins

Published

December 7, 2025

Well, today is Pearl Harbor Day! We should remember all those folks that died on that infamous day, and all the soldiers and civilians that allowed us to win the day during World War II.

Now, on to the subject of Vincenty’s formulae as started in this and this. I looked around recently and discovered formulae that does exactly what I am trying to do, so I decided not to reinvent the wheel anymore, and use what already exists. And, “what is that?” you ask.

“That” is a set of header files written for C++ that allows for calculating inverse distances and azimuths, among many other functions! Charles F. Karney wrote GeographicLib1, with the documentation found here. It also is written for many other programming languages as well, such as Python, Java, Fortran, etc. There is also a DMS header for conversion from and to degrees, minutes, seconds / decimal degrees, which I don’t go into here, as the input is in string format, and I have already written this file using double variables, so don’t wish to rewrite for strings.

Disclaimer: I am not a C++ programmer, so the defined program could be made better and more efficient, I’m sure! For me this is just a hobby.

So now, I am using these in my program to determine the inverse distance rather than Vincenty. Additionally, I am using Cmake for compiling, to increase the portability of the program across different operating systems (OSs), although in this case, I haven’t included different flags for different OSs. This example is done using the normal compilation procedure, (i.e., mkdir build && cd build, cmake .., make). This pulls in the required libraries so make will succeed. The CMakeLists.txt file looks like this:

Code 
cmake_minimum_required (VERSION 3.17.0)
project (karney)

find_package (GeographicLib REQUIRED)

if (NOT CMAKE_CONFIGURATION_TYPES AND NOT CMAKE_BUILD_TYPE)
  # Set a default build type for single-configuration cmake generators
  # if no build type is set.
  set (CMAKE_BUILD_TYPE "Release")
endif ()

add_executable (${PROJECT_NAME} karney.cpp)
target_link_libraries (${PROJECT_NAME} ${GeographicLib_LIBRARIES})

The inportant lines are the find_package() and target_link_libraries(). For Linux, this assumes you have installed the GeographicLib package. So, after performing the below steps, we end with the program, karney, which is executed from the command line.

Code 
mkdir build && cd build
cmake ..
make

The top of the karney.cpp files contains:

Code 
#include <iostream>
#include <GeographicLib/Geodesic.hpp>
#include <iomanip> // for setprecision()

using namespace std;
using namespace GeographicLib;

void menu();
double dms2dd(double D, double M, double S);
void dd2dms(double DD);
int inverse();
int direct();
void display(double s12, double azi1, double azi2);

Note the third statement where I include iomanip for setprecision() as it makes a smaller file than using bits/stdc++.h. I also define the functions used. Next, we have the main body of the program:

Code 
int main() {
    int choice {};
    do {
        //cout << "\033[2J\033[1;1H"; // clear screen
        menu();
        cout << "Enter your choice: ";
        if( ! (cin >> choice)) {
            cout << "Entry Error!" << endl;
            return -1;
        }
        switch (choice) {
            case 1:
                inverse();
                break;
            case 2:
                direct();
                break;
            case 0:
                exit(0);
                break;
            default:
                cout << "Invalid Choice." << endl;
        }
    } while (choice != 0);

    return 0;
}

Then we finish by including the actual functions as defined above:

Code 
void menu() {
    cout << "\033[2J\033[1;1H"; // clear screen
    cout << "==================================" << endl;
    cout << "=  Determine Geodesic Distances  =" << endl;
    cout << "==================================" << endl;
    cout << "Using GeographicLib by C.F.Karney.\n" << endl;
    cout << "NOTE: For location declination, use:" << endl;
    cout << "echo 2025-12-01 D:m:sN D:m:sW elev" << endl;
    cout << "  | MagneticField -r" << endl;
    cout << "\nSelect function:" << endl;
    cout << "1) Inverse." << endl;
    cout << "2) Direct." << endl;
    cout << "0) Quit." << endl;
}
int inverse() {
    const Geodesic& geod = Geodesic::WGS84();
    double lat1D{}, lat1M{}, lat1S{}, lon1D{}, lon1M{}, lon1S{};
    double lat2D{}, lat2M{}, lat2S{}, lon2D{}, lon2M{}, lon2S{};
    double lat1{}, lon1{}, lat2{}, lon2{}, s12{}, azi1{}, azi2{}, lat2a{}, lon2a{};
    cout << "--- INVERSE ---";
    cout << "\nSelect coordinate entry in decimal\ndegrees (DD.dddd), or D\u2070 M\' S\"." << endl;
    int sel{};
    cout << "\t1) Decimal degrees.\n\t2) DMS." << endl;
    cout << "Enter your choice: ";
    if( ! (cin >> sel)) {
            cout << "Entry Error!" << endl;
            return -1;
    }
    if(sel == 1) {
        cout << "Enter Lat1 (DD.dddd): "; // latitude of coordinate 1
        cin >> lat1;
        cout << "Enter Lon1 (DD.dddd): "; // longitude of coordinate 1
        cin >> lon1;
        cout << "Enter Lat2 (DD.dddd): "; // latitude of coordinate 2
        cin >> lat2;
        cout << "Enter Lon2 (DD.dddd): "; // longitude of coordinate 2
        cin >> lon2;
    } else {
        cout << "\nEnter Lat1 (DD mm ss): "; cin >> lat1D >> lat1M >> lat1S;
        cout << "Enter Lon1 (DD mm ss): "; cin >> lon1D >> lon1M >> lon1S;
        cout << "Enter Lat2 (DD mm ss): "; cin >> lat2D >> lat2M >> lat2S;
        cout << "Enter Lon2 (DD mm ss): "; cin >> lon2D >> lon2M >> lon2S;
        lat1 = dms2dd(lat1D, lat1M, lat1S); lon1 = dms2dd(lon1D, lon1M, lon1S);
        lat2 = dms2dd(lat2D, lat2M, lat2S); lon2 = dms2dd(lon2D, lon2M, lon2S);
    }
    geod.Inverse(lat1, lon1, lat2, lon2, s12, azi1, azi2); // determine distance and bearings
    display(s12, azi1, azi2);
    
    cout << "\nPress enter to continue. ";
    getchar(); // what a kludge,
    cin.get(); //   but works for pause.
    return 0; // Function type defines return value
}
int direct() {
    const Geodesic& geod = Geodesic::WGS84();
    double lat1D{}, lat1M{}, lat1S{}, lon1D{}, lon1M{}, lon1S{};
    double lat2D{}, lat2M{}, lat2S{}, lon2D{}, lon2M{}, lon2S{};
    double lat1{}, lon1{}, lat2{}, lon2{}, s12{}, azi1{}, azi2{}, lat2a{}, lon2a{};
    double azi1D{}, azi1M{}, azi1S{};
    cout << "--- DIRECT ---";
    cout << "\nSelect coordinate entry in decimal\ndegrees (DD.dddd), or D\u2070 M\' S\"." << endl;
    int sel{};
    cout << "\t1) Decimal degrees.\n\t2) DMS." << endl;
    cout << "Enter your choice: ";
    if( ! (cin >> sel)) {
            cout << "Entry Error!" << endl;
            return -1;
    }
    if(sel == 1) {
        cout << "Enter starting latitude (DD.dddd): "; cin >> lat1;
        cout << "Enter starting longitude (DD.dddd): "; cin >> lon1;
        cout << " Enter the bearing (DD.dddd): "; cin >> azi1;
    } else {
        cout << "\nEnter starting latitude (DD mm ss): "; cin >> lat1D >> lat1M >> lat1S;
        cout << "Enter starting longitude (DD mm ss): "; cin >> lon1D >> lon1M >> lon1S;
        cout << " Enter the bearing (DD mm ss): "; cin >> azi1D >> azi1M >> azi1S;
        lat1 = dms2dd(lat1D, lat1M, lat1S); lon1 = dms2dd(lon1D, lon1M, lon1S);
        azi1 = dms2dd(azi1D, azi1M, azi1S);
    }
    cout << " (meters = feet / 3.280833437)" << endl;
    cout << " Input the distance (meters): "; cin >> s12;
    geod.Direct(lat1, lon1, azi1, s12, lat2, lon2, azi2); // determine final coordinates
    //cout << lat2 << endl; // TEST
    //cout << lon2 << endl; // TEST
    cout << "Ending latitude: ";
    dd2dms(lat2); 
    cout << "Ending longitude: ";
    dd2dms(lon2);
    display(s12, azi1, azi2);
    
    cout << "\nPress enter to continue. ";
    getchar(); // what a kludge,
    cin.get(); //   but works for pause.
    return 0; // Function type defines return value
}
void display(double s12, double azi1, double azi2) {
    cout << "Distance:" << endl;
    cout << fixed << setprecision(3);
    cout << "\t" << s12 << " meters." << endl;
    cout << "\t" << s12 / 1000 << " km." << endl;
    cout << "\t" << s12 / 0.3048006 << " feet." << endl;
    cout << "\t" << s12 / 0.914401921 << " yards." << endl;
    cout << "\t" << s12 * 0.000621371 << " miles." << endl; 
    //cout << azi1 << endl; // TEST
    //cout << azi2 << endl; // TEST
    cout << "Initial bearing: ";
    if(azi1 < 0) azi1 += 360;
    dd2dms(azi1);
    cout << "Final bearing: ";
    if(azi2 < 0) azi2 += 360;
    dd2dms(azi2);
}
double dms2dd(double D, double M, double S) {
    if(D < 0) { M = -M; S = -S; } // makes all values negative if first negative.
    return D + (M/60) + (S/3600);
}
void dd2dms(double DD) {
    double Deg{}, Min{}, Sec{}, left1{}, left2{};
    left1 = DD * (3600); //cout << left1 << " ";
    Deg = trunc(left1/(3600)); //cout << Deg << " ";
    left2 = left1 - Deg * (3600); //cout << left2 << " ";
    Min = abs(trunc(left2/60.0)); //cout << Min << " ";
    Sec = abs(round(left2 - Min * 60)); //cout << Sec << endl;
    if(Sec >= 60) { Sec = 0; Min +=1; if(Min >=60) { Min = 0; if(Deg <0) { Deg -=1; } else {Deg +=1; } } }
    cout << fixed << setprecision(0) << Deg << "\u2070 " << Min << "\' " << fixed << setprecision(2) << Sec << "\""<< endl;
}

At the top of both inverse and direct functions, we define a class name geod to use those functions. This makes the code a bit cleaner and easier to understand. Also I have included in the menu() function, a short explanation snippet using one of the included utilities, MagneticField, which determines the up-to-date magnetic declination for the entered location. This is very handy if one wishes to set a particular bearing when doing real angles in the field.

One last thing. If we wished to have a CMakeLists.txt file that could be used on other OSs, we could have something that looked like this:

Code 
cmake_minimum_required( VERSION 3.6.2 )

# For a new project it is sufficient to change only its name in the following line
set( PROJECT_NAME karney  )

project( ${PROJECT_NAME} )

#set( CMAKE_BUILD_TYPE Debug )
set( CMAKE_BUILD_TYPE Release )


#[[ADD_DEFINITIONS(-g++ -O2 -fsigned-char -freg-struct-return -Wall -W -Wshadow -Wstrict-prototypes -Wpointer-arith -Wcast-qual -Winline -Werror)]]

if( WIN32 )
    set( CMAKE_CXX_FLAGS "/DWIN32 /D_WINDOWS /W3 /GR /EHsc /std:c++17 /D_UNICODE /DUNICODE" )
    set( CMAKE_CXX_FLAGS_DEBUG "/MDd /Zi /Ob0 /Od /RTC1 /std:c++17 /D_UNICODE /DUNICODE" )
    message( "Win settings chosen..." )
elseif( ${CMAKE_SYSTEM_NAME} STREQUAL "Darwin" )
    set( CMAKE_CXX_FLAGS "-std=c++17 -Wall" )
    set( CMAKE_CXX_FLAGS_DEBUG "-g -std=c++17 -Wall" )
    message( "Mac settings chosen..." )
elseif( UNIX )
    set( CMAKE_CXX_FLAGS "-std=c++17 -Wall" )
    set( CMAKE_CXX_FLAGS_DEBUG "-g -std=c++17 -Wall" )
    message( "Linux settings chosen..." )
endif()


# Inform CMake where the header files are
include_directories( include )


# Automatically add all *.cpp and *.h files to the project
file ( GLOB SOURCES "./src/*.cpp" "./include/*.h" )
add_executable( ${PROJECT_NAME} ${SOURCES} )


# Set the default project 
set_property( DIRECTORY ${CMAKE_CURRENT_SOURCE_DIR} PROPERTY VS_STARTUP_PROJECT ${PROJECT_NAME} )


message( "CMAKE_BUILD_TYPE is ${CMAKE_BUILD_TYPE}" )

For our purposes, this is way more complicated, so is shown for reference only. Note it will not compile without the find_package() and target_link_libraries() lines added. Have a wonderful day, and may God Bless!

Footnotes

  1. C. F. F. Karney, GeographicLib, Version 2.7 (2025-11-06), https://geographiclib.sourceforge.io/C++/2.7↩︎