Category: Elementary Physics

Vector Analysis by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition© 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Pages 48 and 50

Vector Analysis Chapter 1 by James Pate WilliamsDownload
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider 

#include <iostream>
#include <vector>

static double InnerProduct(
	std::vector<double> A,
	std::vector<double> B,
	int n)
{
	double sum = 0.0;

	for (int i = 0; i < n; i++)
		sum += A[i] * B[i];

	return sum;
}

static void VectorProduct(
	std::vector<double> A,
	std::vector<double> B,
	std::vector<double>&C)
{
	C.resize(3);
	C[0] = A[1] * B[2] - A[2] * B[1];
	C[1] = A[2] * B[0] - A[0] * B[2];
	C[2] = A[0] * B[1] - A[1] * B[0];
}

static double TripleProduct(
	std::vector<double> A,
	std::vector<double> B,
	std::vector<double> C)
{
	double sum0 = 0.0, sum1 = 0.0;

	sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
	sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
	return sum0 - sum1;
}

static void Exercises_Section_1_13_1_Triple_Products()
{
	// triple products (a)
	std::vector<double> A1 = { 2, 0, 0 };
	std::vector<double> B1 = { 0, 3, 0 };
	std::vector<double> C1 = { 0, 0, 5 };
	double tp_1 = TripleProduct(A1, B1, C1);
	std::cout << "Exercise 1 (a)" << '\t';
	std::cout << tp_1 << std::endl;
	// triple products (b)
	std::vector<double> A2 = { 1, 1, 1 };
	std::vector<double> B2 = { 3, 1, 0 };
	std::vector<double> C2 = { 0, -1, 5 };
	std::cout << "Exercise 1 (b)" << '\t';
	double tp_2 = TripleProduct(A2, B2, C2);
	std::cout << tp_2 << std::endl;
	// triple products (c)
	std::vector<double> A3 = { 2, -1, 1 };
	std::vector<double> B3 = { 1, 1, 1 };
	std::vector<double> C3 = { 2, 0, 3 };
	double tp_3 = TripleProduct(A3, B3, C3);
	std::cout << "Exercise 1 (c)" << '\t';
	std::cout << tp_3 << std::endl;
	// triple products (d)
	std::vector<double> A4 = { 0, 0, 1 };
	std::vector<double> B4 = { 1, 0, 0 };
	std::vector<double> C4 = { 0, 1, 0 };
	double tp_4 = TripleProduct(A4, B4, C4);
	std::cout << "Exercise 1 (d)" << '\t';
	std::cout << tp_4 << std::endl;
	// volume of a parallelpipped
	std::vector<double> A5 = { 3, 4, 1 };
	std::vector<double> B5 = { 2, 3, 4 };
	std::vector<double> C5 = { 0, 0, 5 }; 
	double tp_5 = TripleProduct(A5, B5, C5);
	std::cout << "Exercise 2" << '\t';
	std::cout << tp_5 << std::endl;
	// volume of a parallelpipped
	std::vector<double> A6 = { 3, 2, 1 };
	std::vector<double> B6 = { 4, 2, 1 };
	std::vector<double> C6 = { 0, 1, 4 };
	std::vector<double> D6 = { 0, 0, 7 };
	std::vector<double> AB = { -1, 0, 0 };
	std::vector<double> AC = { 3, 1, -3 };
	std::vector<double> AD = { 3, 2, -6 };
	std::cout << "Exercise 3" << '\t';
	double tp_6 = TripleProduct(AB, AC, AD);
	std::cout << tp_6 << std::endl;
	// volume of a tetrahedron
	std::vector<double> AB1 = { 1, 1, 0 };
	std::vector<double> AC1 = { 1, -1, 0 };
	std::vector<double> AD1 = { 0, 0, 2 };
	std::cout << "Exercise 4" << '\t';
	double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
	std::cout << fabs(tp_7) << std::endl;
	std::vector<double> P1 = { 0, 0, 0 };
	std::vector<double> P2 = { 1, 1, 0 };
	std::vector<double> P3 = { 3, 4, 0 };
	std::vector<double> P4 = { 4, 5, 0 };
	std::vector<double> P5 = { 0, 0, 1 };
	std::vector<double> Q1 = { 1, 1, 0 };
	std::vector<double> Q2 = { 2, 3, 0 };
	std::vector<double> Q3 = { 1, 1, 1 };
	std::cout << "Exercise 5" << '\t';
	double tp_8 = TripleProduct(Q1, Q2, Q3);
	std::cout << fabs(tp_8) << std::endl;
	std::vector<double> A10 = { 1, 1, 1 };
	std::vector<double> B10 = { 2, 4, -1 };
	std::vector<double> C10 = { 1, 1, 3 };
	double tp_10 = TripleProduct(A10, B10, C10);
	std::vector<double> D10 = { 0 };
	VectorProduct(A10, B10, D10);
	double magnitude = sqrt(InnerProduct(D10, D10, 3));
	std::cout << "Exercise 10" << '\t';
	std::cout << tp_10 / magnitude << '\t';
	std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
	std::vector<double> A11 = { 1, 1, 1 };
	std::vector<double> B11 = { 2, 4, -1 };
	std::vector<double> C11 = { 1, 1, 3 };
	std::vector<double> D11 = { 3, 2, 1 };
	std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
	VectorProduct(A11, B11, AB11);
	VectorProduct(B11, C11, BC11);
	VectorProduct(C11, A11, CA11);
	VectorProduct(BC11, CA11, BCCA11);
	double Q11 = InnerProduct(AB11, BCCA11, 3);
	double A_x = A11[0], A_y = A11[1], A_z = A11[2];
	double B_x = B11[0], B_y = B11[1], B_z = B11[2];
	double C_x = C11[0], C_y = C11[1], C_z = C11[2];
	double term1 = +(A_y * B_z - A_z * B_y) * (B_z * C_x - B_x * C_z) * (C_x * A_y - C_y * A_x);
	double term2 = -(A_y * B_z - A_z * B_y) * (B_x * C_y - B_y * C_x) * (C_z * A_x - C_x * A_z);
	double term3 = +(A_z * B_x - A_x * B_z) * (B_x * C_y - B_y * C_x) * (C_y * A_z - C_z * A_y);
	double term4 = -(A_z * B_x - A_x * B_z) * (B_y * C_z - B_z * C_y) * (C_x * A_y - C_y * A_x);
	double term5 = +(A_x * B_y - A_y * B_x) * (B_y * C_z - B_z * C_y) * (C_z * A_x - C_x * A_z);
	double term6 = -(A_x * B_y - A_y * B_x) * (B_z * C_x - B_x * C_z) * (C_y * A_z - C_z * A_y);
	double P11 = term1 + term2 + term3 + term4 + term5 + term6;
	std::cout << "Q = (A x B) . (B x C) x (C x A) = " << Q11 << std::endl;
	std::cout << "P = (A x B) . (B x C) x (C x A) = " << P11 << std::endl;
}

