Blog Entry © Wednesday, June 4, 2025, Differential Cryptanalysis of DES by James Pate Williams, Jr.

Below are outputs from three different apps perform differential cryptanalysis on three, four, and six round DES algorithms. Three round is first etc.

E  = 007e0e80 680c0000
E* = bf02ac05 40520000
C' = 965d5b67
E = a0bff415 02f60000
E* = 8a6a5ebf 28aa0000
C' = 9c9c1f56
E = ef15068f 695f0000
E* = 05e9a2bf 56040000
C' = d575db2b

1 0 0 0 0 1 0 1 0 0 0 0 0 0 0 0
0 0 0 0 0 1 1 0 0 0 0 1 1 0 0 0
0 1 0 0 0 1 0 0 1 0 0 0 0 0 0 3
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 1

0 0 0 1 0 3 0 0 1 0 0 1 0 0 0 0
0 1 0 0 0 2 0 0 0 0 0 0 1 0 0 0
0 0 0 0 0 1 0 0 1 0 1 0 0 0 1 0
0 0 1 1 0 0 0 0 1 0 1 0 2 0 0 0

0 0 0 0 1 1 0 0 0 0 0 0 0 0 1 0
0 0 0 3 0 0 0 0 0 0 0 0 0 0 1 1
0 2 0 0 0 0 0 0 0 0 0 0 1 1 0 0
0 0 0 0 0 0 1 0 0 0 0 0 1 0 0 0

3 1 0 0 0 0 0 0 0 0 2 2 0 0 0 0
0 0 0 0 1 1 0 0 0 0 0 0 1 0 1 1
1 1 1 0 1 0 0 0 0 1 1 1 0 0 1 0
0 0 0 0 1 1 0 0 0 0 0 0 0 0 2 1

0 0 0 0 0 0 1 0 0 0 1 0 0 0 0 0
0 0 0 0 2 0 0 0 3 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 2 0 0 0 0 0 0 1 0 0 0 0 2 0

1 0 0 1 1 0 0 3 0 0 0 0 1 0 0 1
0 0 0 0 1 1 0 0 0 0 0 0 0 0 0 0
0 0 0 0 1 1 0 1 0 0 0 0 0 0 0 0
1 0 0 1 1 0 1 1 0 0 0 0 0 0 0 0

0 0 2 1 0 1 0 3 0 0 0 1 1 0 0 0
0 1 0 0 0 0 0 0 0 0 0 1 0 0 0 1
0 0 2 0 0 0 2 0 0 0 0 1 2 1 1 0
0 0 0 0 0 0 0 0 0 0 1 0 0 0 1 1

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 1 0 1 0 0 1 0 1
0 3 0 0 0 0 1 0 0 0 0 0 0 0 0 0

47 05 19 00 24 07 07 49
J = bc54c060 71f10000
key = 0001101P0110001P01?01?0P1?00100P0101001P0000??0P111?11?P?100011P
88
Unknown bit count = 8

key = 0001101001100010010011001000100101010010000011011110110001000110
key = 1a624c89 520dec46
E  = f069a40a a8a70000
E* = c0e9f28f 810f0000
C' = fa719d15

E = dabefbef 01a70000
E* = c03fa95f 29570000
C' = ea5605a1

E = ca7d095a 815f0000
E* = d57d094a 800f0000
C' = 0b00c0cb

E = 7fbf56aa 43050000
E* = 65945eba 55550000
C' = 0923cbd6

0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0
0 1 0 0 0 0 1 0 0 0 0 0 0 0 0 0
0 0 0 0 1 0 0 0 1 0 0 0 0 0 0 0

2 0 2 1 1 0 0 0 2 1 0 0 1 0 0 0
1 0 0 1 1 0 0 2 0 1 0 1 0 0 0 2
4 1 1 0 2 0 1 0 1 1 1 0 1 0 1 0
1 0 1 0 0 1 0 0 0 1 0 0 0 1 0 0

1 1 2 1 1 1 1 2 3 2 1 2 1 1 2 1
1 1 2 1 1 1 1 2 1 2 2 1 1 2 1 2
2 3 1 1 4 2 2 2 1 1 1 1 2 2 1 1
1 2 1 1 1 2 1 1 1 1 1 1 1 1 2 2

1 1 1 1 1 1 1 1 1 1 2 1 1 1 1 1
1 1 1 1 1 1 1 1 2 1 1 1 1 1 1 1
1 1 1 2 2 4 1 1 1 1 2 1 2 2 1 2
1 1 1 2 1 2 1 2 1 2 1 1 2 1 1 1

0 1 1 0 0 0 1 0 0 1 0 0 1 1 1 0
0 1 0 1 0 0 0 1 2 0 1 0 3 0 1 0
4 0 1 1 2 1 1 0 0 0 1 1 1 0 0 2
2 0 0 0 0 0 0 0 0 0 0 0 0 0 0 0

