/*
Author: Pate Willimas c 1995
The following program is a translation of the FORTRAN77
subroutine found in "Elementary Numerical Analysis
An Algorithmic Approach Third Edition" (c) 1980 by S. D.
Conte and Carl de Boor pages 343-344. The program
uses Romberg extrapolation to find the integral of a
function. Also see: https://dlmf.nist.gov/3.5#E10
*/
#include <math.h>
#include <stdio.h>
#include <vector>
#include "Bessel.h"
typedef double real;
long fe;
real t[10][10];
real f(real x)
{
fe++;
return(expl(-x * x));
}
real g(real x)
{
fe++;
std::vector<real> j(4);
Bessel::bessj(x, 0, j);
return exp(-x) * j[0];
}
real romberg(
real a, real b,
int start, int row,
real (*fx)(real))
{
int i, /*j, */k, m;
real h, ratio, sum;
m = start;
h = (b - a) / m;
sum = 0.5 * (fx(a) + fx(b));
if (m > 1)
for (i = 1; i < m; i++)
sum += fx(a + i * h);
t[1][1] = h * sum;
//printf_s("Romberg T-Table\n%9.7lf\n", t[1][1]);
if (row < 2) return(t[1][1]);
for (k = 2; k <= row; k++)
{
h = 0.5 * h;
m *= 2;
sum = 0.0;
for (i = 1; i <= m; i += 2)
sum += fx(a + h * i);
t[k][1] = 0.5 * t[k - 1][1] + sum * h;
for (i = 1; i < k; i++)
{
t[k - 1][i] = t[k][i] - t[k - 1][i];
t[k][i + 1] = t[k][i] - t[k - 1][i] /
(powl(4.0, (real)i) - 1.0);
//for (j = 1; j < k; j++)
//printf_s("%9.7lf ", t[k][j]);
}
//printf_s("\n");
}
if (row < 3) return(t[2][2]);
//printf_s("Table of ratios\n");
for (k = 1; k <= row - 2; k++) {
for (i = 1; i <= k; i++)
{
if (t[k + 1][i] == 0.0)
ratio = 0.0;
else
ratio = t[k][i] / t[k + 1][i];
t[k][i] = ratio;
//printf_s("%5.2lf ", ratio);
}
//printf_s("\n");
}
return(t[row][row - 1]);
}
real CompositeTrapezoidalRule(
real a, real b, int n,
real(*fx)(real)) {
real endPts = 0.5 * (fx(a) + fx(b));
real sum = 0, xk = 0.0;
real h = (b - a) / n;
for (int k = 1; k <= n - 1; k++)
{
xk = a + k * h;
sum += fx(xk);
}
return h * (endPts + sum);
}
double CompositeSimpsonsRule(
double a, double b, int n,
real(*fx)(real))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
void Integrate(int trial, real(*fx)(real))
{
real a = 0.0, b = 0.0;
int row, start;
if (trial == 0)
{
a = 0.0;
b = 1.0;
}
else
{
a = 0.0;
b = 30.0;
}
if (trial == 0)
printf("Romberg integration of f(x) = exp(- x * x)\n");
else
printf("Romberg integration of g(x) = exp(-x) * J0(x)\n");
printf("number of trapezoidal intervals = ");
scanf_s("%d", &start);
printf("number of rows in table (<= 8) = ");
scanf_s("%d", &row);
fe = 0;
printf_s("integral = %14.11lf\tevals = %d\n",
romberg(a, b, start, row, fx), fe);
printf("number of trapezoidal intervals = ");
scanf_s("%d", &row);
fe = 0;
printf_s("integral = %14.11lf\tevals = %d\n",
CompositeTrapezoidalRule(a, b, row, fx), fe);
fe = 0;
printf("number of Simpson's intervals = ");
scanf_s("%d", &row);
printf_s("integral = %14.11lf\tevals = %d\n",
CompositeSimpsonsRule(a, b, row, fx), fe);
if (trial == 0)
printf_s("integral = %14.11lf\n", 0.74682413279);
else
{
real pi = 4.0 * atan(1.0);
real si = sin(pi * 45.0 / 180.0);
printf_s("integral = %14.11lf\n", si);
}
}
int main(void)
{
Integrate(0, f);
Integrate(1, g);
return(0);
}
Blog Entry © Sunday, March 29, 2026, James Pate Williams, Jr. Properties of Determinants
// Reference: "A Course in Computational
// Algebraic Number Theory by Henri Cohen
// Chapter 2 Algorithm 2.2.1 and Algorithm
// 2.2.3
#pragma once
#include <vector>
class Determinants
{
public:
static bool GaussianElimnation(
int n,
std::vector<std::vector<double>>& M,
std::vector<double> B,
std::vector<double> X);
static bool Determinant(
double& det, int n,
std::vector<std::vector<double>>& M);
};
#include "Determinants.h"
bool Determinants::GaussianElimnation(
int n,
std::vector<std::vector<double>>& M,
std::vector<double> B,
std::vector<double> X) {
double d = 0.0;
int i = 0, j = -1, k = 0, l = 0;
std::vector<double> C(n);
Step2:
j++;
if (j == n)
goto Step6;
for (i = j; i < n; i++) {
if (M[i][j] > 0)
break;
}
if (i == n)
return false;
if (i > j) {
for (l = j; j < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
temp = B[i];
B[i] = B[j];
B[j] = temp;
}
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
for (k = j + 1; k < n; k++) {
B[k] -= C[k] * B[j];
}
goto Step2;
Step6:
for (i = n - 1; i >= 0; i--) {
double sum = 0.0;
for (j = i + 1; j < n; j++) {
sum += M[i][j] * X[j];
}
X[i] = (B[i] - sum) / M[i][i];
}
return true;
}
bool Determinants::Determinant(
double& det, int n,
std::vector<std::vector<double>>& M) {
double d = 0.0, x = 1.0;
int i = 0, j = -1, k = 0, l = 0;
std::vector<double> C(n, 0);
Step2:
j++;
if (j == n) {
det = x;
return true;
}
for (i = j; i < n; i++) {
if (M[i][j] != 0.0)
break;
}
if (i == n) {
det = 0.0;
return false;
}
if (i > j) {
for (l = j; j < n; l++) {
double temp = M[i][l];
M[i][l] = M[j][l];
M[j][l] = temp;
}
x = -x;
}
d = 1.0 / M[j][j];
for (k = j + 1; k < n; k++) {
C[k] = d * M[k][j];
}
for (k = j + 1; k < n; k++) {
for (l = j + 1; l < n; l++) {
M[k][l] -= C[k] * M[j][l];
}
}
x *= M[j][j];
goto Step2;
}
// ComputeRealDeterminant.cpp : Defines the entry point for the application.
//
#include "framework.h"
#include "ComputeRealDeterminant.h"
#include "Determinants.h"
#include <string>
#include <vector>
#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], numWstr[128];
WCHAR inpStr[16384], outStr[16384];
// 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 DataInputDialog(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_COMPUTEREALDETERMINANT, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_COMPUTEREALDETERMINANT));
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_COMPUTEREALDETERMINANT));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_COMPUTEREALDETERMINANT);
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_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_START:
DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG), hWnd,
DataInputDialog);
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 DataInputDialog(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam) {
static int n = 0;
static std::wstring wstr;
static std::vector < std::vector<double>> A;
static HFONT hFont = nullptr;
static HWND hEditMultiline1 = nullptr;
static HWND hEditMultiline2 = nullptr;
UNREFERENCED_PARAMETER(lParam);
switch (message) {
case WM_INITDIALOG:
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_N, CB_SETCURSEL, 0, 0);
hFont = CreateFont(
-MulDiv(10, 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")
);
hEditMultiline1 = 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_AUTOVSCROLL,
10, 100, 250, 250, // Position and size
hDlg, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE1, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
hEditMultiline2 = CreateWindowEx(
WS_EX_CLIENTEDGE, // Extended style for sunken border
TEXT("EDIT"), // Class name
TEXT(""), // Initial text (can be blank)
WS_CHILD | WS_VISIBLE | ES_LEFT | ES_READONLY,
120, 360, 120, 30, // Position and size
hDlg, // Parent window handle
(HMENU)IDC_EDIT_MULTILINE2, // Unique control ID
hInst, // Application instance
NULL // Extra parameter
);
break;
case WM_COMMAND:
if (LOWORD(wParam) == IDC_BUTTON_COMPUTE) {
GetDlgItemText(hDlg, IDC_COMBO_N, inpStr, 16384);
WCHAR* endptr = nullptr;
n = (int)wcstol(inpStr, &endptr, 10);
A.resize(n);
GetDlgItemText(hDlg, IDC_EDIT_MULTILINE1, inpStr, 16384);
int k = 0;
for (int i = 0; i < n; i++) {
for (int j = 0; j < n; j++) {
A[i].resize(n);
int l = 0;
while (inpStr[k] != L' ') {
numWstr[l++] = inpStr[k++];
}
k++;
numWstr[l++] = L'\0';
A[i][j] = (double)wcstof(numWstr, &endptr);
}
}
double det = 0.0;
bool result = Determinants::Determinant(det, n, A);
swprintf_s(line, 128, L"%lf", det);
SetDlgItemText(hDlg, IDC_EDIT_MULTILINE2, line);
break;
}
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL) {
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
}
return (INT_PTR)FALSE;
}
//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by ComputeRealDeterminant.rc
#define IDS_APP_TITLE 103
#define IDR_MAINFRAME 128
#define IDD_COMPUTEREALDETERMINANT_DIALOG 102
#define IDD_ABOUTBOX 103
#define IDM_ABOUT 104
#define IDM_EXIT 105
#define IDI_COMPUTEREALDETERMINANT 107
#define IDI_SMALL 108
#define IDC_COMPUTEREALDETERMINANT 109
#define IDC_MYICON 2
#ifndef IDC_STATIC
#define IDC_STATIC -1
#endif
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS
#define _APS_NO_MFC 130
#define _APS_NEXT_RESOURCE_VALUE 129
#define _APS_NEXT_COMMAND_VALUE 32771
#define _APS_NEXT_CONTROL_VALUE 1000
#define _APS_NEXT_SYMED_VALUE 110
#endif
#endif
#define IDD_DATA_INPUT_DIALOG 1000
#define IDC_COMBO_N 1010
#define IDC_EDIT_MULTILINE1 1020
#define IDC_EDIT_MULTILINE2 1030
#define IDC_BUTTON_COMPUTE 1040
#define IDM_START 32771
//Microsoft Visual C++ generated resource script.
//
#include "resource.h"
#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
#ifndef APSTUDIO_INVOKED
#include "targetver.h"
#endif
#define APSTUDIO_HIDDEN_SYMBOLS
#include "windows.h"
#undef APSTUDIO_HIDDEN_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS
#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU)
LANGUAGE 9, 1
/////////////////////////////////////////////////////////////////////////////
//
// Icon
//
// Icon with lowest ID value placed first to ensure application icon
// remains consistent on all systems.
IDI_COMPUTEREALDETERMINANT ICON "ComputeRealDeterminant.ico"
IDI_SMALL ICON "small.ico"
/////////////////////////////////////////////////////////////////////////////
//
// Menu
//
IDC_COMPUTEREALDETERMINANT MENU
BEGIN
POPUP "&Begin"
BEGIN
MENUITEM "&Start", IDM_START
MENUITEM SEPARATOR
MENUITEM "E&xit", IDM_EXIT
END
POPUP "&Help"
BEGIN
MENUITEM "&About ...", IDM_ABOUT
END
END
/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//
IDC_COMPUTEREALDETERMINANT ACCELERATORS
BEGIN
"?", IDM_ABOUT, ASCII, ALT
"/", IDM_ABOUT, ASCII, ALT
END
/////////////////////////////////////////////////////////////////////////////
//
// Dialog
//
IDD_ABOUTBOX DIALOGEX 0, 0, 170, 62
STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "About ComputeRealDeterminant"
FONT 8, "MS Shell Dlg"
BEGIN
ICON IDI_COMPUTEREALDETERMINANT,IDC_STATIC,14,14,21,20
LTEXT "ComputeRealDeterminant, Version 1.0",IDC_STATIC,42,14,114,8,SS_NOPREFIX
LTEXT "Copyright (c) 2026",IDC_STATIC,42,26,114,8
DEFPUSHBUTTON "OK",IDOK,113,41,50,14,WS_GROUP
END
IDD_DATA_INPUT_DIALOG DIALOGEX 0, 0, 280, 260
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog"
FONT 10, "Courier New", 700
BEGIN
LTEXT "n:", IDC_STATIC, 10, 10, 30, 10
COMBOBOX IDC_COMBO_N, 70, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "input matrix:", IDC_STATIC, 10, 30, 200, 10
LTEXT "determinant:", IDC_STATIC, 10, 185, 200, 10
PUSHBUTTON "Compute", IDC_BUTTON_COMPUTE, 10, 220, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 220, 220, 50, 16
END
/////////////////////////////////////////////////////////////////////////////
//
// DESIGNINFO
//
#ifdef APSTUDIO_INVOKED
GUIDELINES DESIGNINFO
BEGIN
IDD_ABOUTBOX, DIALOG
BEGIN
LEFTMARGIN, 7
RIGHTMARGIN, 163
TOPMARGIN, 7
BOTTOMMARGIN, 55
END
END
#endif // APSTUDIO_INVOKED
#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
"resource.h\0"
END
2 TEXTINCLUDE
BEGIN
"#ifndef APSTUDIO_INVOKED\r\n"
"#include ""targetver.h""\r\n"
"#endif\r\n"
"#define APSTUDIO_HIDDEN_SYMBOLS\r\n"
"#include ""windows.h""\r\n"
"#undef APSTUDIO_HIDDEN_SYMBOLS\r\n"
"\0"
END
3 TEXTINCLUDE
BEGIN
"\r\n"
"\0"
END
#endif // APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// String Table
//
STRINGTABLE
BEGIN
IDC_COMPUTEREALDETERMINANT "COMPUTEREALDETERMINANT"
IDS_APP_TITLE "ComputeRealDeterminant"
END
#endif
/////////////////////////////////////////////////////////////////////////////
#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
/////////////////////////////////////////////////////////////////////////////
#endif // not APSTUDIO_INVOKED
Blog Entry © Friday, March 27, 2026, by James Pate Williams, Jr. Orthogonal Polynomial Application
Blog Entry © Tuesday, March 24, 2026, by James Pate Williams, Jr. Cubature on a Disk
// CubatureDisk.cpp (c) Tuesday, March 24, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include <iomanip>
#include <iostream>
#include <vector>
#include "Mueller.h"
static void SelectionSort(
std::vector<double>& x,
std::vector<double>& w,
int n) {
for (int i = 0; i < n - 1; i++) {
for (int j = i + 1; j < n; j++) {
if (x[i] > x[j]) {
double t = x[i];
x[i] = x[j];
x[j] = t;
t = w[i];
w[i] = w[j];
w[j] = t;
}
}
}
}
static double fr(double rho, int i, int j) {
return fabs(rho) * pow(rho, i + j);
}
static double ft(double t, int i, int j) {
double t1 = 1.0 - t * t;
return (1.0 / sqrt(t1)) * pow(t1, i / 2.0) * pow(t, j);
}
static double GaussLegendre(
std::vector<double> x, std::vector<double> w,
double (*fx)(double, int, int), int i, int j, int n) {
double integral = 0.0;
for (int k = 0; k < n; k++) {
integral += fx(x[k], i, j) * w[k];
}
return integral;
}
static double Integrate(
double (*fx)(double, int, int),
int i, int j, int n,
unsigned int seed) {
std::string msg;
std::vector<double> x(n);
std::vector<double> weights(n);
Mueller::OrthogonalPolynomialRoots(
"Legendre", "Zeros", n, n, seed, msg, x);
for (int i = 0; i < n; i++) {
double xi = x[i];
double dpn = Mueller::g("Legendre", "Extrema",
xi, n, 0, x);
weights[i] = 2.0 / ((1.0 - xi * xi) * dpn * dpn);
}
SelectionSort(x, weights, n);
return GaussLegendre(x, weights, fx, i, j, n);
}
static double SimpsonsRule(
int i, int j, int n, double a, double b,
double (*fx)(double, int, int))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int k = 1; k < n; k += 2)
{
s += fx(x, i, j);
x += h2;
}
x = a + h2;
for (int k = 2; k < n; k += 2)
{
t += fx(x, i, j);
x += h2;
}
return h * (fx(a, i, j) + 4 * s + 2 * t + fx(b, i, j)) / 3.0;
}
double MonteCarlo(int i, int j, int N) {
double integral = 0.0;
for (int k = 0; k < N; k++)
{
double x = (double)rand() / RAND_MAX;
double y = (double)rand() / RAND_MAX;
if (x * x + y * y <= 1)
integral += pow(x, i) * pow(y, j);
}
return 4.0 * integral / N;
}
static double Integrand(std::vector<double> x, double n) {
double power = 1.0;
for (int k = 0; k < n; k++) {
power *= sqrt(fabs(x[k] - 1.0 / 3.0));
}
return power;
}
static double MonteCarloPower(int N, int n) {
double integration = 0.0;
std::vector<double> x(n);
for (int i = 0; i < N; i++) {
for (int k = 0; k < n; k++) {
x[k] = (double)rand() / RAND_MAX;
}
integration += Integrand(x, n);
}
return integration / N;
}
static double MonteCarloPi(int N) {
double integral = 0.0;
for (int i = 0; i < N; i++) {
double x = (double)rand() / RAND_MAX;
double y = (double)rand() / RAND_MAX;
if (sqrt(x * x + y * y) <= 1)
integral++;
}
return 4.0 * integral / N;
}
int main() {
int i = 4, j = 4, n = 0, N = 10000000;
unsigned int seed = 1;
std::cout << "Integral # 1" << std::endl;
for (int n = 5; n <= 25; n += 5)
{
double integral1 = Integrate(fr, i, j, n, seed++);
double integral2 = Integrate(ft, i, j, n, seed++);
double integral = integral1 * integral2;
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << n << "\t " << integral << std::endl;
}
for (int n = 50000000; n <= 90000000; n += 10000000)
{
double integral1 = SimpsonsRule(i, j, n, -1 + 1.0e-15, 1 - 1.0e-15, fr);
double integral2 = SimpsonsRule(i, j, n, -1 + 1.0e-15, 1 - 1.0e-15, ft);
double integral = integral1 * integral2;
std::cout << "Simpson ";
if (n < 10000000)
std::cout << n << "\t\t";
else
std::cout << n << '\t';
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral << std::endl;
}
double integral3 = MonteCarlo(i, j, N);
std::cout << "Monte Carlo ";
std::cout << N << "\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral3 << std::endl;
N = 10000000;
n = 6;
double integral4 = MonteCarloPower(N, n);
double actual4 = pow(2.0 / 3.0, n) * pow(pow(2.0 / 3.0, 1.5) +
pow(1.0 / 3.0, 1.5), n);
std::cout << "Integral # 2" << std::endl;
std::cout << "Monte Carlo ";
std::cout << N << "\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral4 << std::endl;
std::cout << "actual\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << actual4 << std::endl;
double perror = 100.0 * fabs(integral4 - actual4) /
fabs(actual4);
std::cout << "% error\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << perror << std::endl;
double actual5 = 4.0 * atan(1.0);
double integral5 = MonteCarloPi(N);
std::cout << "Integral # 3" << std::endl;
std::cout << "Monte Carlo ";
std::cout << N << "\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral5 << std::endl;
std::cout << "actual\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << actual5 << std::endl;
perror = 100.0 * fabs(integral5 - actual5) /
fabs(actual5);
std::cout << "% error\t\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << perror << std::endl;
return 0;
}
#pragma once
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Numerical Analysis: An Algorithmic
* Approach" by S. D. Conte and Carl de Boor
* Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
* Translated to better C++ on Saturday, March 21,
* 2026
*/
#include <string>
#include <vector>
class Mueller
{
private:
static double epsilon1, epsilon2;
public:
static double f(std::string name,
double x, int n, int m,
std::vector<double>& zeros);
static double g(
std::string name,
std::string graph,
double x, int n, int m,
std::vector<double>& zeros);
static double OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros);
static void OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros);
};
#include "Mueller.h"
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
double Mueller::epsilon1 = 1.0e-12;
double Mueller::epsilon2 = 1.0e-15;
double Mueller::f(std::string name,
double x, int n, int m,
std::vector<double>& zeros)
{
double a1n = 1, a2n = 1, a3n = 1, a4n = 1;
double f0 = 1, f1 = x, f2 = 1;
if (name == "Laguerre")
f1 = -x + 1;
for (int i = 2; i <= n; i++)
{
if (name == "Chebyshev")
{
a1n = 1.0;
a2n = 0.0;
a3n = 2.0;
a4n = 1.0;
}
else if (name == "Laguerre")
{
a1n = i - 1 + 1;
a2n = 2 * (i - 1) + 1;
a3n = -1.0;
a4n = i - 1;
}
else if (name == "Legendre")
{
a1n = i + 1 - 1;
a2n = 0.0;
a3n = 2.0 * (i - 1) + 1;
a4n = i - 1;
}
f2 = ((a2n + a3n * x) * f1 - a4n * f0) / a1n;
f0 = f1;
f1 = f2;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) < epsilon2)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
}
return f2;
}
double Mueller::g(std::string name,
std::string graph, double x,
int n, int m, std::vector<double>& zeros)
{
double f0 = 1, f1 = x, f2 = 1;
if (graph != "Extrema")
return f(name, x, n, m, zeros);
else
{
if (name == "Laguerre") {
if (x == 1.0)
x -= 0.0000001;
else if (x == -1.0)
x += 0.0000001;
}
double g0x = 1.0, g1x = 1.0, g2x = 1.0;
if (name == "Chebyshev") {
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
else if (name == "Laguerre") {
g0x = -n;
g1x = n;
g2x = x;
}
else if (name == "Legendre")
{
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
if (fabs(g2x) < epsilon2)
{
std::cout << "g2x is too small!";
return 0.0;
}
f0 = f(name, x, n - 1, 0, zeros);
f1 = f(name, x, n, 0, zeros);
f2 = (g1x * f1 + g0x * f0) / g2x;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) == 0.0)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
return f2;
}
}
double Mueller::OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros)
{
double x3 = 0.0;
int i = 2;
while (true)
{
double h1 = x1 - x0;
if (fabs(h1) < epsilon2)
{
msg = "h1 is too small!";
return 0.0;
}
double h2 = x2 - x1;
if (fabs(h2) < epsilon2)
{
msg = "h2 is too small!";
return 0.0;
}
double f0 = g(name, graph, x0, n, m, zeros);
double f1 = g(name, graph, x1, n, m, zeros);
double f2 = g(name, graph, x2, n, m, zeros);
double g1 = (f1 - f0) / h1;
double g2 = (f2 - f1) / h2;
if (fabs(f2) < epsilon2)
{
x3 = x2;
break;
}
if (fabs(h1 + h2) < epsilon2)
{
msg = "h1 + h2 is too small!";
return 0.0;
}
double k2 = (g2 - g1) / (h1 + h2);
double c1 = g2 + h2 * k2;
double d = c1 * c1 - 4.0 * f2 * k2;
if (d < 0.0)
d = 0.0;
double s1 = sqrt(d);
double m1 = c1 - s1;
double p1 = c1 + s1;
double denom;
if (fabs(m1) > fabs(p1))
denom = m1;
else
denom = p1;
if (fabs(denom) < epsilon2)
{
msg = "denom is too small!";
return 0.0;
}
double h3 = -2.0 * f2 / denom;
x3 = x2 + h3;
double f3 = g(name, graph, x3, n, m, zeros);
if (fabs(f3) < epsilon2)
break;
if (fabs(h3) < epsilon2)
{
msg = "h3 is too small!";
return 0.0;
}
double g3 = (f3 - f2) / h3;
double delta = fabs(x3 - x2);
if (delta < epsilon1)
break;
i++;
if (i > maxIt)
break;
x0 = x1;
x1 = x2;
x2 = x3;
if (fabs(x1 - x0) < epsilon1)
break;
if (fabs(x2 - x1) < epsilon1)
break;
}
return x3;
}
void Mueller::OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros) {
double epsilon1 = 1.0e-10;
epsilon2 = 1.0e-20;
double x0 = -0.5, x1 = 0.0, x2 = 0.5;
double z;
int maxIt = 128;
zeros.resize(nr);
srand(seed);
for (int m = 0; m < nr; m++) {
x0 = (double)rand() / RAND_MAX;
x1 = (double)rand() / RAND_MAX;
x2 = (double)rand() / RAND_MAX;
z = zeros[m] = OrthogonalPolynomialRoot(
name, graph, n, epsilon1, epsilon2,
maxIt, m, x0, x1, x2, msg, zeros);
}
}
Blog Entry © Monday, March 23, 2026, by James Pate Williams, Jr. Gauss-Legendre Quadrature
#pragma once
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Numerical Analysis: An Algorithmic
* Approach" by S. D. Conte and Carl de Boor
* Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
* Translated to better C++ on Saturday, March 21,
* 2026
*/
#include <string>
#include <vector>
class Mueller
{
private:
static double epsilon1, epsilon2;
public:
static double f(std::string name,
double x, int n, int m,
std::vector<double>& zeros);
static double g(
std::string name,
std::string graph,
double x, int n, int m,
std::vector<double>& zeros);
static double OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros);
static void OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros);
};
#include "Mueller.h"
#include <cmath>
#include <iostream>
#include <string>
#include <vector>
double Mueller::epsilon1 = 1.0e-12;
double Mueller::epsilon2 = 1.0e-15;
double Mueller::f(std::string name,
double x, int n, int m,
std::vector<double>& zeros)
{
double a1n = 1, a2n = 1, a3n = 1, a4n = 1;
double f0 = 1, f1 = x, f2 = 1;
if (name == "Laguerre")
f1 = -x + 1;
for (int i = 2; i <= n; i++)
{
if (name == "Chebyshev")
{
a1n = 1.0;
a2n = 0.0;
a3n = 2.0;
a4n = 1.0;
}
else if (name == "Laguerre")
{
a1n = i - 1 + 1;
a2n = 2 * (i - 1) + 1;
a3n = -1.0;
a4n = i - 1;
}
else if (name == "Legendre")
{
a1n = i + 1 - 1;
a2n = 0.0;
a3n = 2.0 * (i - 1) + 1;
a4n = i - 1;
}
f2 = ((a2n + a3n * x) * f1 - a4n * f0) / a1n;
f0 = f1;
f1 = f2;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) < epsilon2)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
}
return f2;
}
double Mueller::g(std::string name,
std::string graph, double x,
int n, int m, std::vector<double>& zeros)
{
double f0 = 1, f1 = x, f2 = 1;
if (graph != "Extrema")
return f(name, x, n, m, zeros);
else
{
if (name == "Laguerre") {
if (x == 1.0)
x -= 0.0000001;
else if (x == -1.0)
x += 0.0000001;
}
double g0x = 1.0, g1x = 1.0, g2x = 1.0;
if (name == "Chebyshev") {
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
else if (name == "Laguerre") {
g0x = -n;
g1x = n;
g2x = x;
}
else if (name == "Legendre")
{
g0x = n;
g1x = -n * x;
g2x = 1.0 - x * x;
}
if (fabs(g2x) < epsilon2)
{
std::cout << "g2x is too small!";
return 0.0;
}
f0 = f(name, x, n - 1, 0, zeros);
f1 = f(name, x, n, 0, zeros);
f2 = (g1x * f1 + g0x * f0) / g2x;
// deflate the function
for (int i = 0; i < m; i++)
{
double z = zeros[i];
double denom = x - z;
if (fabs(denom) == 0.0)
zeros[i] = 1.001 * z;
else
f2 /= denom;
}
return f2;
}
}
double Mueller::OrthogonalPolynomialRoot(
std::string name,
std::string graph, int n,
double epsilon1, double epsilon2,
int maxIt, int m, double x0,
double x1, double x2, std::string& msg,
std::vector<double>& zeros)
{
double x3 = 0.0;
int i = 2;
while (true)
{
double h1 = x1 - x0;
if (fabs(h1) < epsilon2)
{
msg = "h1 is too small!";
return 0.0;
}
double h2 = x2 - x1;
if (fabs(h2) < epsilon2)
{
msg = "h2 is too small!";
return 0.0;
}
double f0 = g(name, graph, x0, n, m, zeros);
double f1 = g(name, graph, x1, n, m, zeros);
double f2 = g(name, graph, x2, n, m, zeros);
double g1 = (f1 - f0) / h1;
double g2 = (f2 - f1) / h2;
if (fabs(f2) < epsilon2)
{
x3 = x2;
break;
}
if (fabs(h1 + h2) < epsilon2)
{
msg = "h1 + h2 is too small!";
return 0.0;
}
double k2 = (g2 - g1) / (h1 + h2);
double c1 = g2 + h2 * k2;
double d = c1 * c1 - 4.0 * f2 * k2;
if (d < 0.0)
d = 0.0;
double s1 = sqrt(d);
double m1 = c1 - s1;
double p1 = c1 + s1;
double denom;
if (fabs(m1) > fabs(p1))
denom = m1;
else
denom = p1;
if (fabs(denom) < epsilon2)
{
msg = "denom is too small!";
return 0.0;
}
double h3 = -2.0 * f2 / denom;
x3 = x2 + h3;
double f3 = g(name, graph, x3, n, m, zeros);
if (fabs(f3) < epsilon2)
break;
if (fabs(h3) < epsilon2)
{
msg = "h3 is too small!";
return 0.0;
}
double g3 = (f3 - f2) / h3;
double delta = fabs(x3 - x2);
if (delta < epsilon1)
break;
i++;
if (i > maxIt)
break;
x0 = x1;
x1 = x2;
x2 = x3;
if (fabs(x1 - x0) < epsilon1)
break;
if (fabs(x2 - x1) < epsilon1)
break;
}
return x3;
}
void Mueller::OrthogonalPolynomialRoots(
std::string name,
std::string graph, int n, int nr,
unsigned int seed, std::string& msg,
std::vector<double>& zeros) {
double epsilon1 = 1.0e-10;
epsilon2 = 1.0e-20;
double x0 = -0.5, x1 = 0.0, x2 = 0.5;
double z;
int maxIt = 128;
zeros.resize(nr);
srand(seed);
for (int m = 0; m < nr; m++) {
x0 = (double)rand() / RAND_MAX;
x1 = (double)rand() / RAND_MAX;
x2 = (double)rand() / RAND_MAX;
z = zeros[m] = OrthogonalPolynomialRoot(
name, graph, n, epsilon1, epsilon2,
maxIt, m, x0, x1, x2, msg, zeros);
}
}
// Quadrature1d.cpp (c) Saturday, March 21, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: https://dlmf.nist.gov/3.5#x
#include <iomanip>
#include <iostream>
#include <vector>
#include "Mueller.h"
static void SelectionSort(
std::vector<double>& x,
std::vector<double>& w,
int n) {
for (int i = 0; i < n - 1; i++) {
for (int j = i + 1; j < n; j++) {
if (x[i] > x[j]) {
double t = x[i];
x[i] = x[j];
x[j] = t;
t = w[i];
w[i] = w[j];
w[j] = t;
}
}
}
}
static double fx(double x) {
return x * x * x * exp(-x);
}
static double GaussLegendre(
std::vector<double> x, std::vector<double> w,
double (*fx)(double), int n) {
double integral = 0.0;
for (int i = 0; i < n; i++) {
integral += fx(x[i]) * w[i];
}
return integral;
}
static double Integrate(double (*fx)(double), int n,
unsigned int seed) {
std::string msg;
std::vector<double> x(n);
std::vector<double> weights(n);
Mueller::OrthogonalPolynomialRoots(
"Legendre", "Zeros", n, n, seed, msg, x);
for (int i = 0; i < n; i++) {
double xi = x[i];
double dpn = Mueller::g("Legendre", "Extrema",
xi, n, 0, x);
weights[i] = 2.0 / ((1.0 - xi * xi) * dpn * dpn);
}
SelectionSort(x, weights, n);
std::cout << "Abscissas (x) and Weights (w)" << std::endl;
std::cout << "for Gauss-Legendre Quadrature";
std::cout << std::endl;
std::cout << "x\t\t\tw" << std::endl;
for (int i = 0; i < n; i++) {
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << x[i] << '\t';
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << weights[i] << std::endl;
}
return GaussLegendre(x, weights, fx, n);
}
static double SimpsonsRule(
int n, double a, double b,
double (*fx)(double))
{
double h = (b - a) / n;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < n; i += 2)
{
s += fx(x);
x += h2;
}
x = a + h2;
for (int i = 2; i < n; i += 2)
{
t += fx(x);
x += h2;
}
return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}
int main() {
unsigned int seed = 1;
for (int n = 10; n <= 40; n += 10)
{
double integral = Integrate(fx, n, seed++);
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << n << '\t' << integral << std::endl;
}
for (int n = 50000; n <= 100000; n += 10000)
{
double integral = SimpsonsRule(n, -1, 1, fx);
std::cout << "Simpson ";
std::cout << n << "\t\t";
std::cout << std::fixed << std::setprecision(15);
std::cout << std::setw(18);
std::cout << integral << std::endl;
}
return 0;
}
Blog Entry © Sunday, March 22, 2026, by James Pate Williams, Jr. Mueller’s Method
/*
* MuellersMethod.c (c) Sunday, July 21, 2024 by
* James Pate Williams, Jr. BA, BS, MSwE, PhD
* Translated from the FORTRAN 77 source code
* found in "Elementary Numerical Analysis: An
* Algorithmic Approach" by S. D. Conte and Carl
* de Boor Originally coded in FORTRAN IV in 1982 then
* into C# in March 2015 Finished Tuesday,
* July 23, 2024 The complex division method is
* from "A Numerical Library in C for Scientists
* and Engineers" by H. T. Lau Chapter 1
*/
#include <complex.h>
#include <math.h>
#include <stdio.h>
#include <stdlib.h>
//#define DEBUG
static _Lcomplex HornersMethod(_Lcomplex coeff[], _Lcomplex z, int degree) {
int i = 0;
_Lcomplex c = coeff[degree];
for (i = degree; i >= 1; i--) {
_Lcomplex product = _LCmulcc(c, z);
c._Val[0] = product._Val[0] + coeff[i - 1]._Val[0];
c._Val[1] = product._Val[1] + coeff[i - 1]._Val[1];
}
#ifdef DEBUG
_Lcomplex sum = { 0 };
for (i = 0; i <= degree; i++) {
_Lcomplex term =
_LCmulcc(cpowl(z, _LCbuild(i, 0.0)), coeff[i]);
sum._Val[0] += term._Val[0];
sum._Val[1] += term._Val[1];
}
long double delta = fabsl(cabsl(c) - cabsl(sum));
if (delta > 1.0e-12)
exit(-5);
#endif
return c;
}
static void comdiv(
long double xr, long double xi,
long double yr, long double yi,
long double* zr, long double* zi)
{
long double h, d;
if (fabs(yi) < fabs(yr)) {
if (yi == 0.0) {
*zr = xr / yr;
*zi = xi / yr;
}
else {
h = yi / yr;
d = h * yi + yr;
*zr = (xr + h * xi) / d;
*zi = (xi - h * xr) / d;
}
}
else {
h = yr / yi;
d = h * yr + yi;
*zr = (xr * h + xi) / d;
*zi = (xi * h - xr) / d;
}
}
#ifdef DEBUG
static _Lcomplex MyComplexDivide(_Lcomplex numer, _Lcomplex denom) {
long double norm2 =
denom._Val[0] * denom._Val[0] +
denom._Val[1] * denom._Val[1];
_Lcomplex result = { 0 };
result._Val[0] = (
numer._Val[0] * denom._Val[0] +
numer._Val[1] * denom._Val[1]) / norm2;
result._Val[1] = (
numer._Val[1] * denom._Val[0] -
numer._Val[0] * denom._Val[1]) / norm2;
return result;
}
#endif
static _Lcomplex ComplexDivide(_Lcomplex numer, _Lcomplex denom) {
_Lcomplex result = { 0 };
comdiv(
numer._Val[0], numer._Val[1],
denom._Val[0], denom._Val[1],
&result._Val[0], &result._Val[1]);
#ifdef DEBUG
_Lcomplex myResult = MyComplexDivide(numer, denom);
long double delta = fabsl(cabsl(result) - cabsl(myResult));
if (delta > 1.0e-12)
exit(-6);
#endif
return result;
}
static int Deflate(
_Lcomplex coeff[], _Lcomplex zero,
_Lcomplex* fzero, _Lcomplex* fzrdfl,
_Lcomplex zeros[], int i, int* count,
int degree) {
_Lcomplex denom = { 0 };
(*count)++;
*fzero = HornersMethod(coeff, zero, degree);
*fzrdfl = *fzero;
if (i < 1)
return 0;
for (int j = 1; j <= i; j++) {
denom._Val[0] = zero._Val[0] - zeros[j - 1]._Val[0];
denom._Val[1] = zero._Val[1] - zeros[j - 1]._Val[1];
if (cabsl(denom) == 0.0) {
zeros[i] = _LCmulcr(zero, 1.001);
return 1;
}
else
*fzrdfl = ComplexDivide(*fzrdfl, denom);
}
return 0;
}
static void Mueller(
_Lcomplex coeff[], _Lcomplex zeros[],
double epsilon1, double epsilon2,
int degree, int fnreal, int maxIts, int n, int nPrev) {
double eps1 = max(epsilon1, 1.0e-12);
double eps2 = max(epsilon2, 1.0e-20);
int count = 0, i = 0;
_Lcomplex c = { 0 };
_Lcomplex denom = { 0 };
_Lcomplex divdf1 = { 0 };
_Lcomplex divdf2 = { 0 };
_Lcomplex divdf1p = { 0 };
_Lcomplex fzero = { 0 };
_Lcomplex fzr = { 0 };
_Lcomplex fzdfl = { 0 };
_Lcomplex fzrdfl = { 0 };
_Lcomplex fzrprv = { 0 };
_Lcomplex four = _LCbuild(4.0, 0.0);
_Lcomplex h = { 0 };
_Lcomplex hprev = { 0 };
_Lcomplex sqr = { 0 };
_Lcomplex zero = { 0 };
_Lcomplex p5 = _LCbuild(0.5, 0.0);
_Lcomplex zeropp5 = { 0 };
_Lcomplex zeromp5 = { 0 };
_Lcomplex diff = { 0 };
_Lcomplex tadd = { 0 };
_Lcomplex tmul = { 0 };
_Lcomplex umul = { 0 };
_Lcomplex vmul = { 0 };
for (i = nPrev; i < n; i++) {
count = 0;
Label1:
zero = zeros[i];
h = p5;
zeropp5._Val[0] = zero._Val[0] + p5._Val[0];
zeropp5._Val[1] = zero._Val[1] + p5._Val[1];
if (Deflate(
coeff, zeropp5, &fzr, &divdf1p,
zeros, i, &count, degree))
goto Label1;
zeromp5._Val[0] = zero._Val[0] - p5._Val[0];
zeromp5._Val[1] = zero._Val[1] - p5._Val[1];
if (Deflate(
coeff, zeromp5, &fzr, &fzrprv,
zeros, i, &count, degree))
goto Label1;
hprev._Val[0] = -1.0;
hprev._Val[1] = 0.0;
diff._Val[0] = fzrprv._Val[0] - divdf1p._Val[0];
diff._Val[1] = fzrprv._Val[1] - divdf1p._Val[1];
if (cabsl(hprev) == 0)
exit(-2);
divdf1p = ComplexDivide(diff, hprev);
if (Deflate(
coeff, zero, &fzr, &fzrdfl,
zeros, i, &count, degree))
goto Label1;
Label2:
diff._Val[0] = fzrdfl._Val[0] - fzrprv._Val[0];
diff._Val[1] = fzrdfl._Val[1] - fzrprv._Val[1];
if (cabsl(h) == 0)
exit(-3);
divdf1 = ComplexDivide(diff, h);
diff._Val[0] = divdf1._Val[0] - divdf1p._Val[0];
diff._Val[1] = divdf1._Val[1] - divdf1p._Val[1];
tadd._Val[0] = h._Val[0] + hprev._Val[0];
tadd._Val[1] = h._Val[1] + hprev._Val[1];
if (cabsl(tadd) == 0)
exit(-3);
divdf2 = ComplexDivide(diff, tadd);
hprev = h;
divdf1p = divdf1;
tmul = _LCmulcc(h, divdf2);
c._Val[0] = divdf1._Val[0] + tmul._Val[0];
c._Val[1] = divdf1._Val[1] + tmul._Val[1];
tmul = _LCmulcc(c, c);
umul = _LCmulcc(four, fzrdfl);
vmul = _LCmulcc(umul, divdf2);
sqr._Val[0] = tmul._Val[0] - vmul._Val[0];
sqr._Val[1] = tmul._Val[1] - vmul._Val[1];
if (fnreal && sqr._Val[0] < 0.0)
{
sqr._Val[0] = 0.0;
sqr._Val[1] = 0.0;
}
sqr = csqrtl(sqr);
if ((c._Val[0] * sqr._Val[0] + c._Val[1] * sqr._Val[1]) < 0.0) {
denom._Val[0] = c._Val[0] - sqr._Val[0];
denom._Val[1] = c._Val[1] - sqr._Val[1];
}
else {
denom._Val[0] = c._Val[0] + sqr._Val[0];
denom._Val[1] = c._Val[1] + sqr._Val[1];
}
if (cabsl(denom) <= 0.0)
{
denom._Val[0] = 1.0;
denom._Val[1] = 0.0;
}
if (cabsl(denom) == 0)
exit(-4);
tmul = _LCmulcr(fzrdfl, -2.0);
h = ComplexDivide(tmul, denom);
fzrprv = fzrdfl;
zero._Val[0] = zero._Val[0] + h._Val[0];
zero._Val[1] = zero._Val[1] + h._Val[1];
if (count > maxIts)
goto Label4;
Label3:
if (Deflate(
coeff, zero, &fzr, &fzrdfl,
zeros, i, &count, degree))
goto Label1;
if (cabsl(h) < eps1 * cabsl(zero))
goto Label4;
if (max(cabsl(fzr), cabsl(fzdfl)) < eps2)
goto Label4;
if (cabsl(fzrdfl) >= 10.0 * cabsl(fzrprv)) {
h = _LCmulcr(h, 0.5);
zero._Val[0] = zero._Val[0] - h._Val[0];
zero._Val[1] = zero._Val[1] - h._Val[1];
goto Label3;
}
else
goto Label2;
Label4:
zeros[i] = zero;
}
}
int main(void)
{
double epsilon1 = 1.0e-15;
double epsilon2 = 1.0e-15;
int degree = 0, fnreal = 0, i = 0, maxIts = 1000;
int n = 0, nPrev = 0;
while (1) {
_Lcomplex* coeff = NULL;
_Lcomplex* zeros = NULL;
printf_s("Degree (0 to quit) = ");
scanf_s("%d", °ree);
if (degree == 0)
break;
n = degree;
coeff = calloc(degree + 1, sizeof(_Lcomplex));
if (coeff == NULL)
exit(-1);
zeros = calloc(n, sizeof(_Lcomplex));
if (zeros == NULL)
exit(-1);
for (i = degree; i >= 0; i--) {
printf_s("coefficient[%d].real = ", i);
scanf_s("%Lf", &coeff[i]._Val[0]);
printf_s("coefficient[%d].imag = ", i);
scanf_s("%Lf", &coeff[i]._Val[1]);
}
Mueller(
coeff, zeros, epsilon1,
epsilon2, degree, fnreal,
maxIts, n, nPrev);
printf_s("\n");
for (i = 0; i < degree; i++) {
printf_s("zero[%d].real = %17.10e\t", i, zeros[i]._Val[0]);
printf_s("zero[%d].imag = %17.10e\n", i, zeros[i]._Val[1]);
}
printf_s("\n");
for (i = 0; i < degree; i++) {
_Lcomplex func = HornersMethod(coeff, zeros[i], degree);
printf_s("func[%d].real = %17.10e\t", i, func._Val[0]);
printf_s("func[%d].imag = %17.10e\n", i, func._Val[1]);
}
printf_s("\n");
free(coeff);
free(zeros);
}
return 0;
}
Blog Input © Saturday, March 21, 2026, by James Pate Williams, Jr., Three Model Universes
#pragma once
class UniversalModels
{
public:
static void deSitterUniverse(
double epsilon, double lambda, double ct,
double A, double B, double& K);
static void RadiationUniverse(
double A, double epsilon,
double deltaT, double& K);
static void FriedmannUniverse(
int epsilon, double M,
double cT, double& K);
};
#include "UniversalModels.h"
#include <math.h>
void UniversalModels::deSitterUniverse(
double epsilon, double lambda,
double ct, double A, double B, double& K)
{
if (lambda > 0)
{
if (epsilon == +1)
{
K = cosh(B * ct) / B;
}
else if (epsilon == -1)
{
K = sinh(B * ct) / B;
}
else
{
K = A * exp(B * ct);
}
}
else if (lambda < 0)
{
if (epsilon == -1)
{
lambda = -3.0 * B * B;
K = cos(B * ct) / B;
}
}
}
void UniversalModels::RadiationUniverse(
double A, double epsilon,
double deltaT, double& K)
{
double c = 1.0, kappa = 1.0;
if (epsilon == 0)
{
K = sqrt(2.0 * c * sqrt(kappa * A / 3.0) * deltaT);
}
else if (epsilon == -1)
{
K = sqrt(c * c * pow(deltaT, 2.0) + 2 * c *
sqrt(kappa * A / 3.0) * deltaT);
}
else if (epsilon == +1)
{
K = sqrt(-c * c * pow(deltaT, 2.0) + 2 * c *
sqrt(kappa * A / 3.0) * deltaT);
}
}
void UniversalModels::FriedmannUniverse(
int epsilon, double M, double cT, double& K)
{
if (epsilon == 0)
{
K = M * cT * cT / 12.0;
}
else if (epsilon == -1)
{
K = M * (cosh(cT) - 1.0) / 6.0;
}
else if (epsilon == +1)
{
K = M * (1.0 - cos(cT)) / 6.0;
}
}
// ModelUniverses.cpp (c) Thursday, March 19, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Reference: "General Relativity An Introduction
// to the theory of the gravitational field "
// (c) 1982 by Hans Stephani translated by Martin
// Pollack and John Stewart
#include "framework.h"
#include "ModelUniverses.h"
#include "UniversalModels.h"
#include <math.h>
#include <string>
#include <vector>
#define MAX_LOADSTRING 100
typedef struct tagPoint2dDeSitter
{
double ct, K;
} Point2dDeSitter, * PPoint2dDeSitter;
typedef struct tagPoint2dRadiation
{
double deltaT, K;
} Point2dRadiation, * PPoint2dRadiation;
typedef struct tagPoint2dFriedmann
{
double cT, K;
} Point2dFriedmann, * PPoint2dFriedmann;
// Global Variables:
HINSTANCE hInst; // current instance
WCHAR szTitle[MAX_LOADSTRING]; // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING]; // the main window class name
bool deSitter, radiation, friedmann; // universal model type
std::wstring fTitle, xTitle, yTitle; // graph titles
std::vector<Point2dDeSitter> pointsDeSitter; // graph points
std::vector<Point2dRadiation> pointsRadiation; // graph points
std::vector<Point2dFriedmann> pointsFriedmann; // graph points
// 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 DataInputDialogDeSitter(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DrawGraphDialogDeSitter(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DataInputDialogRadiation(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DrawGraphDialogRadiation(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DataInputDialogFriedmann(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK DrawGraphDialogFriedmann(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_MODELUNIVERSES, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_MODELUNIVERSES));
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_MODELUNIVERSES));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_MODELUNIVERSES);
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_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_DE_SITTER:
deSitter = true;
fTitle = L"de Sitter Universe";
xTitle = L"Bct";
yTitle = L"K";
DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG_DE_SITTER),
hWnd, DataInputDialogDeSitter);
break;
case IDM_RADIATION:
radiation = true;
fTitle = L"Radiation Universe";
xTitle = L"t - t0";
yTitle = L"K";
DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG_RADIATION), hWnd,
DataInputDialogRadiation);
break;
case IDM_FRIEDMANN:
friedmann = true;
fTitle = L"Friedmann Universe";
xTitle = L"cT";
yTitle = L"K";
DialogBox(hInst, MAKEINTRESOURCE(IDD_DATA_INPUT_DIALOG_FRIEDMANN), hWnd,
DataInputDialogFriedmann);
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 DataInputDialogDeSitter(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
static WCHAR data[128] = L"";
static HWND hCombo_A = GetDlgItem(hDlg, IDC_COMBO_A_1);
static HWND hCombo_B = GetDlgItem(hDlg, IDC_COMBO_B);
static HWND hCombo_ct = GetDlgItem(hDlg, IDC_COMBO_CT);
static HWND hCombo_epsilon = GetDlgItem(hDlg, IDC_COMBO_EPSILON_1);
static HWND hCombo_lambda = GetDlgItem(hDlg, IDC_COMBO_LAMBDA);
switch (message)
{
case WM_INITDIALOG:
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"1");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"9");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_ADDSTRING, 0, (LPARAM)L"10");
SendDlgItemMessage(hDlg, IDC_COMBO_A_1, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"1");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"9");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_ADDSTRING, 0, (LPARAM)L"10");
SendDlgItemMessage(hDlg, IDC_COMBO_B, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"1");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"9");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_ADDSTRING, 0, (LPARAM)L"10");
SendDlgItemMessage(hDlg, IDC_COMBO_CT, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_ADDSTRING, 0, (LPARAM)L"-1");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_ADDSTRING, 0, (LPARAM)L"0");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_ADDSTRING, 0, (LPARAM)L"+1");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_1, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-5");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-4");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-3");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-2");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"-1");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"1");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_LAMBDA, CB_SETCURSEL, 0, 0);
return (INT_PTR)TRUE;
case WM_COMMAND:
{
if (LOWORD(wParam) == IDC_BUTTON_DRAW_DE_SITTER)
{
char line[128] = "";
GetDlgItemTextA(hDlg, IDC_COMBO_A_1, line, 128);
double A = atof(line);
GetDlgItemTextA(hDlg, IDC_COMBO_B, line, 128);
double B = atof(line);
GetDlgItemTextA(hDlg, IDC_COMBO_CT, line, 128);
double CT = atof(line);
int epsilon = GetDlgItemInt(hDlg, IDC_COMBO_EPSILON_1, FALSE, TRUE);
GetDlgItemTextA(hDlg, IDC_COMBO_LAMBDA, line, 128);
double lambda = atof(line);
double ct = 0.0, K = 0.0;
double delta_ct = CT / 256.0;
int index = 0;
pointsDeSitter.clear();
while (ct <= CT)
{
UniversalModels::deSitterUniverse(
epsilon, lambda, ct, A, B, K);
Point2dDeSitter pt = { 0 };
pt.ct = ct;
pt.K = K;
ct += delta_ct;
pointsDeSitter.push_back(pt);
}
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG_DE_SITTER),
hDlg, DrawGraphDialogDeSitter);
return (INT_PTR)TRUE;
}
else if (LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
}
}
return (INT_PTR)FALSE;
}
static void FindMinMaxDeSitter(
double& xmin,
double& xmax,
double& ymin,
double& ymax)
{
xmin = pointsDeSitter[0].ct;
xmax = pointsDeSitter[0].ct;
ymin = pointsDeSitter[0].K;
ymax = pointsDeSitter[0].K;
for (int i = 1; i < pointsDeSitter.size(); i++)
{
Point2dDeSitter pt = pointsDeSitter[i];
if (pt.ct < xmin)
xmin = pt.ct;
if (pt.ct > xmax)
xmax = pt.ct;
if (pt.K < ymin)
ymin = pt.K;
if (pt.K > ymax)
ymax = pt.K;
}
}
static void DrawTitles(
HDC hdc, RECT clientRect,
const std::wstring& fTitle,
const std::wstring& xTitle,
const std::wstring& yTitle,
int sx0, int sx1, int sy0, int sy1)
{
HFONT hCustomFont = CreateFont(
-MulDiv(10, GetDeviceCaps(hdc, LOGPIXELSY), 72), // Height in logical units
0, // Width (0 = default)
0, // Escapement
0, // Orientation
FW_BOLD, // Weight (700 = bold)
FALSE, // Italic
FALSE, // Underline
FALSE, // StrikeOut
DEFAULT_CHARSET, // Charset
OUT_DEFAULT_PRECIS, // Output precision
CLIP_DEFAULT_PRECIS, // Clipping precision
DEFAULT_QUALITY, // Quality
FIXED_PITCH | FF_MODERN,// Pitch and family
L"Courier New" // Typeface name
);
HFONT hOldFont = (HFONT)SelectObject(hdc, hCustomFont);
SIZE sz;
int width = clientRect.right - clientRect.left;
// Title: fTitle
GetTextExtentPoint32W(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
int w = sz.cx;
int h = sz.cy;
TextOutW(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());
// draw x-axis title
GetTextExtentPoint32W(hdc, xTitle.c_str(), (int)xTitle.length(), &sz);
w = sz.cx;
TextOutW(hdc, sx0 + (sx1 - sx0 - w) / 2, sy1 + 2 * h,
xTitle.c_str(), (int)xTitle.length());
// draw y-axis title
GetTextExtentPoint32W(hdc, yTitle.c_str(), (int)yTitle.length(), &sz);
w = sz.cx;
TextOutW(hdc, sx1 + w / 5, sy0 + (sy1 - sy0) / 2 - h / 2,
yTitle.c_str(), (int)yTitle.length());
SelectObject(hdc, hOldFont); // Restore original font
}
static void DrawFormattedText(HDC hdc, char loctext[], RECT rect)
{
// Draw the text with formatting options
DrawTextA(hdc, loctext, (int)strlen(loctext), &rect, DT_SINGLELINE | DT_NOCLIP);
}
static INT_PTR CALLBACK DrawGraphDialogDeSitter(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
char line[256] = "";
switch (message)
{
case WM_INITDIALOG:
strcpy_s(line, 256, "Drawing Dialog deSitter Universe");
SetWindowTextA(hDlg, line);
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
case WM_PAINT:
double xmax = 0, xmin = 0;
double ymax = 0, ymin = 0;
FindMinMaxDeSitter(xmin, xmax, ymin, ymax);
double h = 0, pi = 0, plm = 0, theta = 0;
float xSpan = (float)(xmax - xmin);
float ySpan = (float)(ymax - ymin);
RECT rect = { };
GetClientRect(hDlg, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f;
float sx1 = 14.0f * width / 16.0f;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xmin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * ymax);
float px = 0, py = 0, sx = 0, sy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
int i = 0;
float x = (float)xmin;
float y = (float)ymax;
DrawTitles(hdc, rect, fTitle, xTitle, yTitle,
(int)sx0, (int)sx1, (int)sy0, (int)sy1);
px = (float)pointsDeSitter[0].ct;
py = (float)pointsDeSitter[0].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
while (i <= 8)
{
sx = xSlope * x + xInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
LineTo(hdc, (int)sx, (int)sy1);
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%5.4lf", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { };
textRect.left = (long)(sx - size.cx / 2.0f);
textRect.right = (long)(sx + size.cx / 2.0f);
textRect.top = (long)sy1;
textRect.bottom = (long)(sy1 + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
x += deltaX;
i++;
}
i = 0;
y = (float)ymin;
while (i <= 8)
{
sy = ySlope * y + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
LineTo(hdc, (int)sx, (int)sy);
if (i != 0)
{
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%+5.3lf", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { 0 };
textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
textRect.right = (long)(sx0 - size.cx / 2.0f);
textRect.top = (long)(sy - size.cy / 2.0f);
textRect.bottom = (long)(sy + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
}
y += deltaY;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
HRGN clipRegion = CreateRectRgn(
(int)sx0, (int)sy0, // Left, Top
(int)(sx1), (int)(sy1) // Right, Bottom
);
// Apply clipping region
SelectClipRgn(hdc, clipRegion);
px = (float)pointsDeSitter[0].ct;
py = (float)pointsDeSitter[0].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
for (size_t j = 1; j < pointsDeSitter.size(); j++)
{
px = (float)pointsDeSitter[j].ct;
py = (float)pointsDeSitter[j].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
DeleteObject(bPenNew);
bPenNew = NULL;
SelectObject(hdc, hPenOld);
DeleteObject(clipRegion);
EndPaint(hDlg, &ps);
return (INT_PTR)TRUE;
}
return (INT_PTR)FALSE;
}
INT_PTR CALLBACK DataInputDialogRadiation(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
static WCHAR data[128] = L"";
static HWND hCombo_B = GetDlgItem(hDlg, IDC_COMBO_B);
static HWND hCombo_ct = GetDlgItem(hDlg, IDC_COMBO_DELTA_T);
static HWND hCombo_epsilon = GetDlgItem(hDlg, IDC_COMBO_EPSILON_2);
static HWND hCombo_lambda = GetDlgItem(hDlg, IDC_COMBO_LAMBDA);
switch (message)
{
case WM_INITDIALOG:
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e1");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e2");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e3");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e4");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"1e5");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"9");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_ADDSTRING, 0, (LPARAM)L"10");
SendDlgItemMessage(hDlg, IDC_COMBO_A_2, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-5");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-4");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-3");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-2");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e-1");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+0");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+1");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+2");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+3");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+4");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_ADDSTRING, 0, (LPARAM)L"1.0e+5");
SendDlgItemMessage(hDlg, IDC_COMBO_DELTA_T, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_ADDSTRING, 0, (LPARAM)L"-1");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_ADDSTRING, 0, (LPARAM)L"0");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_ADDSTRING, 0, (LPARAM)L"+1");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_2, CB_SETCURSEL, 0, 0);
return (INT_PTR)TRUE;
case WM_COMMAND:
{
if (LOWORD(wParam) == IDC_BUTTON_DRAW_RADIATION)
{
char line[128] = "";
GetDlgItemTextA(hDlg, IDC_COMBO_A_2, line, 128);
double A = atof(line);
GetDlgItemTextA(hDlg, IDC_COMBO_DELTA_T, line, 128);
double deltaT = atof(line);
int epsilon = GetDlgItemInt(hDlg, IDC_COMBO_EPSILON_2, FALSE, TRUE);
double step = deltaT / 256.0, t = 0.0, K = 0.0;
int index = 0;
pointsRadiation.clear();
while (t <= deltaT)
{
UniversalModels::RadiationUniverse(
A, epsilon, t, K);
Point2dRadiation pt = { 0 };
pt.deltaT = t;
pt.K = K;
t += step;
pointsRadiation.push_back(pt);
}
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG_RADIATION),
hDlg, DrawGraphDialogRadiation);
return (INT_PTR)TRUE;
}
else if (LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
}
}
return (INT_PTR)FALSE;
}
static void FindMinMaxRadiation(
double& xmin,
double& xmax,
double& ymin,
double& ymax)
{
xmin = pointsRadiation[0].deltaT;
xmax = pointsRadiation[0].deltaT;
ymin = pointsRadiation[0].K;
ymax = pointsRadiation[0].K;
for (int i = 1; i < pointsRadiation.size(); i++)
{
Point2dRadiation pt = pointsRadiation[i];
if (pt.deltaT < xmin)
xmin = pt.deltaT;
if (pt.deltaT > xmax)
xmax = pt.deltaT;
if (pt.K < ymin)
ymin = pt.K;
if (pt.K > ymax)
ymax = pt.K;
}
}
static INT_PTR CALLBACK DrawGraphDialogRadiation(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
char line[256] = "";
switch (message)
{
case WM_INITDIALOG:
strcpy_s(line, 256, "Drawing Dialog Radiation Universe");
SetWindowTextA(hDlg, line);
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
case WM_PAINT:
double xmax = 0, xmin = 0;
double ymax = 0, ymin = 0;
FindMinMaxRadiation(xmin, xmax, ymin, ymax);
double h = 0, pi = 0, plm = 0, theta = 0;
float xSpan = (float)(xmax - xmin);
float ySpan = (float)(ymax - ymin);
RECT rect = { };
GetClientRect(hDlg, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f;
float sx1 = 14.0f * width / 16.0f;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xmin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * ymax);
float px = 0, py = 0, sx = 0, sy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
int i = 0;
float x = (float)xmin;
float y = (float)ymax;
DrawTitles(hdc, rect, fTitle, xTitle, yTitle,
(int)sx0, (int)sx1, (int)sy0, (int)sy1);
px = (float)pointsRadiation[0].deltaT;
py = (float)pointsRadiation[0].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
while (i <= 8)
{
sx = xSlope * x + xInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
LineTo(hdc, (int)sx, (int)sy1);
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%5.4lf", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { };
textRect.left = (long)(sx - size.cx / 2.0f);
textRect.right = (long)(sx + size.cx / 2.0f);
textRect.top = (long)sy1;
textRect.bottom = (long)(sy1 + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
x += deltaX;
i++;
}
i = 0;
y = (float)ymin;
while (i <= 8)
{
sy = ySlope * y + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
LineTo(hdc, (int)sx, (int)sy);
if (i != 0)
{
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%+5.3lf", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { 0 };
textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
textRect.right = (long)(sx0 - size.cx / 2.0f);
textRect.top = (long)(sy - size.cy / 2.0f);
textRect.bottom = (long)(sy + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
}
y += deltaY;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
HRGN clipRegion = CreateRectRgn(
(int)sx0, (int)sy0, // Left, Top
(int)(sx1), (int)(sy1) // Right, Bottom
);
// Apply clipping region
SelectClipRgn(hdc, clipRegion);
px = (float)pointsRadiation[0].deltaT;
py = (float)pointsRadiation[0].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
for (size_t j = 1; j < pointsRadiation.size(); j++)
{
px = (float)pointsRadiation[j].deltaT;
py = (float)pointsRadiation[j].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
DeleteObject(bPenNew);
bPenNew = NULL;
SelectObject(hdc, hPenOld);
DeleteObject(clipRegion);
EndPaint(hDlg, &ps);
return (INT_PTR)TRUE;
}
return (INT_PTR)FALSE;
}
INT_PTR CALLBACK DataInputDialogFriedmann(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
static WCHAR data[128] = L"";
static HWND hCombo_M = GetDlgItem(hDlg, IDC_COMBO_M);
static HWND hCombo_CT1 = GetDlgItem(hDlg, IDC_COMBO_CT1);
static HWND hCombo_epsilon = GetDlgItem(hDlg, IDC_COMBO_EPSILON_3);
switch (message)
{
case WM_INITDIALOG:
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"1");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"9");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_ADDSTRING, 0, (LPARAM)L"10");
SendDlgItemMessage(hDlg, IDC_COMBO_M, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"1");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"2");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"3");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"4");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"5");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"6");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"7");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"8");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"9");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_ADDSTRING, 0, (LPARAM)L"10");
SendDlgItemMessage(hDlg, IDC_COMBO_CT1, CB_SETCURSEL, 0, 0);
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_ADDSTRING, 0, (LPARAM)L"-1");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_ADDSTRING, 0, (LPARAM)L"0");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_ADDSTRING, 0, (LPARAM)L"+1");
SendDlgItemMessage(hDlg, IDC_COMBO_EPSILON_3, CB_SETCURSEL, 0, 0);
return (INT_PTR)TRUE;
case WM_COMMAND:
{
if (LOWORD(wParam) == IDC_BUTTON_DRAW_FRIEDMANN)
{
char line[128] = "";
GetDlgItemTextA(hDlg, IDC_COMBO_M, line, 128);
double M = atof(line);
GetDlgItemTextA(hDlg, IDC_COMBO_CT1, line, 128);
double cT = atof(line);
int epsilon = GetDlgItemInt(hDlg, IDC_COMBO_EPSILON_3, FALSE, TRUE);
double step = cT / 256.0, t = 0.0, K = 0.0;
int index = 0;
pointsFriedmann.clear();
while (t <= cT)
{
UniversalModels::FriedmannUniverse(
epsilon, M, t, K);
Point2dFriedmann pt = { 0 };
pt.cT = t;
pt.K = K;
t += step;
pointsFriedmann.push_back(pt);
}
DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_GRAPH_DIALOG_FRIEDMANN),
hDlg, DrawGraphDialogFriedmann);
return (INT_PTR)TRUE;
}
else if (LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
}
}
return (INT_PTR)FALSE;
}
static void FindMinMaxFriedmann(
double& xmin,
double& xmax,
double& ymin,
double& ymax)
{
xmin = pointsFriedmann[0].cT;
xmax = pointsFriedmann[0].cT;
ymin = pointsFriedmann[0].K;
ymax = pointsFriedmann[0].K;
for (int i = 1; i < pointsFriedmann.size(); i++)
{
Point2dFriedmann pt = pointsFriedmann[i];
if (pt.cT < xmin)
xmin = pt.cT;
if (pt.cT > xmax)
xmax = pt.cT;
if (pt.K < ymin)
ymin = pt.K;
if (pt.K > ymax)
ymax = pt.K;
}
}
static INT_PTR CALLBACK DrawGraphDialogFriedmann(
HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
UNREFERENCED_PARAMETER(lParam);
char line[256] = "";
switch (message)
{
case WM_INITDIALOG:
strcpy_s(line, 256, "Drawing Dialog Friedmann Universe");
SetWindowTextA(hDlg, line);
return (INT_PTR)TRUE;
case WM_COMMAND:
if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
{
EndDialog(hDlg, LOWORD(wParam));
return (INT_PTR)TRUE;
}
break;
case WM_PAINT:
double xmax = 0, xmin = 0;
double ymax = 0, ymin = 0;
FindMinMaxFriedmann(xmin, xmax, ymin, ymax);
double h = 0, pi = 0, plm = 0, theta = 0;
float xSpan = (float)(xmax - xmin);
float ySpan = (float)(ymax - ymin);
RECT rect = { };
GetClientRect(hDlg, &rect);
float width = (float)(rect.right - rect.left + 1);
float height = (float)(rect.bottom - rect.top - 32 + 1);
float sx0 = 2.0f * width / 16.0f;
float sx1 = 14.0f * width / 16.0f;
float sy0 = 2.0f * height / 16.0f;
float sy1 = 14.0f * height / 16.0f;
float deltaX = xSpan / 8.0f;
float deltaY = ySpan / 8.0f;
float xSlope = (sx1 - sx0) / xSpan;
float xInter = (float)(sx0 - xSlope * xmin);
float ySlope = (sy0 - sy1) / ySpan;
float yInter = (float)(sy0 - ySlope * ymax);
float px = 0, py = 0, sx = 0, sy = 0;
PAINTSTRUCT ps;
POINT wPt = { };
HDC hdc = BeginPaint(hDlg, &ps);
int i = 0;
float x = (float)xmin;
float y = (float)ymax;
DrawTitles(hdc, rect, fTitle, xTitle, yTitle,
(int)sx0, (int)sx1, (int)sy0, (int)sy1);
px = (float)pointsFriedmann[0].cT;
py = (float)pointsFriedmann[0].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
while (i <= 8)
{
sx = xSlope * x + xInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
LineTo(hdc, (int)sx, (int)sy1);
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%5.4lf", x);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { };
textRect.left = (long)(sx - size.cx / 2.0f);
textRect.right = (long)(sx + size.cx / 2.0f);
textRect.top = (long)sy1;
textRect.bottom = (long)(sy1 + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
x += deltaX;
i++;
}
i = 0;
y = (float)ymin;
while (i <= 8)
{
sy = ySlope * y + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx0, (int)sy, &wPt);
LineTo(hdc, (int)sx, (int)sy);
if (i != 0)
{
char numberStr[128] = "";
sprintf_s(numberStr, 128, "%+5.3lf", y);
SIZE size = { };
GetTextExtentPoint32A(
hdc,
numberStr,
(int)strlen(numberStr),
&size);
RECT textRect = { 0 };
textRect.left = (long)(sx0 - size.cx - size.cx / 2.0f);
textRect.right = (long)(sx0 - size.cx / 2.0f);
textRect.top = (long)(sy - size.cy / 2.0f);
textRect.bottom = (long)(sy + size.cy / 2.0f);
DrawFormattedText(hdc, numberStr, textRect);
}
y += deltaY;
i++;
}
HGDIOBJ bPenNew = NULL;
HGDIOBJ hPenOld = NULL;
bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, bPenNew);
HRGN clipRegion = CreateRectRgn(
(int)sx0, (int)sy0, // Left, Top
(int)(sx1), (int)(sy1) // Right, Bottom
);
// Apply clipping region
SelectClipRgn(hdc, clipRegion);
px = (float)pointsFriedmann[0].cT;
py = (float)pointsFriedmann[0].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
wPt.x = wPt.y = 0;
MoveToEx(hdc, (int)sx, (int)sy, &wPt);
for (size_t j = 1; j < pointsFriedmann.size(); j++)
{
px = (float)pointsFriedmann[j].cT;
py = (float)pointsFriedmann[j].K;
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
DeleteObject(bPenNew);
bPenNew = NULL;
SelectObject(hdc, hPenOld);
DeleteObject(clipRegion);
EndPaint(hDlg, &ps);
return (INT_PTR)TRUE;
}
return (INT_PTR)FALSE;
}
//Microsoft Visual C++ generated resource script.
//
#include "resource.h"
#define APSTUDIO_READONLY_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
#ifndef APSTUDIO_INVOKED
#include "targetver.h"
#endif
#define APSTUDIO_HIDDEN_SYMBOLS
#include "windows.h"
#undef APSTUDIO_HIDDEN_SYMBOLS
/////////////////////////////////////////////////////////////////////////////
#undef APSTUDIO_READONLY_SYMBOLS
#if !defined(AFX_RESOURCE_DLL) || defined(AFX_TARG_ENU)
LANGUAGE 9, 1
/////////////////////////////////////////////////////////////////////////////
//
// Icon
//
// Icon with lowest ID value placed first to ensure application icon
// remains consistent on all systems.
IDI_MODELUNIVERSES ICON "ModelUniverses.ico"
IDI_SMALL ICON "small.ico"
/////////////////////////////////////////////////////////////////////////////
//
// Menu
//
IDC_MODELUNIVERSES MENU
BEGIN
POPUP "&Begin"
BEGIN
MENUITEM "&de Sitter Universe", IDM_DE_SITTER
MENUITEM "&Radiation Universe", IDM_RADIATION
MENUITEM "&Friedmann Universe", IDM_FRIEDMANN
MENUITEM SEPARATOR
MENUITEM "E&xit", IDM_EXIT
END
POPUP "&Help"
BEGIN
MENUITEM "&About ...", IDM_ABOUT
END
END
/////////////////////////////////////////////////////////////////////////////
//
// Accelerator
//
IDC_MODELUNIVERSES ACCELERATORS
BEGIN
"?", IDM_ABOUT, ASCII, ALT
"/", IDM_ABOUT, ASCII, ALT
END
/////////////////////////////////////////////////////////////////////////////
//
// Dialog
//
IDD_ABOUTBOX DIALOGEX 0, 0, 170, 62
STYLE DS_SETFONT | DS_MODALFRAME | DS_FIXEDSYS | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "About ModelUniverses"
FONT 8, "MS Shell Dlg"
BEGIN
ICON IDI_MODELUNIVERSES,IDC_STATIC,14,14,21,20
LTEXT "ModelUniverses, Version 1.0",IDC_STATIC,42,14,114,8,SS_NOPREFIX
LTEXT "Copyright (c) 2026",IDC_STATIC,42,26,114,8
DEFPUSHBUTTON "OK",IDOK,113,41,50,14,WS_GROUP
END
IDD_DATA_INPUT_DIALOG_DE_SITTER DIALOGEX 0, 0, 280, 200
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog de Sitter Universe"
FONT 10, "Courier New", 700
BEGIN
LTEXT "A:", IDC_STATIC, 10, 10, 100, 10
COMBOBOX IDC_COMBO_A_1, 120, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "B:", IDC_STATIC, 10, 40, 100, 10
COMBOBOX IDC_COMBO_B, 120, 40, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "ct:", IDC_STATIC, 10, 70, 100, 10
COMBOBOX IDC_COMBO_CT, 120, 70, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "epsilon:",
IDC_STATIC, 10, 100, 50, 60
COMBOBOX IDC_COMBO_EPSILON_1, 120, 100, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "Lambda:",
IDC_STATIC, 10, 130, 50, 60
COMBOBOX IDC_COMBO_LAMBDA, 120, 130, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
PUSHBUTTON "Draw", IDC_BUTTON_DRAW_DE_SITTER, 10, 160, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 210, 160, 50, 16
END
IDD_DRAW_GRAPH_DIALOG_DE_SITTER DIALOGEX 0, 0, 600, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog de Sitter Universe"
FONT 10, "Courier New", 700
BEGIN
END
IDD_DATA_INPUT_DIALOG_RADIATION DIALOGEX 0, 0, 280, 170
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog Radiation Universe"
FONT 10, "Courier New", 700
BEGIN
LTEXT "A:", IDC_STATIC, 10, 10, 100, 10
COMBOBOX IDC_COMBO_A_2, 120, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "t - t0:", IDC_STATIC, 10, 40, 100, 10
COMBOBOX IDC_COMBO_DELTA_T, 120, 40, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "epsilon:",
IDC_STATIC, 10, 70, 50, 60
COMBOBOX IDC_COMBO_EPSILON_2, 120, 70, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
PUSHBUTTON "Draw", IDC_BUTTON_DRAW_RADIATION, 10, 100, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 210, 100, 50, 16
END
IDD_DRAW_GRAPH_DIALOG_RADIATION DIALOGEX 0, 0, 600, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog Radiation Universe"
FONT 10, "Courier New", 700
BEGIN
END
IDD_DATA_INPUT_DIALOG_FRIEDMANN DIALOGEX 0, 0, 280, 170
STYLE DS_SETFONT | WS_POPUP | WS_CAPTION | WS_SYSMENU
CAPTION "Data Input Dialog Friedmann Universe"
FONT 10, "Courier New", 700
BEGIN
LTEXT "M:", IDC_STATIC, 10, 10, 100, 10
COMBOBOX IDC_COMBO_M, 120, 10, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "cT:", IDC_STATIC, 10, 40, 100, 10
COMBOBOX IDC_COMBO_CT1, 120, 40, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
LTEXT "epsilon:",
IDC_STATIC, 10, 70, 50, 60
COMBOBOX IDC_COMBO_EPSILON_3, 120, 70, 50, 60, CBS_DROPDOWNLIST
| WS_VSCROLL | WS_TABSTOP
PUSHBUTTON "Draw", IDC_BUTTON_DRAW_FRIEDMANN, 10, 100, 50, 16
PUSHBUTTON "Cancel", IDCANCEL, 210, 100, 50, 16
END
IDD_DRAW_GRAPH_DIALOG_FRIEDMANN DIALOGEX 0, 0, 600, 310
STYLE DS_SETFONT | WS_POPUP | WS_VISIBLE | WS_CAPTION | WS_SYSMENU
CAPTION "Draw Graph Dialog Friedmann Universe"
FONT 10, "Courier New", 700
BEGIN
END
/////////////////////////////////////////////////////////////////////////////
//
// DESIGNINFO
//
#ifdef APSTUDIO_INVOKED
GUIDELINES DESIGNINFO
BEGIN
IDD_ABOUTBOX, DIALOG
BEGIN
LEFTMARGIN, 7
RIGHTMARGIN, 163
TOPMARGIN, 7
BOTTOMMARGIN, 55
END
END
#endif // APSTUDIO_INVOKED
#ifdef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// TEXTINCLUDE
//
1 TEXTINCLUDE
BEGIN
"resource.h\0"
END
2 TEXTINCLUDE
BEGIN
"#ifndef APSTUDIO_INVOKED\r\n"
"#include ""targetver.h""\r\n"
"#endif\r\n"
"#define APSTUDIO_HIDDEN_SYMBOLS\r\n"
"#include ""windows.h""\r\n"
"#undef APSTUDIO_HIDDEN_SYMBOLS\r\n"
"\0"
END
3 TEXTINCLUDE
BEGIN
"\r\n"
"\0"
END
#endif // APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// String Table
//
STRINGTABLE
BEGIN
IDC_MODELUNIVERSES "MODELUNIVERSES"
IDS_APP_TITLE "ModelUniverses"
END
#endif
/////////////////////////////////////////////////////////////////////////////
#ifndef APSTUDIO_INVOKED
/////////////////////////////////////////////////////////////////////////////
//
// Generated from the TEXTINCLUDE resource.
//
/////////////////////////////////////////////////////////////////////////////
#endif // not APSTUDIO_INVOKED
//{{NO_DEPENDENCIES}}
// Microsoft Visual C++ generated include file.
// Used by ModelUniverses.rc
#define IDS_APP_TITLE 103
#define IDR_MAINFRAME 128
#define IDD_MODELUNIVERSES_DIALOG 102
#define IDD_ABOUTBOX 103
#define IDM_ABOUT 104
#define IDM_EXIT 105
#define IDI_MODELUNIVERSES 107
#define IDI_SMALL 108
#define IDC_MODELUNIVERSES 109
#define IDC_MYICON 2
#ifndef IDC_STATIC
#define IDC_STATIC -1
#endif
// Next default values for new objects
//
#ifdef APSTUDIO_INVOKED
#ifndef APSTUDIO_READONLY_SYMBOLS
#define _APS_NO_MFC 130
#define _APS_NEXT_RESOURCE_VALUE 129
#define _APS_NEXT_COMMAND_VALUE 32771
#define _APS_NEXT_CONTROL_VALUE 1000
#define _APS_NEXT_SYMED_VALUE 110
#endif
#endif
#define IDD_DATA_INPUT_DIALOG_DE_SITTER 1000
#define IDC_COMBO_A_1 1010
#define IDC_COMBO_B 1020
#define IDC_COMBO_CT 1030
#define IDC_COMBO_EPSILON_1 1040
#define IDC_COMBO_LAMBDA 1050
#define IDC_BUTTON_DRAW_DE_SITTER 1060
#define IDD_DRAW_GRAPH_DIALOG_DE_SITTER 2000
#define IDD_DATA_INPUT_DIALOG_RADIATION 3000
#define IDC_COMBO_A_2 3010
#define IDC_COMBO_DELTA_T 3020
#define IDC_COMBO_EPSILON_2 3030
#define IDC_BUTTON_DRAW_RADIATION 3040
#define IDD_DRAW_GRAPH_DIALOG_RADIATION 4000
#define IDD_DATA_INPUT_DIALOG_FRIEDMANN 5000
#define IDC_COMBO_M 5010
#define IDC_COMBO_CT1 5020
#define IDC_COMBO_EPSILON_3 5030
#define IDC_BUTTON_DRAW_FRIEDMANN 5040
#define IDD_DRAW_GRAPH_DIALOG_FRIEDMANN 6000
#define IDM_DE_SITTER 32771
#define IDM_RADIATION 32772
#define IDM_FRIEDMANN 32773
Blog Entry © Thursday, March 19, 2026, by James Pate Williams, Jr. Solution of Laplace’s Two-Dimensional Equation on a Square Boundary
#pragma once
class Laplace
{
private:
static int N;
public:
Laplace(int N);
double** Solve(
double L,
double epsilon,
int& iterations,
double (*f)(double));
void PrintSolution(
double L,
double epsilon,
double** u,
double** w,
int& iterations);
};
#include "Laplace.h"
#include <float.h>
#include <math.h>
#include <iomanip>
#include <iostream>
int Laplace::N = 0;
static double fxy(double x, double y)
{
return x * y * exp(-0.5 * (x * x + y * y));
}
static double gxy(double x, double y)
{
return x * sin(y);
}
Laplace::Laplace(int N)
{
this->N = N;
}
double** Laplace::Solve(
double L, double epsilon, int& iterations,
double (*f)(double))
{
double h = L / (N + 1.0), uv = 0, max = 0;
double** u = new double*[N + 2];
if (u == nullptr)
exit(-1);
for (int i = 0; i < N + 2; i++)
{
u[i] = new double[N + 2];
if (u[i] == nullptr)
exit(-2);
for (int j = 0; j < N + 2; j++)
u[i][j] = 0.0;
}
double** v = new double* [N + 2];
if (v == nullptr)
exit(-3);
for (int i = 0; i < N + 2; i++)
{
v[i] = new double[N + 2];
if (v[i] == nullptr)
exit(-4);
for (int j = 0; j < N + 2; j++)
v[i][j] = 0.0;
}
iterations = 0;
for (int j = 0; j <= N + 1; j++)
u[0][j] = 0.0;
for (int j = 0; j <= N + 1; j++)
u[N + 1][j] = 0.0;
for (int i = 0; i <= N + 1; i++)
u[i][0] = f(i * h);
for (int i = 0; i <= N + 1; i++)
u[i][N + 1] = 0.0;
do
{
iterations++;
for (int i = 1; i <= N; i++)
for (int j = 1; j <= N; j++)
u[i][j] = 0.25 * (u[i - 1][j] + u[i + 1][j] + u[i][j - 1] +
u[i][j + 1]);
for (int i = 0; i <= N + 1; i++)
for (int j = 0; j <= N + 1; j++)
v[i][j] = u[i][j];
for (int i = 1; i <= N; i++)
for (int j = 1; j <= N; j++)
v[i][j] = 0.25 * (v[i - 1][j] + v[i + 1][j] + v[i][j - 1] +
v[i][j + 1]);
max = DBL_MIN;
for (int i = 1; i <= N; i++)
{
for (int j = 1; j <= N; j++)
{
uv = fabs(v[i][j] - u[i][j]);
max = fmax(uv, max);
}
}
for (int i = 0; i <= N + 1; i++)
for (int j = 0; j <= N + 1; j++)
u[i][j] = v[i][j];
} while (max > epsilon);
for (int i = 0; i < N + 2; i++)
delete v[i];
return u;
}
void Laplace::PrintSolution(
double L,
double epsilon,
double** u,
double** w,
int& iterations)
{
double h = L / (N + 1.0);
std::cout << std::endl;
std::cout << "x\t\ty\t\tu(x, y)\t\tw(x, y)\t\tpercent error\r\n\r\n";
for (int i = 0; i <= N + 1; i++)
{
double x = i * h;
for (int j = 0; j <= N + 1; j++)
{
double y = j * h;
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8);
std::cout << x << '\t';
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8);
std::cout << y << '\t';
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8);
std::cout << u[i][j] << '\t';
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8);
std::cout << w[i][j] << '\t';
double error = 0.0;
if (w[i][j] != 0.0)
error = 100.0 * fabs(u[i][j] - w[i][j]) /
fabs(w[i][j]);
std::cout << std::fixed << std::setprecision(5);
std::cout << std::setw(8);
std::cout << error << std::endl;
}
}
std::cout << std::endl;
std::cout << "iterations = " << iterations;
std::cout << std::endl;
}
// LaplaceEquationFD2d.cpp (c) Tuesday, March 17, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// by James Pate Williams, Jr., BA, ,BS, MSwE, PhD
// Reference: Numerical Computation 2 Methods,
// Software And Analysis (c) 1995 by Chrisoph E. Ueberhuber
// Pages 392 - 393 Poisson Matrices
// https://math.libretexts.org/Bookshelves/Differential_Equations/Introduction_to_Partial_Differential_Equations_(Herman)/06%3A_Problems_in_Higher_Dimensions/6.03%3A_Laplaces_Equation_in_2D
#include <math.h>
#include <iostream>
#include "Laplace.h"
double** w;
double SimpsonsRule(double L, double a, double b,
int Nsr, int n, double (*f)(double, double, int))
{
double h = (b - a) / Nsr;
double h2 = 2.0 * h;
double s = 0.0;
double t = 0.0;
double x = a + h;
for (int i = 1; i < Nsr; i += 2)
{
s += f(L, x, n);
x += h2;
}
x = a + h2;
for (int i = 2; i < Nsr; i += 2)
{
t += f(L, x, n);
x += h2;
}
return
h * (f(L, a, n) + 4 * s + 2 * t + f(L, b, n)) / 3.0;
}
double f(double x)
{
return x * exp(-x) * cos(x);
}
double g(double L, double x, int n)
{
double pi = 4.0 * atan(1.0);
double fx = f(x);
return fx * sin(n * pi * x / L);
}
double bn(double L, int Nsr, int n)
{
double cs = 2.0 / L;
double sr = SimpsonsRule(L, 0.0, L, Nsr, n, g);
return cs * sr;
}
void ComputeW(double L, int N)
{
double h = L / (N + 1.0);
w = new double* [N + 2];
if (w == nullptr)
exit(-5);
for (int i = 0; i < N + 2; i++)
{
w[i] = new double[N + 2];
if (w[i] == nullptr)
exit(-6);
for (int j = 0; j < N + 2; j++)
w[i][j] = 0.0;
}
for (int j = 0; j <= N; j++)
w[0][j] = 0.0;
for (int j = 0; j <= N; j++)
w[N + 1][j] = 0.0;
for (int i = 0; i <= N; i++)
w[i][0] = f(i * h);
for (int i = 0; i <= N; i++)
w[i][N + 1] = 0.0;
double pi = 4.0 * atan(1.0);
double x = h, y = 0.0;
for (int i = 1; i <= N; i++)
{
y = h;
for (int j = 1; j <= N; j++)
{
double s = 0.0;
for (int n = 1; n <= 64; n++)
{
double npi = n * pi;
double an = bn(L, 512, n) / sinh(npi);
s += an * sin(npi * x / L) * sinh(npi * (L - y) / L);
}
w[i][j] = s;
y += h;
}
x += h;
}
}
int main()
{
while (true)
{
char line[128] = "";
std::cout << "N = ";
std::cin.getline(line, 128);
int N = atoi(line);
if (N < 4 || N > 128)
{
std::cout << "N < 4 || N > 128, try again";
std::cout << std::endl;
}
std::cout << "epsilon = ";
std::cin.getline(line, 128);
double epsilon = atof(line);
Laplace l(N);
double L = 1.0;
int iterations = 0;
ComputeW(L, N);
double** u = l.Solve(L, epsilon, iterations, f);
l.PrintSolution(L, epsilon, u, w, iterations);
for (int i = 0; i < N + 2; i++)
{
if (u[i] != nullptr)
delete u[i];
if (w[i] != nullptr)
delete w[i];
}
break;
}
return 0;
}
Blog Entry © Tuesday, March 17, 2026, by James Pate Williams, Jr., Comparison of Power Series Arctangent and Arcsine Functions with the C++ Built-In Functions
#pragma once
//https://scipp-legacy.pbsci.ucsc.edu/~haber/ph116A/taylor11.pdf
class Functions
{
public:
static void initialize();
static double factorial(int n);
static double arccosine(double x, int n);
static double arccosecant(double x);
static double arccotangent(double x);
static double arcsecant(double x);
static double arcsine(double x, int n);
static double arctangent(double x, int n);
};
#include "Functions.h"
#include <math.h>
double pi2 = 0.0;
double Functions::factorial(int n)
{
double nfactorial = 1.0;
for (int i = 2; i <= n; i++)
nfactorial *= i;
return nfactorial;
}
void Functions::initialize()
{
pi2 = 2.0 * atan(1.0);
}
double Functions::arccosine(double x, int n)
{
double sum = pi2 - arcsine(x, n);
return sum;
}
double Functions::arccosecant(double x)
{
double sum = 0;
return sum;
}
double Functions::arccotangent(double x)
{
double sum = 0;
return sum;
}
double Functions::arcsecant(double x)
{
double sum = 0;
return sum;
}
double Functions::arcsine(double x, int n)
{
double sum = 0.0;
if (fabs(x) <= 1.0)
{
for (int i = n; i >= 0; i--)
{
double ifact = factorial(i);
double i2 = 2.0 * i, i21 = 2 * i + 1;
double coeff = factorial(i2) /
(pow(2, i2) * ifact * ifact * (i21));
sum += coeff * pow(x, i21);
}
}
return sum;
}
double Functions::arctangent(double x, int n)
{
double sum = 0.0;
if (fabs(x) <= 1.0)
{
double one = 0.0;
if (n % 2 == 0)
one = 1.0;
else
one = -1.0;
for (int i = 2 * n + 1; i >= 0; i--)
{
one = -one;
double i21 = 2 * i + 1.0;
sum += one * pow(x, i21) / i21;
}
}
return sum;
}
// InvTrigConsoleCPP.cpp (c) Monday, March 16, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
#include "Functions.h"
#include <math.h>
#include <iomanip>
#include <iostream>
static void CreateTable(
char title[],
double a, double b, int n, int nPts,
double(*fx)(double, int), double(*gx)(double))
{
double x = a;
double deltaX = (b - a) / nPts;
std::cout << title << std::endl;
std::cout << "# Terms = " << (2 * n + 2) << std::endl;
std::cout << "x" << '\t' << "fx" << "\t\t";
std::cout << "sx" << "\t\t" << "error" << std::endl;
for (int i = 0; i < nPts; i++)
{
double f = fx(x, n);
double s = gx(x);
double e = 0.0;
if (fabs(s) != 0)
e = 100.0 * fabs(f - s) / fabs(s);
std::cout << std::fixed << std::setprecision(4);
std::cout << std::setw(4);
std::cout << x << '\t';
std::cout << std::fixed << std::setprecision(11);
std::cout << std::setw(12);
std::cout << f << '\t';
std::cout << std::fixed << std::setprecision(11);
std::cout << std::setw(12);
std::cout << s << '\t';
std::cout << e << std::endl;
x += deltaX;
}
}
int main()
{
Functions::initialize();
while (true)
{
char line[128] = "";
std::cout << "== Menu == " << std::endl;
std::cout << "1 arcsine" << std::endl;
std::cout << "2 arctangent" << std::endl;
std::cout << "3 exit" << std::endl;
std::cout << "option # : ";
std::cin.getline(line, 128);
char option = line[0];
if (option == '3')
break;
if (option < '1' || option > '2')
{
std::cout << "invalid option" << std::endl;
continue;
}
std::cout << "# points: ";
std::cin.getline(line, 128);
int nPts = atoi(line);
std::cout << "# terms: ";
std::cin.getline(line, 128);
int n = atoi(line);
if (option == '1')
{
char title[] = "Series arcsin Versus C++ asin";
CreateTable(title, 0.0, 1.0, n, nPts,
Functions::arcsine, asin);
}
else if (option == '2')
{
char title[] = "Series arctan Versus C++ atan";
CreateTable(title, 0.0, 1.0, n, nPts,
Functions::arctangent, atan);
}
}
return 0;
}
Numerical Explorations 