static void Exercises_Section_1_14_Vector_Identities()
{
	std::vector<double> A11 = { 1, 1, 1 };
	std::vector<double> B11 = { 2, 4, -1 };
	std::vector<double> C11 = { 1, 1, 3 };
	std::vector<double> D11 = { 3, 2, 1 };
	std::cout << "Section 1.14 page 50 Exercises Exercise 1" << std::endl;
	std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
	std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
	std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
	std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
	std::cout << "TPI1 = (A x B) x (C x D) = [A, B, D]C - [A, B, C]D = " << std::endl;
	double TP1411a = TripleProduct(A11, B11, D11);
	double TP1411b = TripleProduct(A11, B11, C11);
	std::cout << "[A, B, D] = " << TP1411a << std::endl;
	std::cout << "[A, B, C] = " << TP1411b << std::endl;
	std::cout << "TPI1_x = " << TP1411a * C11[0] << std::endl;
	std::cout << "TPI1_y = " << TP1411a * C11[1] << std::endl;
	std::cout << "TPI1_z = " << TP1411a * C11[2] << std::endl;
	std::cout << "TPI2_x = " << TP1411b * D11[0] << std::endl;
	std::cout << "TPI2_y = " << TP1411b * D11[1] << std::endl;
	std::cout << "TPI2_z = " << TP1411b * D11[2] << std::endl;
	std::cout << "RHS1 = [A, B, D]C - [A, B, C]D = " << std::endl;
	std::vector<double> RHS1(3);
	RHS1[0] = TP1411a * C11[0] - TP1411b * D11[0];
	RHS1[1] = TP1411a * C11[1] - TP1411b * D11[1];
	RHS1[2] = TP1411a * C11[2] - TP1411b * D11[2];
	std::cout << "RHS1_x = " << RHS1[0] << std::endl;
	std::cout << "RHS1_y = " << RHS1[1] << std::endl;
	std::cout << "RHS1_z = " << RHS1[2] << std::endl;
	std::vector<double> CD11(3), TD11(3);
	std::vector<double> AB11(3), BC11(3), CA11(3), BCCA11(3);
	VectorProduct(A11, B11, AB11);
	VectorProduct(B11, C11, BC11);
	VectorProduct(C11, A11, CA11);
	VectorProduct(BC11, CA11, BCCA11);
	VectorProduct(A11, B11, AB11);
	VectorProduct(C11, D11, CD11);
	VectorProduct(AB11, CD11, TD11);
	std::cout << "A = " << A11[0] << '\t' << A11[1] << '\t' << A11[2] << std::endl;
	std::cout << "B = " << B11[0] << '\t' << B11[1] << '\t' << B11[2] << std::endl;
	std::cout << "C = " << C11[0] << '\t' << C11[1] << '\t' << C11[2] << std::endl;
	std::cout << "D = " << D11[0] << '\t' << D11[1] << '\t' << D11[2] << std::endl;
	std::cout << "A x B = " << AB11[0] << '\t' << AB11[1] << '\t' << AB11[2] << std::endl;
	std::cout << "C x D = " << CD11[0] << '\t' << CD11[1] << '\t' << CD11[2] << std::endl;
	std::cout << "TD11 = (A x B) x (C x D) = " << std::endl;
	std::cout << "TD11_x = " << TD11[0] << std::endl;
	std::cout << "TD11_y = " << TD11[1] << std::endl;
	std::cout << "TD11_z = " << TD11[2] << std::endl;
	VectorProduct(B11, C11, BC11);
	VectorProduct(C11, A11, CA11);
	VectorProduct(BC11, CA11, D11);
	VectorProduct(A11, B11, AB11);
	double ip12 = InnerProduct(AB11, D11, 3);
	std::cout << "2. Inner Product = " << ip12 << std::endl;
	double tp12 = TripleProduct(A11, B11, C11);
	std::cout << "2. Triple Product ^ 2 = " << tp12 * tp12 << std::endl;
	std::vector<double> ABC11(3), BAC11(3), CAB11(3);
	VectorProduct(A11, BC11, ABC11);
	VectorProduct(B11, CA11, BAC11);
	VectorProduct(C11, AB11, CAB11);
	double zx = ABC11[0] + BAC11[0] + CAB11[0];
	double zy = ABC11[1] + BAC11[1] + CAB11[1];
	double zz = ABC11[2] + BAC11[2] + CAB11[2];
	std::cout << "3. Zero Vector = " << zx << ' ' << zy << ' ' << zz;
	std::cout << std::endl;
}

