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;
}
Blog Entry © Sunday March 15, 2026, by James Pate Williams, Jr. Graphing Integer Order Bessel Functions of the First Kind and the real Gamma Function
Blog Entry © Saturday, March 14, 2026, by James Pate Williams, Jr. Power Series Solution of a Second Order Linear Ordinary Differential Equation
// SecondOrderLinearODE.cpp (c) Friday March 13, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solves the Bessel Equation of the First Kind for Order 0
// Using the power series found in the reference below -
// References: https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/17%3A_Second-Order_Differential_Equations/17.04%3A_Series_Solutions_of_Differential_Equations
// https://mathworld.wolfram.com/BesselDifferentialEquation.html
#include "framework.h"
#include <Windows.h>
#include "SecondOrderLinearODE.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 fTitle[128];
double chi[1025], eta[1025], coeff[65];
// Forward declarations of functions included in this code module:
ATOM MyRegisterClass(HINSTANCE hInstance);
BOOL InitInstance(HINSTANCE, int);
LRESULT CALLBACK WndProc(HWND, UINT, WPARAM, LPARAM);
INT_PTR CALLBACK About(HWND, UINT, WPARAM, LPARAM);
int APIENTRY 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_SECONDORDERLINEARODE, szWindowClass, MAX_LOADSTRING);
MyRegisterClass(hInstance);
// Perform application initialization:
if (!InitInstance (hInstance, nCmdShow))
{
return FALSE;
}
HACCEL hAccelTable = LoadAccelerators(hInstance, MAKEINTRESOURCE(IDC_SECONDORDERLINEARODE));
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_SECONDORDERLINEARODE));
wcex.hCursor = LoadCursor(nullptr, IDC_ARROW);
wcex.hbrBackground = (HBRUSH)(COLOR_WINDOW+1);
wcex.lpszMenuName = MAKEINTRESOURCEW(IDC_SECONDORDERLINEARODE);
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;
}
static double Factorial(int n)
{
double factorial = 1.0;
for (int i = 2; i <= n; i++)
factorial *= i;
return factorial;
}
static void ComputeCoefficients()
{
double a0 = 1.0;
for (int k = 0; k <= 64; k++)
{
double kFactorial = Factorial(k);
coeff[k] = a0 * pow(-1.0, k) /
(pow(2.0, 2.0 * k) * pow(kFactorial, 2.0));
}
}
static void ComputePoints()
{
double chi0 = -15.0;
double deltaChi = 30.0 / 1024;
for (int i = 0; i <= 1024; i++)
{
double sum = 0.0;
for (int k = 64; k >= 0; k--)
sum += pow(chi0, 2 * k) * coeff[k];
chi[i] = chi0;
eta[i] = sum;
chi0 += deltaChi;
}
}
static void FindMinMax(
double x[],
double& min,
double& max)
{
min = x[0];
max = x[0];
for (int i = 1; i <= 1024; i++)
{
if (x[i] < min)
min = x[i];
if (x[i] > max)
max = x[i];
}
}
static void DrawTitles(HDC hdc, RECT clientRect, const std::wstring& fTitle,
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
GetTextExtentPoint32(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
int w = sz.cx;
int h = sz.cy;
TextOut(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());
// x-axis title: "x"
std::wstring xTitle = L"x";
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());
// y-axis title: "p(x)"
std::wstring yTitle = L"J0(x)";
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);
}
//
// FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
// PURPOSE: Processes messages for the main window.
//
// WM_COMMAND - process the application menu
// WM_PAINT - Paint the main window
// WM_DESTROY - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
switch (message)
{
case WM_CREATE:
ComputeCoefficients();
ComputePoints();
break;
case WM_COMMAND:
{
int wmId = LOWORD(wParam);
// Parse the menu selections:
switch (wmId)
{
case IDM_ABOUT:
DialogBox(hInst, MAKEINTRESOURCE(IDD_ABOUTBOX), hWnd, About);
break;
case IDM_EXIT:
DestroyWindow(hWnd);
break;
default:
return DefWindowProc(hWnd, message, wParam, lParam);
}
}
break;
case WM_PAINT:
{
PAINTSTRUCT ps;
HDC hdc = BeginPaint(hWnd, &ps);
ComputePoints();
double xMin = 0, yMin = 0;
double xMax = 0, yMax = 0;
FindMinMax(chi, xMin, xMax);
FindMinMax(eta, yMin, yMax);
double h = 0, pi = 0, plm = 0, theta = 0;
float xSpan = (float)(xMax - xMin);
float ySpan = (float)(yMax - yMin);
RECT rect = { };
GetClientRect(hWnd, &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;
POINT wPt = { 0 };
int i = 0;
float x = (float)xMin;
float y = (float)yMax;
wcscpy_s(fTitle, 128, L"Bessel Function of the First Kind Order 0");
DrawTitles(hdc, rect, fTitle, (int)sx0, (int)sx1, (int)sy0, (int)sy1);
px = (float)chi[0];
py = (float)eta[0];
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 = { 0 };
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 nPenNew = NULL;
HGDIOBJ hPenOld = NULL;
nPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
hPenOld = SelectObject(hdc, nPenNew);
HRGN clipRegion = CreateRectRgn(
(int)sx0, (int)sy0, // Left, Top
(int)(sx1), (int)(sy1) // Right, Bottom
);
// Apply clipping region
SelectClipRgn(hdc, clipRegion);
px = (float)chi[0];
py = (float)eta[0];
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 <= 1024; j++)
{
px = (float)chi[j];
py = (float)eta[j];
sx = xSlope * px + xInter;
sy = ySlope * py + yInter;
LineTo(hdc, (int)sx, (int)sy);
}
DeleteObject(nPenNew);
nPenNew = NULL;
SelectObject(hdc, hPenOld);
DeleteObject(clipRegion);
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;
}
Blog Entry © Thursday and Friday, March 12 – 13, 2026, Naive Power Series Computations of the Cosine and Sine Trigonometric Functions and Polynomial Fit to the Trig Functions Using Least Squares
// TrigLSPolyFit.cpp (c) Thursday, March 12, 2026
// by James Pate Williams, Jr.
// Least Squares Curve Fitting to the Trigonometric
// functions Cosine and Sine
#include <iomanip>
#include <iostream>
#include <vector>
#include "PolyLSFit.h"
static const int HornerN = 50;
double even_factors[HornerN], odd_factors[HornerN + 1];
static double factorial(int n)
{
double nfactorial = 1.0;
for (int i = 2; i <= n; i++)
nfactorial *= i;
return nfactorial;
}
static void computeEvenFactors()
{
int one = 1;
for (int n = 0; n < 50; n += 2)
{
even_factors[n] = one / factorial(n);
one = -one;
}
}
static void computeOddFactors()
{
int one = 1;
for (int n = 1; n < 51; n += 2)
{
odd_factors[n] = one / factorial(n);
one = -one;
}
}
static double fx_cos(double x, int n)
{
double sum = even_factors[n];
for (int i = n - 1; i >= 0; i--)
{
sum = sum * x + even_factors[i];
}
return sum;
}
static double fx_sin(double x, int n)
{
double sum = odd_factors[n];
for (int i = n - 1; i >= 0; i--)
{
sum = sum * x + odd_factors[i];
}
return sum;
}
static const int M = 33, N = 100, N21 = N + N + 1;
std::vector<double> cosX(N21), cosY(N21), cosCoeff;
std::vector<double> sinX(N21), sinY(N21), sinCoeff;
std::vector<double> pX(M), pY(M);
static void LSPolyFit(int HornerN, int degree)
{
double pi = 4.0 * atan(1.0);
double deltaX = 2.0 * pi / N;
double X = -pi;
for (int i = 0; i < N21; i++)
{
cosX[i] = X;
cosY[i] = fx_cos(X, HornerN);
sinX[i] = X;
sinY[i] = fx_sin(X, HornerN + 1);
X += deltaX;
}
PolyLSFit::leastSquares(cosX, cosY, cosCoeff, degree, N21);
PolyLSFit::leastSquares(sinX, sinY, sinCoeff, degree, N21);
std::cout << std::setw(6) << "x";
std::cout << std::setw(17) << "fit cos";
std::cout << std::setw(17) << "C++ cos";
std::cout << std::setw(17) << "% Error";
std::cout << std::endl;
deltaX = 2.0 * pi / (M - 1);
X = -pi;
for (int i = 0; i < M; i++)
{
double Y = cosCoeff[degree];
for (int j = degree - 1; j >= 0; j--)
{
Y = Y * X + cosCoeff[j];
}
pX[i] = X;
pY[i] = Y;
std::cout << std::fixed << std::setprecision(3);
std::cout << std::setw(6);
std::cout << X << ' ';
if (Y >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::fixed << std::setprecision(16);
std::cout << std::setw(16);
std::cout << fabs(Y) << ' ';
double cx = cos(X);
if (cx >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::fixed << std::setprecision(16);
std::cout << std::setw(16);
std::cout << fabs(cx) << ' ';
double pcDifference = 0.0;
if (fabs(cx) >= 1.0e-12)
pcDifference = 100.0 * fabs(Y - cx) / fabs(cx);
else
pcDifference = -1.0;
std::cout << std::fixed << std::setprecision(16);
std::cout << std::setw(16);
if (pcDifference != -1.0)
std::cout << pcDifference << std::endl;
else
std::cout << "not a number" << std::endl;
X += deltaX;
}
std::cout << std::endl;
std::cout << std::setw(6) << "x";
std::cout << std::setw(17) << "fit sin";
std::cout << std::setw(17) << "C++ sin";
std::cout << std::setw(17) << "% Error";
std::cout << std::endl;
deltaX = 2.0 * pi / (M - 1);
X = -pi;
for (int i = 0; i < M; i++)
{
double Y = sinCoeff[degree];
for (int j = degree - 1; j >= 0; j--)
{
Y = Y * X + sinCoeff[j];
}
pX[i] = X;
pY[i] = Y;
std::cout << std::fixed << std::setprecision(3);
std::cout << std::setw(6);
std::cout << X << ' ';
if (Y >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::fixed << std::setprecision(16);
std::cout << std::setw(16);
std::cout << fabs(Y) << ' ';
double sx = sin(X);
if (sx >= 0)
std::cout << '+';
else
std::cout << '-';
std::cout << std::fixed << std::setprecision(16);
std::cout << std::setw(16);
std::cout << fabs(sx) << ' ';
double pcDifference = 0.0;
if (fabs(sx) >= 1.0e-12)
pcDifference = 100.0 * fabs(Y - sx) / fabs(sx);
else
pcDifference = -1.0;
std::cout << std::fixed << std::setprecision(16);
std::cout << std::setw(16);
if (pcDifference != -1.0)
std::cout << pcDifference << std::endl;
else
std::cout << "not a number" << std::endl;
X += deltaX;
}
std::cout << std::endl;
}
int main()
{
computeEvenFactors();
computeOddFactors();
std::cout << "# terms 16" << std::endl;
LSPolyFit(HornerN, 16);
std::cout << "# terms 20" << std::endl;
LSPolyFit(HornerN, 20);
std::cout << "# terms 24" << std::endl;
LSPolyFit(HornerN, 24);
return 0;
}
Numerical Explorations 