Category: Exterior Ballistics

// CurvatureOfTheEarth.cpp : Defines the entry point for the application.
//

#include "pch.h"
#include "framework.h"
#include "CurvatureOfTheEarth.h"
#include "GreatCircleDistance.h"
#include "Vincenty.h"

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
WCHAR line[128];                                // general purpose buffer
std::wstring outputText;                        // output wide character text

// Forward declarations of functions included in this code module:
ATOM                MyRegisterClass(HINSTANCE hInstance);
BOOL                InitInstance(HINSTANCE, int);
LRESULT CALLBACK    WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    About(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK    TableDialog(HWND, UINT, WPARAM, LPARAM);

int APIENTRY wWinMain(_In_ HINSTANCE hInstance,
                     _In_opt_ HINSTANCE hPrevInstance,
                     _In_ LPWSTR    lpCmdLine,
                     _In_ int       nCmdShow)
{
    UNREFERENCED_PARAMETER(hPrevInstance);
    UNREFERENCED_PARAMETER(lpCmdLine);

    // TODO: Place code here.

    // Initialize global strings
    LoadStringW(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
    LoadStringW(hInstance, IDC_CURVATUREOFTHEEARTH, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

    // Perform application initialization:
    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

    HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_CURVATUREOFTHEEARTH));

    MSG msg;

    // Main message loop:
    while (GetMessage(&msg, nullptr, 0, 0))
    {
        if (!TranslateAccelerator(msg.hwnd, hAccelTable, &msg))
        {
            TranslateMessage(&msg);
            DispatchMessage(&msg);
        }
    }

    return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
    WNDCLASSEXW wcex = { 0 };

    wcex.cbSize = sizeof(WNDCLASSEX);

    wcex.style          = CS_HREDRAW | CS_VREDRAW;
    wcex.lpfnWndProc    = WndProc;
    wcex.cbClsExtra     = 0;
    wcex.cbWndExtra     = 0;
    wcex.hInstance      = hInstance;
    wcex.hIcon          = LoadIcon(hInstance, MAKEINTRESOURCE(IDI_CURVATUREOFTHEEARTH));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_CURVATUREOFTHEEARTH);
    wcex.lpszClassName  = szWindowClass;
    wcex.hIconSm        = LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

    return RegisterClassExW(&wcex);
}

//
//   FUNCTION: InitInstance(HINSTANCE, int)
//
//   PURPOSE: Saves instance handle and creates main window
//
//   COMMENTS:
//
//        In this function, we save the instance handle in a global variable and
//        create and display the main program window.
//
BOOL InitInstance(HINSTANCE hInstance, int nCmdShow)
{
   hInst = hInstance; // Store instance handle in our global variable

   HWND hWnd = CreateWindowW(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, nullptr, nullptr, hInstance, nullptr);

   if (!hWnd)
   {
      return FALSE;
   }

   ShowWindow(hWnd, nCmdShow);
   UpdateWindow(hWnd);

   return TRUE;
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_CREATE:
        DialogBox(hInst, MAKEINTRESOURCE(IDD_TABLE_DIALOG), hWnd, TableDialog);
        break;
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);
            // Parse the menu selections:
            switch (wmId)
            {
            case IDM_ABOUT:
                DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
                break;
            case IDM_EXIT:
                DestroyWindow(hWnd);
                break;
            default:
                return DefWindowProc(hWnd, message, wParam, lParam);
            }
        }
        break;
    case WM_PAINT:
        {
            PAINTSTRUCT ps;
            HDC hdc = BeginPaint(hWnd, &ps);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

// Message handler for about box.
INT_PTR CALLBACK About(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    }
    return (INT_PTR)FALSE;
}