int main()
{
	Exercises_Section_1_13_1_Triple_Products();
	Exercises_Section_1_14_Vector_Identities();
	return 0;
}

Blog Entry © Tuesday, April 14, 2026, by James Pate Williams, Jr. Exercises and Supplementary Problems from Introduction to Vector Analysis Fourth Edition © 1979 by Harry F. Davis and Arthur David Snider Selected Exercises from Chapter 1 Page 48

Vector Analysis Chapter 1 by James Pate WilliamsDownload
// VectorAnalysis.cpp © Tuesday, April 14, 2026
// by James Pate Williams, Jr.
// Reference: Introduction to Vector Analysis Fourth Edition
// © 1979 by Harry F. Davis and Arthur David Snider 

#include <iostream>
#include <vector>

static double InnerProduct(
	std::vector<double> A,
	std::vector<double> B,
	int n)
{
	double sum = 0.0;

	for (int i = 0; i < n; i++)
		sum += A[i] * B[i];

	return sum;
}

static void VectorProduct(
	std::vector<double> A,
	std::vector<double> B,
	std::vector<double>&C)
{
	C.resize(3);
	C[0] = A[1] * B[2] - A[2] * B[1];
	C[1] = A[0] * B[2] - A[2] * B[0];
	C[2] = A[0] * B[1] - A[1] * B[0];
}