1 2 1 2 1 1 1 1 1 1 1 1 2 4 1 2
1 1 1 1 1 1 1 1 1 1 2 2 1 2 1 3
1 1 1 1 2 1 2 1 1 1 1 1 1 1 1 1
1 2 1 2 3 1 3 1 1 1 1 1 1 1 1 1

1 0 2 2 0 0 0 1 0 2 1 0 2 0 0 0
1 1 0 0 0 1 0 0 0 1 4 1 0 0 0 1
1 1 1 0 0 2 0 1 0 0 1 0 1 0 0 0
0 0 1 1 0 0 0 0 1 1 0 0 1 0 0 0

0 0 0 0 0 0 0 0 1 0 4 2 1 1 1 0
0 0 0 0 0 0 0 1 0 0 3 2 0 1 0 1
1 0 1 1 1 0 1 0 0 0 2 0 0 1 0 0
0 0 0 1 0 0 0 1 0 0 2 0 0 1 0 1

00 32 36 37 32 13 26 10
J = 02092580 d68a0000
key = 10?10?0P???1101P0011000P0000??0P00??01?P??00100P?101000P11?0101P

Unknown bit Count = 14

7774

key = 1001010011111011001100010000000100100110110010001101000011001011
key = 94fb3101 26c8d0cb
Plaintext-Ciphertext Pairs:
0 748502cd 38451097 03c70306 d8a09f10
0 38747564 38451097 78560a09 60e6d4cb
1 48691102 6acdff31 45fa285b e5adc730
1 375bd31f 6acdff31 134f7915 ac253457
2 357418da 013fec86 d8a31b2f 28bbc5cf
2 12549847 013fec86 0f317ac2 b23cb944

Possible J Values:
J # J Cnt J-Value(s):
0 1 47
1 1 5
2 1 19
3 1 0
4 1 24
5 1 7
6 1 7
7 1 49

key = 0001101001100010010011001000100101010010000011011110110001000110
tst = 0001101001100010010011001000100101010010000011011110110001000110
key = 1a624c89 520dec46
tst = 1a624c89 520dec46

J # J-Value True J-Value:
0 47 47
1 5 5
2 19 19
3 0 0
4 24 24
5 7 7
6 7 7
7 49 49

Tested Keys # = 1073
Runtime (s) = 0

Runtime per test (s) = 0.000028
Runtime max estimated (s) = 0.029784

Blog Entry (c) Tuesday, June 3, 2025, Sorting Algorithms by James Pate Williams, Jr.

First, we make sure the sorts actually work as expected and then we do some timing tests.