INT_PTR CALLBACK TableDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    static double deltaTime = 0;
    static int height = 20, width = 80;
    static size_t delta = 0;
    static HFONT hFont = NULL;
    static HWND hEditMultiline = NULL;
    static GreatCircleDistance gcd;
    static Vincenty vincenty;

    switch (message)
    {
    case WM_INITDIALOG:
        hFont = CreateFont(
            -MulDiv(7, GetDeviceCaps(GetDC(hDlg), LOGPIXELSY), 72),
            0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
            DEFAULT_CHARSET, OUT_DEFAULT_PRECIS,
            CLIP_DEFAULT_PRECIS, DEFAULT_QUALITY,
            FIXED_PITCH | FF_MODERN,
            TEXT("Courier New")
        );

        hEditMultiline = CreateWindowEx(
            WS_EX_CLIENTEDGE,                       // Extended style for sunken border
            TEXT("EDIT"),                           // Class name
            TEXT(""),                               // Initial text (can be blank)
            WS_CHILD | WS_VISIBLE | WS_VSCROLL | ES_LEFT | ES_MULTILINE | ES_READONLY,
            10, 0, 10 * width, 30 * height,         // Position and size
            hDlg,                                   // Parent window handle
            (HMENU)IDC_EDIT_MULTILINE,              // Unique control ID
            hInst,                                  // Application instance
            NULL                                    // Extra parameter
        );

        SendMessage(hEditMultiline, WM_SETFONT, (WPARAM)hFont, TRUE);
        return (INT_PTR)TRUE;

    case WM_COMMAND:
        if (LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }

        if (LOWORD(wParam) == IDC_BUTTON_COMPUTE)
        {
            outputText += L"Range in Yards Versus Curvature of the Earth in Feet\r\n\r\n";
            outputText += L"Yards\t 0\t 100\t 200\t 300\t 400\t 500\t 600\t 700\t 800\t 900\r\n";
                
            for (double row = 1000.0; row <= 40000.0; row += 1000.0)
            {
                swprintf_s(line, L"%6.1lf\t", row);
                outputText += line;

                for (double col = 0.0; col <= 900.0; col += 100.0)
                {
                    double r = row + col;
                    double distance = 3.0 * gcd.Interpolate(r);

                    swprintf_s(line, L"%4.1lf\t", distance);
                    outputText += line;
                }

                outputText += L"\r\n";
            }

            outputText += L"Yards\t 0\t 100\t 200\t 300\t 400\t 500\t 600\t 700\t 800\t 900\r\n";
            SetWindowText(hEditMultiline, outputText.c_str());
            return (INT_PTR)TRUE;
        }
    }

    return (INT_PTR)FALSE;
}

#pragma once
#include "Vincenty.h"

struct PointRaR
{
    double ra;   // curvature of the Earth correction in yards
    double r;    // flat Earth distance (chord) in yards

    PointRaR(double ra, double r)
    {
        this->ra = ra;
        this->r = r;
    }
};

class GreatCircleDistance
{

private:

    Vincenty vincenty;
    std::vector<PointRaR> pts;

    bool binarySearch(double x, int& lt, int& rt);

public:

    friend int compare(
        const PointRaR& lt,
        const PointRaR& rt)
    {
        if (lt.r < rt.r)
            return -1;
        if (lt.r > rt.r)
            return +1;

        return 0;
    };

    // R is in yards, returns in yards
    // Curvature of the Earth correction

    double Interpolate(double R);

    // construction of the curvature of the Earth table

    GreatCircleDistance();
};

#pragma once

class Vincenty
{

public:

    static double Re; // Radius of Earth in meters
    double deltaSigma(
        double phi1, double lambda1,
        double phi2, double lambda2);
    double distance(
        double phi1, double lambda1,
        double phi2, double lambda2);
    double distanceConstantLambda(
        double phi1, double phi2);
    double distanceConstantXY(double z);
    double x(double phi, double lambda);
    double y(double phi, double lambda);
    double z(double phi);
    double phi(double z);
    double lambda(double x, double y);
};

#include "pch.h"
#include "GreatCircleDistance.h"