static double TripleProduct(
	std::vector<double> A,
	std::vector<double> B,
	std::vector<double> C)
{
	double sum0 = 0.0, sum1 = 0.0;

	sum0 += A[0] * B[1] * C[2] + A[1] * B[2] * C[0] + A[2] * B[0] * C[1];
	sum1 += A[1] * B[0] * C[2] + A[0] * B[2] * C[1] + A[2] * B[1] * C[0];
	return sum0 - sum1;
}

static void Exercises_Section_1_13_1_Triple_Products()
{
	// triple products (a)
	std::vector<double> A1 = { 2, 0, 0 };
	std::vector<double> B1 = { 0, 3, 0 };
	std::vector<double> C1 = { 0, 0, 5 };
	double tp_1 = TripleProduct(A1, B1, C1);
	std::cout << "Exercise 1 (a)" << '\t';
	std::cout << tp_1 << std::endl;
	// triple products (b)
	std::vector<double> A2 = { 1, 1, 1 };
	std::vector<double> B2 = { 3, 1, 0 };
	std::vector<double> C2 = { 0, -1, 5 };
	std::cout << "Exercise 1 (b)" << '\t';
	double tp_2 = TripleProduct(A2, B2, C2);
	std::cout << tp_2 << std::endl;
	// triple products (c)
	std::vector<double> A3 = { 2, -1, 1 };
	std::vector<double> B3 = { 1, 1, 1 };
	std::vector<double> C3 = { 2, 0, 3 };
	double tp_3 = TripleProduct(A3, B3, C3);
	std::cout << "Exercise 1 (c)" << '\t';
	std::cout << tp_3 << std::endl;
	// triple products (d)
	std::vector<double> A4 = { 0, 0, 1 };
	std::vector<double> B4 = { 1, 0, 0 };
	std::vector<double> C4 = { 0, 1, 0 };
	double tp_4 = TripleProduct(A4, B4, C4);
	std::cout << "Exercise 1 (d)" << '\t';
	std::cout << tp_4 << std::endl;
	// volume of a parallelpipped
	std::vector<double> A5 = { 3, 4, 1 };
	std::vector<double> B5 = { 2, 3, 4 };
	std::vector<double> C5 = { 0, 0, 5 }; 
	double tp_5 = TripleProduct(A5, B5, C5);
	std::cout << "Exercise 2" << '\t';
	std::cout << tp_5 << std::endl;
	// volume of a parallelpipped
	std::vector<double> A6 = { 3, 2, 1 };
	std::vector<double> B6 = { 4, 2, 1 };
	std::vector<double> C6 = { 0, 1, 4 };
	std::vector<double> D6 = { 0, 0, 7 };
	std::vector<double> AB = { -1, 0, 0 };
	std::vector<double> AC = { 3, 1, -3 };
	std::vector<double> AD = { 3, 2, -6 };
	std::cout << "Exercise 3" << '\t';
	double tp_6 = TripleProduct(AB, AC, AD);
	std::cout << tp_6 << std::endl;
	// volume of a tetrahedron
	std::vector<double> AB1 = { 1, 1, 0 };
	std::vector<double> AC1 = { 1, -1, 0 };
	std::vector<double> AD1 = { 0, 0, 2 };
	std::cout << "Exercise 4" << '\t';
	double tp_7 = TripleProduct(AB1, AC1, AD1) / 6.0;
	std::cout << fabs(tp_7) << std::endl;
	std::vector<double> P1 = { 0, 0, 0 };
	std::vector<double> P2 = { 1, 1, 0 };
	std::vector<double> P3 = { 3, 4, 0 };
	std::vector<double> P4 = { 4, 5, 0 };
	std::vector<double> P5 = { 0, 0, 1 };
	std::vector<double> Q1 = { 1, 1, 0 };
	std::vector<double> Q2 = { 2, 3, 0 };
	std::vector<double> Q3 = { 1, 1, 1 };
	std::cout << "Exercise 5" << '\t';
	double tp_8 = TripleProduct(Q1, Q2, Q3);
	std::cout << fabs(tp_8) << std::endl;
	std::vector<double> A10 = { 1, 1, 1 };
	std::vector<double> B10 = { 2, 4, -1 };
	std::vector<double> C10 = { 1, 1, 3 };
	double tp_10 = TripleProduct(A10, B10, C10);
	std::vector<double> D10 = { 0 };
	VectorProduct(A10, B10, D10);
	double magnitude = sqrt(InnerProduct(D10, D10, 3));
	std::cout << "Exercise 10" << '\t';
	std::cout << tp_10 / magnitude << '\t';
	std::cout << 2.0 * sqrt(38.0) / 19.0 << std::endl;
}