^^ Integer Data ^^
** Options Menu **
1 Single Sorting Tests
2 Statistical Tests
3 Create log file and log events
4 Exit
Option number: 1
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 1
Insertion Sort
PRNG Seed 1
Number of Samples 20
Maximum Sample Value 20
11 -17
-16 -16
9 -14
-17 -12
-14 -7
17 0
10 0
0 3
7 6
3 7
8 7
0 8
-7 9
11 10
-12 11
13 11
7 12
12 13
16 16
6 17
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 2
std::qsort
PRNG Seed 1
Number of Samples 20
Maximum Sample Value 20
11 -17
-16 -16
9 -14
-17 -12
-14 -7
17 0
10 0
0 3
7 6
3 7
8 7
0 8
-7 9
11 10
-12 11
13 11
7 12
12 13
16 16
6 17
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 3
Singleton's FORTRAN Sort
PRNG Seed 1
Number of Samples 20
Maximum Sample Value 20
11 -17
-16 -16
9 -14
-17 -12
-14 -7
17 0
10 0
0 3
7 6
3 7
8 7
0 8
-7 9
11 10
-12 11
13 11
7 12
12 13
16 16
6 17
-- Single Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number:
^^ Integer Data ^^
** Options Menu **
1 Single Sorting Tests
2 Statistical Tests
3 Create log file and log events
4 Exit
Option number: 2
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 1
Insertion Sort
PRNG Seed 1
Number of Samples 1000
Maximum Sample Value 1000
Number of Experiments 100
Runtimes in microseconds
Minimum runtime 524
Maximum runtime 1751
Mean runtime 802
Median runtime 668
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 2
std::qsort
PRNG Seed 1
Number of Samples 1000
Maximum Sample Value 1000
Number of Experiments 100
Runtimes in microseconds
Minimum runtime 115
Maximum runtime 1751
Mean runtime 481
Median runtime 391
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number: 3
Singleton's FORTRAN Sort
PRNG Seed 1
Number of Samples 1000
Maximum Sample Value 1000
Number of Experiments 100
Runtimes in microseconds
Minimum runtime 93
Maximum runtime 1751
Mean runtime 363
Median runtime 174
-- Statistical Sorting Tests --
1 Insertion Sort
2 std::qsort
3 Singleton's FORTRAN Sort
4 Exit Submenu
5 Exit Program
Sort option number:
; Copilot and James Pate Williams, Jr.
; 2/8/2025 - 2/9/2025
; We use the eax register for array indices
; The array base in register ecx
; The register ebx is general purpose;
;
;class SortingCPP {
;public:
;	static void InsertionSort(std::vector<T>& a)
;	{
;		for (size_t j = 1; j < a.size(); j++)
;		{
;			T key = a[j];
;			int i = j - 1;
;
;			while (i >= 0 && a[i] > key)
;			{
;				a[i + 1] = a[i];
;				i--;
;			}
;
;			a[i + 1] = key;
;		}
;	};

.MODEL FLAT, C
.STACK 4096

.DATA
    ; Allocate space for uninitialized variables
    i DWORD ?
    j DWORD ?
    key DWORD ?
    n DWORD ?
    t DwORD ?

.CODE
InsertionSortASM PROC

    ; Parameters:
    ; array = [esp + 8]
    ; n = [esp + 12]

    push ebp
    mov ebp, esp
    sub esp, 16                 ; Allocate space for local variables

    mov ecx, [ebp + 8]          ; base of array
    mov eax, [ebp + 12]         ; n number of array elements
    mov [n], eax                ; store n

    ; Initialize variables
    mov dword ptr [i], 0        ; i = 0
    mov dword ptr [j], 1        ; j = 1

for_loop:

    mov eax, [j]                ; load j into register
    mov ebx, [n]                ; load n
    cmp eax, ebx                ; compare j to n
    je  Exit                    ; we are done

    mov ebx, [ecx + eax * 4]    ; ebx = a[j]
    mov [key], ebx              ; key = a[j]
    dec eax                     ; j = j - 1 
    mov [i], eax;               ; i = j - 1
    inc eax                     ; increment
    inc eax                     ; j = j + 1
    mov [j], eax                ; store j

while_loop:

    mov eax, [i]                ; load i into register
    cmp eax, -1                 ; is i == -1 ?
    jz  end_while               ; end the while loop
    mov ebx, [ecx + eax * 4]    ; load a[i]
    mov eax, [key]              ; load key into register
    cmp ebx, eax                ; compare a[i] to key
    jle end_while               ; end the while loop
    mov eax, [i]                ; load i
    mov ebx, [ecx + eax * 4]    ; load a[i]
    inc eax                     ; eax = i + 1
    mov edx, [ecx + eax * 4]    ; load a[i + 1]
    ;mov [t], ebx               ; t = a[i]
    mov edx, ebx                ; edx = a[i]
    mov eax, [i]                ; load i again
    inc eax                     ; i + 1
    mov [ecx + eax * 4], edx    ; a[i + 1] = a[i]
    dec eax                     ; i--
    dec eax                     ; i--
    mov [i], eax                ; store updated i
    jmp while_loop              ; continue while

end_while:

    mov eax, [i]                ; load i
    inc eax                     ; eax = i + 1
    mov ebx, [key]              ; ebx = key
    mov [ecx + eax * 4], ebx    ; a[i + 1] = key

    jmp for_loop                ; continue for loop

Exit:

    mov esp, ebp
    pop ebp
    ret

InsertionSortASM ENDP
END

Blog Entry (c) Monday, June 2, 2025, Exercises from “Exterior ballistics, 1935” by James Pate Williams, Jr.

Below are some Exercises from the textbook Exterior ballistics, 1935 by former Lieutenant Commander Ernest Edward Herrman of the United States Navy Naval Academy in Annapolis, Maryland, see Chapter Four:

Density         Published Density

1.183473 1.183472
1.306584 1.306582
1.202410 1.202408
0.694163 0.694162

Exercises from Exterior ballistics, 1935
By Lieutenant Commander Ernest Edward Herrmann
Exercise 1 Page 41

Book Mine

0.157100 0.157563
0.492680 0.492810
0.517340 0.517829
0.675240 0.675580
0.844370 0.844597
1.018110 1.018197
1.014010 1.014027
1.094820 1.094910

Percentage Differences Exercise 1 Page 41

Percentage Difference [1] = 0.294283%
Percentage Difference [2] = 0.026383%
Percentage Difference [3] = 0.094477%
Percentage Difference [4] = 0.050340%
Percentage Difference [5] = 0.026880%
Percentage Difference [6] = 0.008545%
Percentage Difference [7] = 0.001676%
Percentage Difference [8] = 0.008220%

Three Altitude Related Density Calculations
Rho1 is from Exterior ballistics, 1935
Rho2 is from NASA
Rho3 is from Wikipedia

h ft Rho1 SI Rho2 SI Rho3

0 1.000000 1.226614 1.208993
1000 0.968902 1.215727 1.174010
2000 0.938772 1.204914 1.139806
3000 0.909578 1.194175 1.106369
4000 0.881292 1.183510 1.073689
5000 0.853886 1.172917 1.041752
6000 0.827332 1.162397 1.010547
7000 0.801604 1.151949 0.980062
8000 0.776676 1.141574 0.950286

Three Altitude Related Ballistic Density Calculations
Rho1 is from Exterior ballistics, 1935
Rho2 is from NASA
Rho3 is from Wikipedia
STP = Temperature = 59 F Pressure = 29.53 Humidity = 78%

Temp Press h ft Rho1 SI Rho2 SI Rho3

65 29.60 1000 0.968902 1.203134 1.161848
85 29.75 18000 0.566291 0.992154 0.656245
57 30.25 8000 0.776676 1.174418 0.977626
69 29.80 13000 0.663193 1.077573 0.801819
32 30.15 15000 0.622587 1.157366 0.822138

Exercise 3 Problems 1 - 4 m below means mine
Ra Ram are the Mayevski Retardations
Rf Rfm are the forces of Air Resistence
V Vm are the Applicable Velocities

Ra Ram Rf Rfm V Vm

567.20 567.22 229.28 229.19 2626 2626
323.85 323.86 503.50 503.29 3114 3114
91.58 91.58 2477.50 2476.50 2862 2862
62.02 62.02 4049.40 4047.83 2584 2584

Exercise 4 - Problems 1 to 5
Ra Ram BC Density

293.82 306.84 3.110358 0.989712
160.92 185.47 4.737834 1.009200
119.06 123.90 6.991926 0.952504
73.58 84.81 10.328256 0.987977
60.52 69.59 12.442559 1.080786

Exercise 5 - Problems 1 to 4
Ra Gv i-book i-mine

511.76 803.10 1.001500 0.959553
334.04 1071.14 0.602720 0.578983
90.62 922.06 0.601870 0.577421
57.30 799.71 0.616870 0.593688

Test Case from Exterior ballistics, 1935
Temperature 84 F Pressure 29.90 In Hg
Density SI = 0.959428 Book Density = 0.960
Altitude 18,000 Feet
Density SI = 0.999439 Book Density = 0.991

The differences above are probably due to the fact that Lieutenant Commander Herrmann used logarithms and logarithm tables and perhaps a slide-rule. I use double precision C/C++ real numbers.

Blog Entry © Saturday, May 31, 2025, Air Density and Ballistic Density Computations by James Pate Williams, Jr.

Online references: https://www.grc.nasa.gov/WWW/K-12/airplane/atmosmet.html and https://www.nist.gov/system/files/documents/calibrations/metv29i1p67-2.pdf

Density         Published Density

1.183473 1.183472
1.306584 1.306582
1.202410 1.202408
0.694163 0.694162

Exercises from Exterior ballistics, 1935
By Lieutenant Commander Ernest Edward Herrmann
Percentage Differences Exercise 1 Page 41

Percentage Difference [1] = 0.294283%
Percentage Difference [2] = 0.026383%
Percentage Difference [3] = 0.094477%
Percentage Difference [4] = 0.050340%
Percentage Difference [5] = 0.026880%
Percentage Difference [6] = 0.008545%
Percentage Difference [7] = 0.001676%
Percentage Difference [8] = 0.008220%

Three Altitude Related Density Calculations
Rho1 is from Exterior ballistics, 1935
Rho2 is from NASA
Rho3 is from Wikipedia

h ft Rho1 SI Rho2 SI Rho3

0 1.000000 1.226614 1.208993
1000 0.968902 1.215727 1.174010
2000 0.938772 1.204914 1.139806
3000 0.909578 1.194175 1.106369
4000 0.881292 1.183510 1.073689
5000 0.853886 1.172917 1.041752
6000 0.827332 1.162397 1.010547
7000 0.801604 1.151949 0.980062
8000 0.776676 1.141574 0.950286

Three Altitude Related Ballistic Density Calculations
Rho1 is from Exterior ballistics, 1935
Rho2 is from NASA Rho3 is from Wikipedia
STP = Temperature = 59 F Pressure = 29.53 Humidity = 78%

Temp Press h ft Rho1 SI Rho2 SI Rho3

65 29.60 1000 0.968902 1.203134 1.161848
85 29.75 18000 0.566291 0.992155 0.656245
57 30.25 8000 0.776676 1.174418 0.977626
69 29.80 13000 0.663193 1.077573 0.801819
32 30.15 15000 0.622587 1.157366 0.822138

Test Case from Exterior ballistics, 1935
Temperature 84 F Pressure 29.90 In Hg
Density SI = 0.959428 Book Density = 0.960
Altitude 18,000 Feet
Density SI = 0.999439 Book Density = 0.991

Blog Entry © Friday, May 30, 2025, Ballistic Coefficient Exercises from Exterior ballistics, 1935 Written by Ernest Edward Herrmann, C/C++ Code by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science

#include <iomanip>
#include <iostream>

double BookAnswers[] = {
    0.15710, 0.49268, 0.51734, 0.67524,
    0.84437, 1.01811, 1.01401, 1.09482 };
double MyAnswers[] = {
    0.157563, 0.492810, 0.517829, 0.675580,
    0.844597, 1.018197, 1.014027, 1.094910 };
double PD[8];

int main()
{
    std::cout << "Percentage Differences" << std::endl;
    std::cout << std::fixed << std::setprecision(6);

    for (int i = 0; i < 8; i++)
    {
        PD[i] = 100.0 * fabs(BookAnswers[i] - MyAnswers[i]) /
            (0.5 * (BookAnswers[i] + MyAnswers[i]));
        std::cout << "Percentage Difference [" << i + 1 << "] = ";
        std::cout << PD[i] << '%' << std::endl;
    }
}

Blog Entry © Wednesday, May 28, 2025, More C/C++ Siacci Method Ballistic Screenshots by James Pate Williams, Jr., BA, BS, MSWE, PhD

Blog Entry © Tuesday, May 27, 2025, Translation of a Siacci’s Method C# App to C/C++ by James Pate Williams, Jr., BA, BS, MSWE, PhD

Blog Entry © Saturday, May 25, 2025, Return to Baseball Ballistics by James Pate Williams, Jr. BA, BS, Master of Software Engineering, PhD Computer Science Return to Baseball Ballistics

Blog Entry © Friday, May 23, 2025, NASA Geopotential Altitude and the Variation of the Gravitational Acceleration by James Pate Williams, Jr., BA, BS, Master of Software Engineering, PhD Computer Science

More source code is available upon request.
Reference: 19770009539.pdf
Win32 C/C++ application output to a Win32 window:
H Z1 Z2 Z3 g
80000 79005.7 79005.7 79005.7 9.56735
80500 79493.3 79493.3 79493.3 9.56590
81000 79980.9 79980.9 79980.9 9.56446
81500 80468.3 80468.3 80468.3 9.56301
82000 80955.7 80955.7 80955.7 9.56156
82500 81443.0 81443.0 81443.0 9.56011
83000 81930.2 81930.2 81930.2 9.55867
83500 82417.4 82417.4 82417.4 9.55722
84000 82904.5 82904.5 82904.5 9.55577
84500 83391.5 83391.5 83391.5 9.55433
85000 83878.4 83878.4 83878.4 9.55288
85500 84365.3 84365.3 84365.3 9.55144
H Geopotential Altitude
Z1 Geopotential Gauss-Legendre 128 Steps
Z2 Geopotential Simpson's Rule 256 Steps
Z3 Geopotential Trapezoidal Rule 512 Steps
g gravitational acceleration in m/s/s
https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf
See Table 8. Page 9
// NASAGeopotentialAltitude.cpp : Defines the entry point for the application.
// Copyright (c) Friday, May 23, 2025 by James Pate Williams, Jr.
// BA, BS, MSWE, PhD All Applicable Rights Reserved

#include "stdafx.h"
#include "GeopotentialAltitude.h"
#include "Integration1d.h"
#include "NASAGeopotentialAltitude.h"
#include <math.h>
#include <stdio.h>
#include <tchar.h>
#include <string.h>
#include <vector>
using namespace std;

#define MAX_LOADSTRING 100

// Global Variables:
HINSTANCE hInst;								// current instance
TCHAR szTitle[MAX_LOADSTRING];					// The title bar text
TCHAR szWindowClass[MAX_LOADSTRING];			// the main window class name

// 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 APIENTRY _tWinMain(HINSTANCE hInstance,
                     HINSTANCE hPrevInstance,
                     LPTSTR    lpCmdLine,
                     int       nCmdShow)
{
	UNREFERENCED_PARAMETER(hPrevInstance);
	UNREFERENCED_PARAMETER(lpCmdLine);

 	// TODO: Place code here.
	MSG msg;
	HACCEL hAccelTable;

	// Initialize global strings
	LoadString(hInstance, IDS_APP_TITLE, szTitle, MAX_LOADSTRING);
	LoadString(hInstance, IDC_NASAGEOPOTENTIALALTITUDE, szWindowClass, MAX_LOADSTRING);
	MyRegisterClass(hInstance);

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

	hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_NASAGEOPOTENTIALALTITUDE));

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

	return (int) msg.wParam;
}