bool GreatCircleDistance::binarySearch(
    double x, int& lt, int& rt)
{
    int n = static_cast<int>(pts.size()), L = 0, M, R = n - 1;

Label10:

    if (x == pts[L].r)
    {
        lt = rt = L;
        return true;
    }

    if (x == pts[R].r)
    {
        lt = rt = R;
        return true;
    }

    if (x > pts[L].r && x < pts[R].r && R - L == 1)
    {
        lt = L;
        rt = R;
        return true;
    }

    if (x > pts[L].r && x < pts[R].r)
    {
        M = (L + R) / 2;

        if (x > pts[M].r)
        {
            L = M;
            goto Label10;
        }

        if (x < pts[M].r)
        {
            R = M;
            goto Label10;
        }
    }

    lt = rt = -1;
    return false;
}

double GreatCircleDistance::Interpolate(double R)
{
    int lt, rt;

    if (binarySearch(R, lt, rt))
    {
        double x0 = pts[lt].ra, x1 = pts[rt].ra;
        double y0 = pts[lt].r, y1 = pts[rt].r;
        double deltaX = x1 - x0, deltaY = y1 - y0;
        double ra = deltaX * (R - y0) / deltaY + x0;

        return ra;
    }

    return -1.0;
}

GreatCircleDistance::GreatCircleDistance()
{
    double deltaPhi = 0.000001, phi1 = 0.0, phi2 = deltaPhi, delta;
    double deltaRa, d0 = vincenty.z(phi1), d1, r, ra;
    int cnt = 0;

    pts.push_back(PointRaR(0.0, 0.0));

    while (cnt < 10000)
    {
        d1 = vincenty.z(phi2);

        if (d0 >= d1)
        {
            delta = d0 - d1;
            deltaRa = d0 * d0 - d1 * d1;
        }

        else
        {
            delta = d1 - d0;
            deltaRa = d1 * d1 - d0 * d0;
        }

        r = delta;
        ra = sqrt(deltaRa);
        ra = r >= ra ? r - ra : ra - r;
        pts.push_back(PointRaR(1.0936 * r, 1.0936 * ra));
        phi2 += deltaPhi;
        cnt++;
    }
}

#include "pch.h"
#include "Vincenty.h"

double Vincenty::Re = 6378137.0;	// radius of Earth in meters

double Vincenty::deltaSigma(
    double phi1, double lambda1,
    double phi2, double lambda2)
{
    double deltaPhi = phi1 - phi2, deltaLambda = lambda1 - lambda2;
    double cosPhi1 = cos(phi1), cosPhi2 = cos(phi2);
    double sinPhi1 = sin(phi1), sinPhi2 = sin(phi2);
    double cosDeltaLambda = cos(deltaLambda), sinDeltaLambda = sin(deltaLambda);
    double numer1 = cosPhi2 * sinDeltaLambda;
    double numer2 = cosPhi1 * sinPhi2 - sinPhi1 * cosPhi2 * cosDeltaLambda;
    double numer = sqrt(numer1 * numer1 + numer2 * numer2);
    double denom = sinPhi1 * sinPhi2 + cosPhi1 * cosPhi2 * cosDeltaLambda;

    return atan2(numer, denom);
}

double Vincenty::distance(
    double phi1, double lambda1,
    double phi2, double lambda2)
{
    return Re * deltaSigma(phi1, lambda1, phi2, lambda2);
}

double Vincenty::distanceConstantLambda(
    double phi1, double phi2)
{
    return Re * distance(phi1, 0.0, phi2, 0.0);
}

double Vincenty::distanceConstantXY(double z)
{
    return Re * distance(phi(0.0), 0.0, phi(z), 0.0);
}

double Vincenty::x(double phi, double lambda)
{
    return Re * sin(phi) * cos(lambda);
}

double Vincenty::y(double phi, double lambda)
{
    return Re * sin(phi) * sin(lambda);
}

double Vincenty::z(double phi)
{
    return Re * cos(phi);
}

double Vincenty::phi(double z)
{
    return acos(z / sqrt(Re));
}

double Vincenty::lambda(double x, double y)
{
    return asin(y / sqrt(x * x + y * y));
}