This is my project for Professor Kai Chang’s Comp 640 Advanced Computer Graphics at Auburn University in 1999. We used the Open Graphics Library also known as OpenGL. It is a hybrid C/C++ computer application.
Education
LaGrange College (1971 – 1972, 1976 – 79, BA Chemistry, GPA 3.08)
LaGrange College (1990 – 1994, BS Computer Science, GPA 4.00)
Georgia Institute of Technology (1980 – 1983, No Degree, GPA 2.79)
Auburn University (1998 – 2000, Master of Software Engineering, GPA 3.880)
Auburn University (2000 – 2005, Doctor of Philosophy, Computer Science, GPA 3.871)
Georgia Climate Data – Temperatures in Degrees Fahrenheit by James Pate Williams, Jr.
I became interested in attempting to predict average annual temperatures in the state of Georgia way back in the day. I found a neat website for temperatures for January 1895 to December 2001. Unfortunately, the website no longer exists online, however, I saved the data to an extinct PC of mine and a USB solid state drive. I used several methods to attempt predictions of the annual temperatures from 2002 to 2025. The first algorithm is polynomial least squares utilizing a 32-degree polynomial and a 75-degree polynomial. The application was written in one of my favorite programming languages, namely, C#.
Mathieu Functions by James Pate Williams, Jr.
Old adage or aphorism: “A picture is worth a thousand words.”
Parabolic Cylinder Functions by James Pate Williams, Jr.
Solution to the differential equation: d2y/dx2 = -0.25 * x * x * y – a * y which is valid for all real and complex numbers. We examine the real solutions to the second order ordinary differential equations. See “Handbook of Mathematical Functions” by Milton Abramowitz and Irene A. Stegun, Chapter 19. Parabolic Cylinder functions. Page 686 Equation 19.2.7 Recurrence formula and 19.2.5 Series solutions.
New Knight’s Tour Application by James Pate Williams, Jr.
I created a new Win32 desktop C/C++ application to solve the Knight’s Tour. The app is limited to a six-by-six chessboard with any desired starting column and row. The user can choose between two algorithms to solve an instance of the problem, namely chronological backtracking and Warnsdorff method. A complete circuit or tour is found by both algorithms in between 30 and 40 seconds for backtracking and 2 to 3 seconds for the Warnsdorff method. The starting row is 3 and starting column 4 for to find a cyclic solution.
n-Body Solar System Software by James Pate Williams, Jr.
I have solved two planetary models of our solar system. The first was by approximating the Taylor series expansion for a n-body gravitational system. The approximate equations are given below:
The second computation used a fifth order Runge-Kutta method to solve the system of 3n second order non-linear ordinary differential equations. See the code below:
using System;
using System.Collections.Generic;
using System.ComponentModel;
namespace nBodyProblem
{
public struct Vector3
{
public double x, y, z;
}
class nBodyTaylorSeries
{
private double G = 6.67384e-11; // gravitational constant
private int n; // number of masses
private List<double> Mass; // masses
private List<Vector3> X0; // initial positions
private List<Vector3> V0; // initial velocities
public nBodyTaylorSeries(
int n,
List<double> Mass,
List<Vector3> X0,
List<Vector3> V0)
{
this.n = n;
this.Mass = Mass;
this.X0 = X0;
this.V0 = V0;
}
private double Distance(Vector3 u, Vector3 v)
{
// Euclidean distance
double uvxs = Math.Pow(u.x - v.x, 2);
double uvys = Math.Pow(u.y - v.y, 2);
double uvzs = Math.Pow(u.z - v.z, 2);
return Math.Sqrt(uvxs + uvys + uvzs);
}
private Vector3 D2(
int i,
List<Vector3> X)
{
// second derivative at t0
double sumX = 0, sumY = 0, sumZ = 0;
Vector3 d2 = new Vector3();
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 3);
double m = Mass[k];
sumX += m * (X[k].x - X[i].x) / d;
sumY += m * (X[k].y - X[i].y) / d;
sumZ += m * (X[k].z - X[i].z) / d;
}
}
d2.x = G * sumX;
d2.y = G * sumY;
d2.z = G * sumZ;
return d2;
}
private Vector3 D3(
int i,
List<Vector3> X,
List<Vector3> V)
{
// third derivative at t0
double sumX1 = 0, sumY1 = 0, sumZ1 = 0;
double sumX2 = 0, sumY2 = 0, sumZ2 = 0;
Vector3 d3 = new Vector3();
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 3);
double m = Mass[k];
sumX1 += m * (V[k].x - V[i].x) / d;
sumY1 += m * (V[k].y - V[i].y) / d;
sumZ1 += m * (V[k].z - V[i].z) / d;
}
}
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 5);
double m = Mass[k];
sumX2 += m * (X[k].x - X[i].x) * (V[k].x - V[i].x) / d;
sumY2 += m * (X[k].y - X[i].y) * (V[k].y - V[i].y) / d;
sumZ2 += m * (X[k].z - X[i].z) * (V[k].z - V[i].z) / d;
}
}
d3.x = G * (sumX1 - 1.5 * sumX2);
d3.y = G * (sumY1 - 1.5 * sumY2);
d3.z = G * (sumZ1 - 1.5 * sumZ2);
return d3;
}
private Vector3 D4(
int i,
List<Vector3> X,
List<Vector3> V,
List<Vector3> A)
{
// fourth derivative at t0
double sumX1 = 0, sumY1 = 0, sumZ1 = 0;
double sumX2 = 0, sumY2 = 0, sumZ2 = 0;
double sumX3 = 0, sumY3 = 0, sumZ3 = 0;
double sumX4 = 0, sumY4 = 0, sumZ4 = 0;
Vector3 d4 = new Vector3();
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 3);
double m = Mass[k];
sumX1 += m * (A[k].x - A[i].x) / d;
sumY1 += m * (A[k].y - A[i].y) / d;
sumZ1 += m * (A[k].z - A[i].z) / d;
}
}
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 5);
double m = Mass[k];
sumX2 += m * (X[k].x - X[i].x) * (A[k].x - A[i].x) / d;
sumY2 += m * (X[k].y - X[i].y) * (A[k].y - A[i].y) / d;
sumZ2 += m * (X[k].z - X[i].z) * (A[k].z - A[i].z) / d;
}
}
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 5);
double m = Mass[k];
sumX3 += m * Math.Pow(V[k].x - V[i].x, 2) / d;
sumY3 += m * Math.Pow(V[k].y - V[i].y, 2) / d;
sumZ3 += m * Math.Pow(V[k].z - V[i].z, 2) / d;
}
}
for (int k = 0; k < n; k++)
{
if (k != i)
{
double d = Math.Pow(Distance(X[k], X[i]), 7);
double m = Mass[k];
sumX4 += m * (X[k].x - X[i].x) * Math.Pow(V[k].x - V[i].x, 2) / d;
sumY4 += m * (X[k].y - X[i].y) * Math.Pow(V[k].y - V[i].y, 2) / d;
sumZ4 += m * (X[k].z - X[i].z) * Math.Pow(V[k].z - V[i].z, 2) / d;
}
}
d4.x = G * (sumX1 - 1.5 * sumX2 - 3 * sumX3 + 15 * sumX4 / 4);
d4.y = G * (sumY1 - 1.5 * sumY2 - 3 * sumY3 + 15 * sumY4 / 4);
d4.z = G * (sumZ1 - 1.5 * sumZ2 - 3 * sumZ3 + 15 * sumZ4 / 4);
return d4;
}
public Vector3 Position(
int day0,
int day1,
int i,
BackgroundWorker worker,
DoWorkEventArgs e)
{
double t0 = day0 * 24 * 3600;
double t1 = day1 * 24 * 3600;
double delta = 60, t = t0;
List<Vector3> Q0 = new List<Vector3>();
List<Vector3> R0 = new List<Vector3>();
for (int j = 0; j < n; j++)
{
Q0.Add(X0[j]);
R0.Add(V0[j]);
}
while (t <= t1)
{
if (worker.CancellationPending)
{
e.Cancel = true;
t = t1 + delta;
}
else
{
List<Vector3> Q1 = new List<Vector3>();
List<Vector3> R1 = new List<Vector3>();
List<Vector3> d2 = new List<Vector3>();
List<Vector3> d3 = new List<Vector3>();
List<Vector3> d4 = new List<Vector3>();
for (int j = 0; j < n; j++)
{
// update derivatives
d2.Add(D2(j, Q0));
d3.Add(D3(j, Q0, R0));
Q1.Add(new Vector3());
R1.Add(new Vector3());
}
for (int j = 0; j < n; j++)
d4.Add(D4(i, Q0, R0, d2));
for (int j = 0; j < n; j++)
{
double m = Mass[j];
Vector3 Qj = new Vector3();
// update positions
Qj.x = Q0[j].x + delta * (R0[j].x + (delta * (d2[j].x / 2 + delta * (d3[j].x / 6 + delta * d4[j].x / 24)) / m));
Qj.y = Q0[j].y + delta * (R0[j].y + (delta * (d2[j].y / 2 + delta * (d3[j].y / 6 + delta * d4[j].y / 24)) / m));
Qj.z = Q0[j].z + delta * (R0[j].z + (delta * (d2[j].z / 2 + delta * (d3[j].z / 6 + delta * d4[j].z / 24)) / m));
Q1[j] = Qj;
}
for (int j = 0; j < n; j++)
{
double m = Mass[j];
Vector3 Rj = new Vector3();
// update velocities
Rj.x = R0[j].x + delta * (d2[j].x + delta * (d3[j].x / 2 + delta * d4[j].x / 6)) / m;
Rj.y = R0[j].y + delta * (d2[j].y + delta * (d3[j].y / 2 + delta * d4[j].y / 6)) / m;
Rj.z = R0[j].z + delta * (d2[j].z + delta * (d3[j].z / 2 + delta * d4[j].z / 6)) / m;
R1[j] = Rj;
}
for (int j = 0; j < n; j++)
{
Q0[j] = Q1[j];
R0[j] = R1[j];
}
t += delta;
int percentProgress = (int)(100 * t / t1);
worker.ReportProgress(percentProgress);
}
}
// update positions and velocities
X0 = new List<Vector3>();
V0 = new List<Vector3>();
for (int j = 0; j < n; j++)
{
X0.Add(Q0[j]);
V0.Add(R0[j]);
}
return Q0[i];
}
}
}
using System;
using System.Collections.Generic;
namespace nBodyProblem
{
public struct Vector3
{
public double x, y, z;
}
class RK3N
{
// function rk3n from "A Numerical Library in C for
// Scientists and Engineers" by J.T. Lau PhD p. 465
private void rk3n(ref double x, double a, double b,
double[] y, double[] ya, double[] z, double[] za,
Func<int, int, double, double[], double> fxyj,
double[] e, double[] d, bool fi, int n)
{
bool first, last=false, reject, test, ta, tb;
int j, jj;
double xl, h, hmin, ind, hl=0, absh, fhm, discry, discrz, toly, tolz, mu, mu1=0, fhy, fhz;
double[] yl = new double[n + 1];
double[] zl = new double[n + 1];
double[] k0 = new double[n + 1];
double[] k1 = new double[n + 1];
double[] k2 = new double[n + 1];
double[] k3 = new double[n + 1];
double[] k4 = new double[n + 1];
double[] k5 = new double[n + 1];
double[] ee = new double[4 * n + 1];
if (fi)
{
d[3] = a;
for (jj = 1; jj <= n; jj++)
{
d[jj + 3] = ya[jj];
d[n + jj + 3] = za[jj];
}
}
d[1] = 0.0;
xl = d[3];
for (jj = 1; jj <= n; jj++)
{
yl[jj] = d[jj + 3];
zl[jj] = d[n + jj + 3];
}
if (fi) d[2] = b - d[3];
absh = h = Math.Abs(d[2]);
if (b - xl < 0.0) h = -h;
ind = Math.Abs(b - xl);
hmin = ind * e[1] + e[2];
for (jj = 2; jj <= 2 * n; jj++)
{
hl = ind * e[2 * jj - 1] + e[2 * jj];
if (hl < hmin) hmin = hl;
}
for (jj = 1; jj <= 4 * n; jj++) ee[jj] = e[jj] / ind;
first = reject = true;
test = true;
if (fi)
{
last = true;
test = false;
}
while (true)
{
if (test)
{
absh = Math.Abs(h);
if (absh < hmin)
{
h = (h > 0.0) ? hmin : -hmin;
absh = hmin;
}
ta = (h >= b - xl);
tb = (h >= 0.0);
if ((ta && tb) || (!(ta || tb)))
{
d[2] = h;
last = true;
h = b - xl;
absh = Math.Abs(h);
}
else
last = false;
}
test = true;
if (reject)
{
x = xl;
for (jj = 1; jj <= n; jj++) y[jj] = yl[jj];
for (j = 1; j <= n; j++) k0[j] = fxyj(n, j, x, y) * h;
}
else
{
fhy = h / hl;
for (jj = 1; jj <= n; jj++) k0[jj] = k5[jj] * fhy;
}
x = xl + 0.276393202250021 * h;
for (jj = 1; jj <= n; jj++)
y[jj] = yl[jj] + (zl[jj] * 0.276393202250021 +
k0[jj] * 0.038196601125011) * h;
for (j = 1; j <= n; j++) k1[j] = fxyj(n, j, x, y) * h;
x = xl + 0.723606797749979 * h;
for (jj = 1; jj <= n; jj++)
y[jj] = yl[jj] + (zl[jj] * 0.723606797749979 +
k1[jj] * 0.261803398874989) * h;
for (j = 1; j <= n; j++) k2[j] = fxyj(n, j, x, y) * h;
x = xl + h * 0.5;
for (jj = 1; jj <= n; jj++)
y[jj] = yl[jj] + (zl[jj] * 0.5 + k0[jj] * 0.046875 + k1[jj] *
0.079824155839840 - k2[jj] * 0.001699155839840) * h;
for (j = 1; j <= n; j++) k4[j] = fxyj(n, j, x, y) * h;
x = (last ? b : xl + h);
for (jj = 1; jj <= n; jj++)
y[jj] = yl[jj] + (zl[jj] + k0[jj] * 0.309016994374947 +
k2[jj] * 0.190983005625053) * h;
for (j = 1; j <= n; j++) k3[j] = fxyj(n, j, x, y) * h;
for (jj = 1; jj <= n; jj++)
y[jj] = yl[jj] + (zl[jj] + k0[jj] * 0.083333333333333 + k1[jj] *
0.301502832395825 + k2[jj] * 0.115163834270842) * h;
for (j = 1; j <= n; j++) k5[j] = fxyj(n, j, x, y) * h;
reject =false;
fhm = 0.0;
for (jj = 1; jj <= n; jj++)
{
discry = Math.Abs((-k0[jj] * 0.5 + k1[jj] * 1.809016994374947 +
k2[jj] * 0.690983005625053 - k4[jj] * 2.0) * h);
discrz = Math.Abs((k0[jj] - k3[jj]) * 2.0 - (k1[jj] + k2[jj]) * 10.0 +
k4[jj] * 16.0 + k5[jj] * 4.0);
toly = absh * (Math.Abs(zl[jj]) * ee[2 * jj - 1] + ee[2 * jj]);
tolz = Math.Abs(k0[jj]) * ee[2 * (jj + n) - 1] + absh * ee[2 * (jj + n)];
reject = ((discry > toly) || (discrz > tolz) || reject);
fhy = discry / toly;
fhz = discrz / tolz;
if (fhz > fhy) fhy = fhz;
if (fhy > fhm) fhm = fhy;
}
mu = 1.0 / (1.0 + fhm) + 0.45;
if (reject)
{
if (absh <= hmin)
{
d[1] += 1.0;
for (jj = 1; jj <= n; jj++)
{
y[jj] = yl[jj];
z[jj] = zl[jj];
}
first = true;
if (b == x) break;
xl = x;
for (jj = 1; jj <= n; jj++)
{
yl[jj] = y[jj];
zl[jj] = z[jj];
}
}
else
h *= mu;
}
else
{
if (first)
{
first = false;
hl = h;
h *= mu;
}
else
{
fhy = mu * h / hl + mu - mu1;
hl = h;
h *= fhy;
}
mu1 = mu;
for (jj = 1; jj <= n; jj++)
z[jj] = zl[jj] + (k0[jj] + k3[jj]) * 0.083333333333333 +
(k1[jj] + k2[jj]) * 0.416666666666667;
if (b == x) break;
xl = x;
for (jj = 1; jj <= n; jj++)
{
yl[jj] = y[jj];
zl[jj] = z[jj];
}
}
}
if (!last) d[2] = h;
d[3] = x;
for (jj = 1; jj <= n; jj++)
{
d[jj + 3] = y[jj];
d[n + jj + 3] = z[jj];
}
}
// Euclidean distance
private double distance(
double x1,
double y1,
double z1,
double x2,
double y2,
double z2)
{
double xs = Math.Pow(x1 - x2, 2);
double ys = Math.Pow(y1 - y2, 2);
double zs = Math.Pow(z1 - z2, 2);
return Math.Sqrt(xs + ys + zs);
}
// n-body gravitational force
private double ftxj(int n, int j, double t, double[] x)
{
double G = 6.67384e-11; // gravitational constant
double sum = 0;
int j3 = j / 3 + 1;
int m3 = j % 3;
// remember we are now working with 3n equations
for (int k = 1; k <= n / 3; k++)
{
int k3 = 3 * (k - 1) + 1;
double d = distance(
x[k3],
x[k3 + 1],
x[k3 + 2],
x[j3],
x[j3 + 1],
x[j3 + 2]);
double denom = d * d * d;
if (denom != 0)
sum += (x[k3 + m3] - x[j3 + m3]) / denom;
}
return G * sum;
}
public Vector3 nBody(
int day0,
int day1,
int i,
int n,
List<double> Mass,
ref List<Vector3> X0,
ref List<Vector3> V0)
{
int n3 = 3 * n; // 3n equations n (x, y, z)
double a = day0 * 24 * 3600; // convert to seconds
double b = day1 * 24 * 3600; // convert to seconds
double x = a;
double[] y = new double[n3 + 1];
double[] ya = new double[n3 + 1];
double[] z = new double[n3 + 1];
double[] za = new double[n3 + 1];
double[] e = new double[4 * n3 + 1];
double[] d = new double[2 * n3 + 4];
// initialize error constraints
for (int j = 0; j <= 4 * n3; j++)
e[j] = 1.0e-12;
// create a system of 3n equations
for (int j = 0; j < n; j++)
{
int j3 = 3 * j;
ya[j3 + 1] = X0[j].x;
ya[j3 + 2] = X0[j].y;
ya[j3 + 3] = X0[j].z;
za[j3 + 1] = V0[j].x;
za[j3 + 2] = V0[j].y;
za[j3 + 3] = V0[j].z;
}
// solve the system of 3n ordinary differential equations
rk3n(ref x, a, b, y, ya, z, za, ftxj, e, d, true, n3);
// we return the position of the desired ith body
Vector3 result = new Vector3();
int i3 = 3 * i;
result.x = y[i3 + 1];
result.y = y[i3 + 2];
result.z = y[i3 + 3];
// update the positions and velocities
// wipe the slates clean
X0 = new List<Vector3>();
V0 = new List<Vector3>();
for (int j = 0; j < n; j++)
{
int j3 = 3 * j;
Vector3 X = new Vector3(); // new position
Vector3 V = new Vector3(); // new velocity
X.x = y[j3 + 1]; // 3n position components
X.y = y[j3 + 2];
X.z = y[j3 + 3];
V.x = z[j3 + 1]; // 3n velocity components
V.y = z[j3 + 2];
V.z = z[j3 + 3];
// update vectors
X0.Add(X);
V0.Add(V);
}
return result;
}
}
}
Results will be added at a later date.
Win32 C Orthogonal Polynomials
I learned about the Laguerre and Legendre polynomials when I first read “Introduction to Quantum Mechanics” by Pauling and Wilson way back in the early 1970s. I later learned about the Chebyshev and other orthogonal polynomials. Beginning on March 30, 2015, I created yet another application to graph various orthogonal polynomials in C#. A few days ago, I wrote a Win32 C application to graph the Chebyshev, Laguerre, and Legendre orthogonal polynomials. At a later date, I will probably add Hermite and Jacobi polynomials.
Win32 C Romberg Extrapolation by James Pate Williams, Jr.
// Romberg.cpp : Defines the entry point for the application.
// https://learn.microsoft.com/en-us/windows/win32/controls/create-a-simple-combo-box
#include "stdafx.h"
#include <CommCtrl.h>
#include <math.h>
#include <stdio.h>
#include "Romberg.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
// 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 RombergDialog(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_ROMBERG, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_ROMBERG));
MSG msg;
// 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.
//
ATOM MyRegisterClass(HINSTANCE hInstance)
{
WNDCLASSEXW 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_ROMBERG));
wcex.hCursor = LoadCursor(NULL, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_ROMBERG);
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, 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_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_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case ID_OPTION_ROMBERGINTEGRATION:
DialogBox(hInst, MAKEINTRESOURCE(IDD_DIALOG1), hWnd, RombergDialog);
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;
}
double f0(double x)
{
return exp(-x * x);
}
double f1(double x)
{
return x * x;
}
double f2(double x)
{
double pi = 4.0 * atan(1.0);
return sin(101.0 * pi * x);
}
double f3(double x)
{
double pi = 4.0 * atan(1.0);
return 1.0 + sin(10.0 * pi * x);
}
double f4(double x)
{
return fabs(x - 1.0 / 3.0);
}
double f5(double x)
{
return sqrt(x);
}
double f6(double x)
{
double result = 1.0;
if (x != 0.0)
result = sin(x) / x;
return result;
}
/* https://en.wikipedia.org/wiki/Romberg%27s_method#Implementation */
char table[8192];
double R1[2048], R2[2048];
void print_row(size_t i, double* R)
{
char buffer[128];
sprintf_s(buffer, 128, "R[%2zu] = ", i);
strcat_s(table, 8192, buffer);
for (size_t j = 0; j <= i; ++j)
{
sprintf_s(buffer, 128, "%.*lf ", 16, R[j]);
strcat_s(table, 8192, buffer);
}
strcat_s(table, 8192, "\r\n");
}
double RombergExtrapolation(
double(*f)(double),
double a, double b,
size_t max_steps,
double acc)
{
double *Rp = &R1[0], *Rc = &R2[0];
double h = (b - a);
Rp[0] = (f(a) + f(b))*h*.5;
print_row(0, Rp);
for (size_t i = 1; i < max_steps; ++i)
{
h /= 2.;
double c = 0;
size_t ep = 1 << (i - 1);
for (size_t j = 1; j <= ep; ++j)
{
c += f(a + (2 * j - 1)*h);
}
Rc[0] = h*c + .5*Rp[0];
for (size_t j = 1; j <= i; ++j)
{
double n_k = pow(4, j);
Rc[j] = (n_k*Rc[j - 1] - Rp[j - 1]) / (n_k - 1);
}
print_row(i, Rc);
if (i > 1 && fabs(Rp[i - 1] - Rc[i]) < acc)
{
return Rc[i];
}
double *rt = Rp;
Rp = Rc;
Rc = rt;
}
return Rp[max_steps - 1];
}
INT_PTR CALLBACK RombergDialog(HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
WORD word = LOWORD(wParam);
static TCHAR functions[7][128] =
{
TEXT("f(x) = exp(-x * x)"),
TEXT("f(x) = x * x"),
TEXT("f(x) = sin(101 * pi * x)"),
TEXT("f(x) = 1 + sin(10 * pi * x)"),
TEXT("f(x) = abs(x - 1.0 / 3.0)"),
TEXT("f(x) = sqrt(x)"),
TEXT("f(x) = sin(x) / x or 1 if x = 0")
};
static TCHAR A[1024];
static HWND cbWindow;
static int k = 0;
static int itemIndex;
switch (message)
{
case WM_INITDIALOG:
cbWindow = GetDlgItem(hDlg, IDC_COMBO1);
memset(&A, 0, sizeof(A));
for (k = 0; k < 6; k++)
{
// Add string to combobox.
SendMessage(cbWindow, CB_ADDSTRING, (WPARAM)k, (LPARAM)functions[k]);
}
// Send the CB_SETCURSEL message to display an initial item
// in the selection field
SendMessage(cbWindow, CB_SETCURSEL, (WPARAM)0, (LPARAM)0);
return (INT_PTR)TRUE;
case WM_COMMAND:
if (word == IDOK || word == IDCANCEL)
{
EndDialog(hDlg, word);
return (INT_PTR)TRUE;
}
if (HIWORD(wParam) == CBN_SELCHANGE)
// If the user makes a selection from the list:
// Send CB_GETCURSEL message to get the index of the selected list item.
// Send CB_GETLBTEXT message to get the item.
// Display the item in a messagebox.
{
itemIndex = SendMessage((HWND)lParam, (UINT)CB_GETCURSEL,
(WPARAM)0, (LPARAM)0);
TCHAR listItem[256];
(TCHAR)SendMessage((HWND)lParam, (UINT)CB_GETLBTEXT,
(WPARAM)itemIndex, (LPARAM)listItem);
}
else if (word == IDC_BUTTON1)
{
char buffer[128];
GetDlgItemTextA(hDlg, IDC_EDIT1, buffer, 128);
double a = atof(buffer);
GetDlgItemTextA(hDlg, IDC_EDIT2, buffer, 128);
double b = atof(buffer);
double integral = 0.0;
GetDlgItemTextA(hDlg, IDC_EDIT4, buffer, 128);
double acc = atof(buffer);
int max_steps = GetDlgItemInt(hDlg, IDC_EDIT5, FALSE, FALSE);
switch (itemIndex)
{
case 0:
integral = RombergExtrapolation(f0, a, b, max_steps, acc);
break;
case 1:
integral = RombergExtrapolation(f1, a, b, max_steps, acc);
break;
case 2:
integral = RombergExtrapolation(f2, a, b, max_steps, acc);
break;
case 3:
integral = RombergExtrapolation(f3, a, b, max_steps, acc);
break;
case 4:
integral = RombergExtrapolation(f4, a, b, max_steps, acc);
break;
case 5:
integral = RombergExtrapolation(f5, a, b, max_steps, acc);
break;
case 6:
integral = RombergExtrapolation(f6, a, b, max_steps, acc);
break;
}
sprintf_s(buffer, 128, "Integral = %.*lf\r\n", 16, integral);
strcat_s(table, 8192, buffer);
SetDlgItemTextA(hDlg, IDC_EDIT7, table);
}
else if (word == IDC_BUTTON3)
{
table[0] = '\0';
SetDlgItemTextA(hDlg, IDC_EDIT7, table);
}
}
return (INT_PTR)FALSE;
}
Numerical Explorations 