int main()
{
	Exercises_Section_1_13_1_Triple_Products();
	return 0;
}

Blog Entry © Tuesday, April 7, 2026, by James Pate Williams, Jr., Hydrogen-like Atom Polar and Azimuthal Wavefunctions

Blog Entry Tuesday April 7 2026Download

Blog Entry © Saturday April 4, 2026, by James Pate Williams, Jr., Hydrogen-like Radial Electron Distribution Functions in CPP

Blog Entry Saturday April 4 2026Download
#pragma once

class RadialWaveFunction
{
private:
	static double Factorial(int n);
	static double Laguerre(double rho, int n, int l);
public:
	static double R(double r, int Z, int n, int l);
};


#include <math.h>
#include "RadialWaveFunction.h"

double RadialWaveFunction::Factorial(int n)
{
	double factorial = 1.0;

	for (int i = 2; i <= n; i++)
		factorial *= i;

	return factorial;
}

double RadialWaveFunction::Laguerre(double rho, int n, int l)
{
	double sum = 0.0;

	for (int k = 0; k <= n - l - 1; k++)
	{
		double factor1 = pow(-1, k + 2 * l + 1);
		double factor2 = pow(Factorial(n + l), 2.0);
		double factor3 = pow(rho, k);
		double numer = factor1 * factor2 * factor3;
		double factor4 = Factorial(n - l - 1 - k);
		double factor5 = Factorial(2 * l + 1 + k);
		double factor6 = Factorial(k);
		double denom = factor4 * factor5 * factor6;

		sum += numer / denom;
	}

	return sum;
}

double RadialWaveFunction::R(double r, int Z, int n, int l)
{
	double rho = 2.0 * Z * r / n;
	double numer1 = pow(2.0 * Z / n, 3.0);
	double numer2 = Factorial(n - l - 1);
	double denom1 = Factorial(n + n);
	double denom2 = pow(Factorial(n + l), 3.0);
	double numer3 = -sqrt(numer1 * numer2 /
		(denom1 * denom2));
	double exp2 = exp(-0.5 * rho);
	double rhol = pow(rho, l);
	return numer3 * exp2 * rhol * Laguerre(rho, n, l);
}

Blog Input © Saturday, March 21, 2026, by James Pate Williams, Jr., Three Model Universes

Blog Entry Saturday March 21 2026Download
#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 5, 2026, by James Pate Williams, Jr., Simple Filtering of a Noisy Signal in the Programming Language C

Blog Entry Thursday March 5 2026Download

Blog Entry (c) Tuesday February 3, 2026, by James Pate Williams, Jr. Derivation of the Classical Kinetic Energy Formula