//
//  FUNCTION: MyRegisterClass()
//
//  PURPOSE: Registers the window class.
//
//  COMMENTS:
//
//    This function and its usage are only necessary if you want this code
//    to be compatible with Win32 systems prior to the 'RegisterClassEx'
//    function that was added to Windows 95. It is important to call this function
//    so that the application will get 'well formed' small icons associated
//    with it.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
	WNDCLASSEX wcex;

	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_NASAGEOPOTENTIALALTITUDE));
	wcex.hCursor		= LoadCursor(NULL, IDC_ARROW);
	wcex.hbrBackground	= (HBRUSH)(COLOR_WINDOW+1);
	wcex.lpszMenuName	= MAKEINTRESOURCE(IDC_NASAGEOPOTENTIALALTITUDE);
	wcex.lpszClassName	= szWindowClass;
	wcex.hIconSm		= LoadIcon(wcex.hInstance, MAKEINTRESOURCE(IDI_SMALL));

	return RegisterClassEx(&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)
{
   HWND hWnd;

   hInst = hInstance; // Store instance handle in our global variable

   hWnd = CreateWindow(szWindowClass, szTitle, WS_OVERLAPPEDWINDOW,
      CW_USEDEFAULT, 0, CW_USEDEFAULT, 0, NULL, NULL, hInstance, NULL);

   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_CREATE   - creates a multiple line edit control
//  WM_COMMAND	- process the application menu
//  WM_PAINT	- Paint the main window
//  WM_DESTROY	- post a quit message and return
//
//

#define ID_EDITCHILD 100
char asciiData[16384];
char asciiLine[16384];
TCHAR textData[16384];

void ConvertCharToTChar(const char* charArray, TCHAR* tcharArray, size_t tcharSize) {
#ifdef _UNICODE
	MultiByteToWideChar(CP_ACP, 0, charArray, -1, tcharArray, (int)tcharSize);
#else
    strncpy(tcharArray, charArray, tcharSize - 1);
    tcharArray[tcharSize - 1] = '\0'; // Ensure null-termination
#endif
}

LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
	int wmId, wmEvent;
	PAINTSTRUCT ps;
	HDC hdc;
	static HWND hwndEdit;
	static double zArray[] = {
		80000, 80500, 81000, 81500, 82000, 82500,
		83000, 83500, 84000, 84500, 85000, 85500 };
	static vector<double> z(12);
	static vector<double> potential1(12);
	static vector<double> potential2(12);
	static vector<double> potential3(12);
	static vector<double> potential4(12);
	size_t numberPoints = z.size();

	for (size_t i = 0; i < z.size(); i++)
	{
		z[i] = zArray[i];
	}

	switch (message)
	{
	case WM_CREATE:
	GeopotentialAltitude::GaussianLegendre(
		numberPoints, z, potential1);
	GeopotentialAltitude::SimpsonsRule(
		numberPoints, z, potential2);
	GeopotentialAltitude::TrapezoidalRule(
		numberPoints, z, potential3);

	strcpy_s(asciiData, 16384, "H\tZ1\tZ2\tZ3\tg\r\n");

	for (size_t i = 0; i < z.size(); i++)
	{
		double zi = z[i];
		double p1 = potential1[i];
		double p2 = potential2[i];
		double p3 = potential3[i];
		double g = g0 * pow(r0 / (r0 + p1), 2.0);

		sprintf_s(
			asciiLine, 16384,
			"%5.0lf\t%5.1lf\t%5.1lf\t%5.1lf\t%7.5lf\r\n",
			zi, p1, p2, p3, g);
		strcat_s(asciiData, 16384, asciiLine);
	}

	strcat_s(asciiData, "H Geopotential Altitude\r\n");
	strcat_s(asciiData, "Z1 Geopotential Gauss-Legendre 128 Steps\r\n");
	strcat_s(asciiData, "Z2 Geopotential Simpson's Rule 256 Steps\r\n");
	strcat_s(asciiData, "Z3 Geopotential Trapezoidal Rule 512 Steps\r\n");
	strcat_s(asciiData, "g gravitational acceleration in m/s/s\r\n");
	strcat_s(asciiData, "https://ntrs.nasa.gov/api/citations/19770009539/downloads/19770009539.pdf\r\n");
	strcat_s(asciiData, "See Table 8. Page 9\r\n");

	ConvertCharToTChar(asciiData, textData, sizeof(textData));

	// https://learn.microsoft.com/en-us/windows/win32/controls/use-a-multiline-edit-control

	hwndEdit = CreateWindowEx(
		0, L"EDIT",   // predefined class 
        NULL,         // no window title 
        WS_CHILD | WS_VISIBLE | WS_VSCROLL | 
        ES_LEFT | ES_MULTILINE | ES_AUTOVSCROLL, 
        0, 0, 0, 0,   // set size in WM_SIZE message 
        hWnd,         // parent window 
        (HMENU) ID_EDITCHILD,   // edit control ID 
		hInst,
        /*(HINSTANCE) GetWindowLongPtr(hWnd, GWLP_HINSTANCE),*/ 
		NULL);        // pointer not needed 
    // Add text to the window. 
	SendMessage(hwndEdit, WM_SETTEXT, 0, (LPARAM) textData); 
    return 0; 
	case WM_COMMAND:
		wmId    = LOWORD(wParam);
		wmEvent = HIWORD(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:
		hdc = BeginPaint(hWnd, &ps);
		// TODO: Add any drawing code here...
		EndPaint(hWnd, &ps);
		break;
	case WM_SIZE: 
            // Make the edit control the size of the window's client area. 
		MoveWindow(hwndEdit, 
			0, 0,                  // starting x- and y-coordinates 
            LOWORD(lParam),        // width of client area 
            HIWORD(lParam),        // height of client area 
            TRUE);                 // repaint window 
        return 0; 
	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;
}

Blog Entry © Thursday, May 22, 2025, Baseball Ballistics with Simple Drag Computation by James Pate Williams, Jr.

// https://phys.libretexts.org/Bookshelves/University_Physics/Physics_(Boundless)/3%3A_Two-Dimensional_Kinematics/3.3%3A_Projectile_Motion

#pragma once

class ClassicalBallistics
{
public:
	ClassicalBallistics(void);
	~ClassicalBallistics(void);
	// analytical formula for maximum height
	double H(double V0, double theta0, double g);
	// analytical formula for flight time
	double T(double H, double g);
	// analytical formula for velocity at apex
	double Va(double V0, double theta0);
	// analytical formula for the range
	double L(double Va, double T);
	// analytical formula for time to apex
	double ta(double T);
	// analytical formula for range to apex
	double xa(double L);
	double xa(double H, double L, double theta0);
	double Theta1(double theta0);
	double Theta1(double H, double L, double xa);
	double ClassicalBallistics::V1(double V0);
};

#include "StdAfx.h"
#include "ClassicalBallistics.h"
#include <math.h>

ClassicalBallistics::ClassicalBallistics(void)
{
}
ClassicalBallistics::~ClassicalBallistics(void)
{
}
// analytical formula for maximum height
double ClassicalBallistics::H(double V0, double theta0, double g)
{
	double V02 = V0 * V0, sinTheta0 = sin(theta0);
	double sinTheta02 = sinTheta0 * sinTheta0;

	return V02 * sinTheta02 / (g * 2.0);
}
// analytical formula for flight time
double ClassicalBallistics::T(double H, double g)
{
	return 2.0 * sqrt(2.0 * H / g);
}
// analytical formula for velocity at apex
double ClassicalBallistics::Va(double V0, double theta0)
{
	double cosTheta0 = cos(theta0);
	double V02 = V0 * V0;

	return V0 * cosTheta0;
}
// analytical formula for the range
double ClassicalBallistics::L(double Va, double T)
{
	return Va * T;
}
// analytical formula for time to apex
double ClassicalBallistics::ta(double T)
{
	return 0.5 * T;
}
// analytical formula for range to apex
double ClassicalBallistics::xa(double L)
{
	return 0.5 * L;
}
double ClassicalBallistics::xa(double H, double L, double theta0)
{
	return sqrt(L * H / tan(theta0));
}
double ClassicalBallistics::Theta1(double theta0)
{
	return -theta0;
}
double ClassicalBallistics::Theta1(double H, double L, double xa)
{
	return atan(L * H / pow(L - xa, 2.0));
}
double ClassicalBallistics::V1(double V0)
{ 
	return V0;
}

// http://www.lajpe.org/sep13/04-LAJPE-782_Chudinov.pdf

#pragma once

const double g = 9.81;
// accleration due to gravity in Metric (SI) Units
const double k = 0.000625;
// k is typical for a baseball

class DragBallistics
{

public:

	DragBallistics(void);
	~DragBallistics(void);
    double f(double theta);
    // analytical formula for maximum height
    double H(double V0, double theta0, double g, double k);
    // analytical formula for flight time
    double T(double H, double g);
    // analytical formula for velocity at apex
    double Va(double V0, double theta0, double k);
    // analytical formula for the range
	double L(double Va, double T);
    double ta(double H, double T, double Va, double k);
    double xa(double H, double L, double theta0);
    double Theta1(double H, double L, double xa);
    double V1(double V0, double theta0, double theta1, double k);
};

#include "StdAfx.h"
#include "DragBallistics.h"
#include <math.h>

DragBallistics::DragBallistics(void)
{
}
DragBallistics::~DragBallistics(void)
{
}
double DragBallistics::f(double theta)
{
	double pi = 4.0 * atan(1.0);
	double sinTheta = sin(theta), cosTheta = cos(theta),
		cosTheta2 = cosTheta * cosTheta;
    return (sinTheta / cosTheta2) + log(tan(0.5 * theta + pi / 4.0));
}
// analytical formula for maximum height
double DragBallistics::H(double V0, double theta0, double g, double k)
{
	double V02 = V0 * V0, sinTheta0 = sin(theta0),
		sinTheta02 = sinTheta0 * sinTheta0;
    return V02 * sinTheta02 / (g * (2.0 + k * V02 * sinTheta0));
}
// analytical formula for flight time
double DragBallistics::T(double H, double g)
{
	return 2.0 * sqrt(2.0 * H / g);
}
// analytical formula for velocity at apex
double DragBallistics::Va(double V0, double theta0, double k)
{
	double pi = 4.0 * atan(1.0);
	double cosTheta0 = cos(theta0), cosTheta02 = cosTheta0 * cosTheta0;
    double V02 = V0 * V0, sinTheta0 = sin(theta0);
    return V0 * cosTheta0 / sqrt(1.0 + k * V02 * (sinTheta0 + 
		cosTheta02 * log(tan(0.5 * theta0 + pi / 4.0))));
}
// analytical formula for the range
double DragBallistics::L(double Va, double T)
{
	return Va * T;
}
double DragBallistics::ta(double H, double T, double Va, double k)
{
	return 0.5 * (T - k * H * Va);
}
double DragBallistics::xa(double H, double L, double theta0)
{
	return sqrt(L * H / tan(theta0));
}
double DragBallistics::Theta1(double H, double L, double xa)
{
	return -atan(L * H / pow(L - xa, 2.0));
}
double DragBallistics::V1(double V0, double theta0, double theta1, double k)
{
	double V02 = V0 * V0;
    double sinTheta0 = sin(theta0), cosTheta0 = cos(theta0),
		cosTheta02 = cosTheta0 * cosTheta0;
    return V0 * cosTheta0 / (cos(theta1) * sqrt(1.0 + k * 
		V02 * cosTheta02 * (f(theta0) - f(theta1))));
}

// BaseballBallisticsWin32Console.cpp
// Translated from August 2017 C# application
// May 21, 2025 (c) James Pate Williams, Jr.

#include "stdafx.h"
#include "ClassicalBallistics.h"
#include "DragBallistics.h"
#include <iomanip>
#include <iostream>
#include <math.h>

void PrintResults(char title[],
				  double H, double T, double Va,
				  double L, double Ta, double xa,
				  double Theta1, double V1)
{
	std::cout << title << std::endl;
	std::cout << std::fixed << std::setprecision(2);
	std::cout << "H      = " << H << std::endl;
	std::cout << "T      = " << T << std::endl;
	std::cout << "Va     = " << Va << std::endl;
	std::cout << "L      = " << L << std::endl;
	std::cout << "Ta     = " << Ta << std::endl;
	std::cout << "xa     = " << xa << std::endl;
	std::cout << "theta1 = " << Theta1 << std::endl;
	std::cout << "V1     = " << V1 << std::endl;
}

int _tmain(int argc, _TCHAR* argv[])
{
	while (true)
	{
		char line[128] = { };
		std::cout << "V0 (m / s) or 0 to quit: ";
		std::cin.getline(line, 128);
		double V0 = atof(line);
		if (V0 == 0)
		{
			break;
		}
		std::cout << "Enter angle in degrees: ";
		std::cin.getline(line, 128);
		double theta0 = atof(line);
		double pi = 4.0 * atan(1.0);
		theta0 *= pi / 180.0;
		ClassicalBallistics cBall;
		DragBallistics dBall;
		double cH = cBall.H(V0, theta0, g);
        double dH = dBall.H(V0, theta0, g, k);
        double cT = cBall.T(cH, g);
        double dT = dBall.T(dH, g);
        double cVa = cBall.Va(V0, theta0);
        double dVa = dBall.Va(V0, theta0, k);
        double cL = cBall.L(cVa, cT);
        double dL = dBall.L(dVa, dT);
        double cta = cBall.ta(cT);
        double dta = dBall.ta(dH, dT, dVa, k);
        double cxa = cBall.xa(cL);
        double dxa = dBall.xa(dH, dL, theta0);
        double cTheta1 = 180.0 * cBall.Theta1(theta0) / pi;
        double dTheta1 = 180.0 * dBall.Theta1(dH, dL, dxa) / pi;
        double cV1 = cBall.V1(V0);
        double dV1 = dBall.V1(V0, theta0, pi * dTheta1 / 180.0, k);
		PrintResults("Classical Ballistics",
			cH, cT, cVa, cL, cta, cxa, cTheta1, cV1);
		PrintResults("Drag Ballistics",
			dH, dT, dVa, dL, dta, dxa, dTheta1, dV1);
	}
	return 0;
}