Derivation of the Classical Kinetic Energy Equation from EinsteinDownload
// KineticEnergy.cpp
// Author: James Pate Williams, Jr. and NIST
// Copyright Monday February 2, 2026

#include <iomanip>
#include <iostream>
#include <vector>

const double c = 2.99792458E8;
const double mass0 = 1.67492750056E-27;

static double KineticEnergy(
    double v, size_t number,
    double& kilotonsTNT,
    std::vector<double>& terms)
{
    double c2 = c * c, ke = 0;;
    
    if (number >= 0)
        terms[0] = mass0 * c2;
    if (number >= 1)
        terms[1] = 0.5 * mass0 * v * v;
    if (number >= 2)
        terms[2] = 3.0 * mass0 * pow(v, 4.0) / (8.0 * c2);
    if (number >= 3)
        terms[3] = 225.0 * mass0 * pow(v, 6.0) / (720.0 * c2 * c2);

    for (size_t i = 0; i <= number; i++)
        ke += terms[i];

    kilotonsTNT = 2.3900573613767E-13 * ke;
    return ke;
}

int main()
{
    char line[128] = "";

    while (true)
    {
        double ke = 0.0, kilotonsTNT = 0.0, v = 0.0;
        size_t number = 0;

        std::cout << "Enter v = ";
        std::cin.getline(line, 128);
        v = atof(line);

        if (v == 0)
            break;

        std::cout << "# Terms = ";
        std::cin.getline(line, 128);
        number = static_cast<size_t>(atoi(line));
        std::vector<double> terms(number + 1);
        ke = KineticEnergy(v, number, kilotonsTNT, terms);
        std::cout << "KE = " << ke << " joules" << std::endl;
        std::cout << "KE = " << kilotonsTNT << " kilotons TNT";
        std::cout << std::endl;
        std::cout << std::scientific;

        for (size_t i = 0; i <= number; i++)
        {
            std::cout << "KE[" << i << "] = " << terms[i];
            std::cout << " joules" << std::endl;
            kilotonsTNT = 2.3900573613767E-13 * terms[i];
            std::cout << "KE[" << i << "] = " << kilotonsTNT;
            std::cout << " kilotons TNT" << std::endl;
        }
    }
}

Blog Entry © Sunday, January 25, 2026, by James Pate Williams, Jr., Schwarzschild Solution of Einstein’s General Relativity Gravitational Field Equation

Blog Entry Sunday January 25 2026Download

Blog Entry © Friday, January 16, 2026, by James Pate Williams, Jr., Another Update of My Iowa Class Battleship Artillery Exterior Ballistics Application

Blog Entry Friday January 16 2026Download

Blog Entry © Wednesday, January 14, 2026, by James Pate Williams, Jr. Curvature of the Earth Table

Blog Entry Wednesday January 14 2026Download
// CurvatureOfTheEarth.cpp : Defines the entry point for the application.
//

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

#define MAX_LOADSTRING 100

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

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

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

    // TODO: Place code here.

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

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

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

    MSG msg;

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

    return (int) msg.wParam;
}

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

    wcex.cbSize = sizeof(WNDCLASSEX);

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

    return RegisterClassExW(&wcex);
}

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

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

   if (!hWnd)
   {
      return FALSE;
   }

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

   return TRUE;
}

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

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

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

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

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

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

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

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

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

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

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

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

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

    return (INT_PTR)FALSE;
}


#pragma once
#include "Vincenty.h"

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

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

class GreatCircleDistance
{

private:

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

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

public:

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

        return 0;
    };

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

    double Interpolate(double R);

    // construction of the curvature of the Earth table

    GreatCircleDistance();
};


#pragma once

class Vincenty
{

public:

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


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

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

Label10:

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

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

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

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

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

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

    lt = rt = -1;
    return false;
}

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

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

        return ra;
    }

    return -1.0;
}

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

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

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

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

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

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


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

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

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

    return atan2(numer, denom);
}

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

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

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

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

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

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

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

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