Category: Elementary Physics

Blog Entry © Thursday, May 28, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Solution of the Potential Equation in Rectangle using Fixed Point Iteration in Python

# NOTE:
# This implementation prioritizes clarity and correctness over optimization.
# Further performance improvements can be made if needed.
# (c) May 26, 2026 by James Pate Willims, Jr.
# I had some help from the Microsoft Copilot
# to calculate runtimes and define matrices
# Computes the potential in a rectangle
# Reference: "Boundary Value Problems
# Second Edition" by David L. Powers
# See pages 179 to 182 for the analytic
# solution of this Laplace Equation
# Stand alone application using
# Microsoft Visual Studio 2022
# Community Version

import math
import time

xi = yi = 10
u = [[0.0 for _ in range(xi + 2)] for _ in range(yi + 2)]
v = [[0.0 for _ in range(xi + 2)] for _ in range(yi + 2)]

def ComputeBoundaryValues(x, y):
    if x == 0:
        return 0
    if x == 1:
        return 0
    if y == 0 or y == 1:
        if x > 0.0 and x < 0.5:
            return 2.0 * x
        elif x >= 0.5 and x < 1.0:
            return 2.0 - 2.0 * x
                 
    return 0.0

def ComputeParams(its, norm, params):
    params['iterations'] = its
    params['norm'] = norm

def Compute(h, k, xi, yi, maxIts, params):
    # Use a simple fixed-point iteration to
    # compute an approximate solution
    for i in range(0, xi + 1):
        for j in range(0, yi + 1):
            u[i][j] = ComputeBoundaryValues(i * h, j * k)

    for its in range(1, maxIts + 1):
        for i in range(1, xi):
            for j in range(1, yi): 
                u[i][j] = 0.25 * (u[i + 1][j] + u[i - 1][j] + u[i][j + 1] + u[i][j - 1]);

    norm = 0

    for i in range(0, xi + 1):
        for j in range(0, yi + 1):
            norm += math.fabs(u[i][j] * u[i][j])

    norm = math.sqrt(norm)
    params['iterations'] = its
    params['norm'] = norm

def f(x, y):
    # Analytic solution series expansion n = 1 to 100 
    sum = 0.0

    for n in range(1, 101):
        factor1 = math.sin(n * math.pi / 2.0) / (n * n)
        factor2 = math.sinh(n * math.pi * y)
        factor3 = math.sinh(n * math.pi * (1 - y))
        factor4 = math.sin(n * math.pi * x)
        term = (factor2 + factor3) / math.sinh(n * math.pi)
        sum += factor1 * term * factor4
    return 8.0 * sum / (math.pi * math.pi)

avgPE = 0
deltaX = 1.0 / xi
deltaY = 1.0 / yi
maxIts = 50
start_time = time.perf_counter()

for i in range(0, xi + 1):
    for j in range(0, yi + 1):
        v[i][j] = f(i * deltaX, j * deltaY)

minPE = +1000000000
maxPE = -1000000000
params = {}

Compute(deltaX, deltaY, xi, yi, maxIts, params)
print("Approximate\tAnalytic\tPercent Error")

for i in range(0, xi + 1):
    for j in range(0, yi + 1):
        if (math.fabs(u[i][j]) > 1.0e-12 and
            math.fabs(v[i][j]) > 1.0e-12):
            pe = 100.0 * math.fabs((v[i][j] - u[i][j]) / v[i][j])
        else:
            pe = 0.0

        avgPE += pe

        if (pe < minPE):
            minPE = pe

        if (pe > maxPE):
            maxPE = pe

        if math.fabs(pe) != 0.0:
            print("{:10.8f}".format(u[i][j]), "\t", "{:10.8f}".format(v[i][j]), "\t", "{:10.8f}".format(pe))

avgPE /= (xi * yi)
end_time = time.perf_counter()
# Calculate elapsed time in milliseconds
elapsed_ms = (end_time - start_time) * 1000

print("Iterations = ", params['iterations'])
print("Frobenius Norm = ", params['norm'])
print("Minimum Percent Error = ", "{:10.8f}".format(minPE))
print("Average Percent Error = ", "{:10.8f}".format(avgPE))
print("Maximum Percent Error = ", "{:10.8f}".format(maxPE))
print("Elapsed Milliseconds  = ", "{:10.8f}".format(elapsed_ms))


Approximate     Analytic        Percent Error
0.20000000       0.19999972      0.00013831
0.16633455       0.16663592      0.18085704
0.13739427       0.13768928      0.21425509
0.11591159       0.11605132      0.12040154
0.10292732       0.10291871      0.00836375
0.09864668       0.09854114      0.10710198
0.10305612       0.10291871      0.13351975
0.11613103       0.11605132      0.06868445
0.13763395       0.13768928      0.04018144
0.16650309       0.16663592      0.07971563
0.20000000       0.19999972      0.00013831
0.40000000       0.39999927      0.00018169
0.32813882       0.32854798      0.12453472
0.26763492       0.26776105      0.04710489
0.22369849       0.22344522      0.11334607
0.19755216       0.19705465      0.25246951
0.18899096       0.18834090      0.34515051
0.19778524       0.19705465      0.37075402
0.22409557       0.22344522      0.29105301
0.26806863       0.26776105      0.11486989
0.32844379       0.32854798      0.03171204
0.40000000       0.39999927      0.00018169
0.60000000       0.59999811      0.00031532
0.47888974       0.47875999      0.02710181
0.38176991       0.38059768      0.30799710
0.31425594       0.31267225      0.50650231
0.27518847       0.27354681      0.60013876
0.26255192       0.26082096      0.66365837
0.27549367       0.27354681      0.71170992
0.31477587       0.31267225      0.67278721
0.38233779       0.38059768      0.45720485
0.47928905       0.47875999      0.11050679
0.60000000       0.59999811      0.00031532
0.80000000       0.79999202      0.00099701
0.60602379       0.60222488      0.63081180
0.46685956       0.46199176      1.05365615
0.37704317       0.37308607      1.06063957
0.32711029       0.32392135      0.98448025
0.31121899       0.30818168      0.98555933
0.32745161       0.32392135      1.08985035
0.37762462       0.37308607      1.21648870
0.46749464       0.46199176      1.19112076
0.60647034       0.60222488      0.70496220
0.80000000       0.79999202      0.00099701
1.00000000       0.99594729      0.40692036
0.67874673       0.65811281      3.13531720
0.50319768       0.49282441      2.10486107
0.40066316       0.39470092      1.51057096
0.34574699       0.34157202      1.22228048
0.32848272       0.32468552      1.16950228
0.34608839       0.34157202      1.32223096
0.40124475       0.39470092      1.65792212
0.50383291       0.49282441      2.23375660
0.67919339       0.65811281      3.20318626
1.00000000       0.99594729      0.40692036
0.80000000       0.79999202      0.00099701
0.60615279       0.60222488      0.65223206
0.46709299       0.46199176      1.10418350
0.37734885       0.37308607      1.14257161
0.32745218       0.32392135      1.09002599
0.31156100       0.30818168      1.09653477
0.32776105       0.32392135      1.18538026
0.37787503       0.37308607      1.28360546
0.46766769       0.46199176      1.22857864
0.60655688       0.60222488      0.71933113
0.80000000       0.79999202      0.00099701
0.60000000       0.59999811      0.00031532
0.47910949       0.47875999      0.07300025
0.38216756       0.38059768      0.41247566
0.31477665       0.31267225      0.67303766
0.27577086       0.27354681      0.81304189
0.26313452       0.26082096      0.88702848
0.27602079       0.27354681      0.90440981
0.31520243       0.31267225      0.80920943
0.38263258       0.38059768      0.53465917
0.47943646       0.47875999      0.14129609
0.60000000       0.59999811      0.00031532
0.40000000       0.39999927      0.00018169
0.32837881       0.32854798      0.05149022
0.26806920       0.26776105      0.11508237
0.22426717       0.22344522      0.36785094
0.19818820       0.19705465      0.57524397
0.18962723       0.18834090      0.68297865
0.19836093       0.19705465      0.66290001
0.22456142       0.22344522      0.49953914
0.26839057       0.26776105      0.23510691
0.32860478       0.32854798      0.01728752
0.40000000       0.39999927      0.00018169
0.20000000       0.19999972      0.00013831
0.16650327       0.16663592      0.07960469
0.13769959       0.13768928      0.00748883
0.11631140       0.11605132      0.22411146
0.10337449       0.10291871      0.44285504
0.09909401       0.09854114      0.56105764
0.10346087       0.10291871      0.52678322
0.11645855       0.11605132      0.35090580
0.13786030       0.13768928      0.12420922
0.16661627       0.16663592      0.01179315
0.20000000       0.19999972      0.00013831
Iterations =  50
Frobenius Norm =  4.028216200275417
Minimum Percent Error =  0.00000000
Average Percent Error =  0.54286140
Maximum Percent Error =  3.20318626
Elapsed Milliseconds  =  36.23520000
Press any key to continue . . .

Blog Entry (c) Tuesday, May 26, 2026, by James Pate Williams, Jr. and Microsoft’s Copilot Hydrogen-Like Atomic Radial Wave Functions

Included a downloadable PDF and Microsoft Excel Workbook.

Hydrogen Like Radial Wave FunctionsDownload
Hydrogen Like Atomic Wave FunctonsDownload

Blog Entry (c) Monday, May 25, 2026, by James Pate Williams, Jr. Quantum Mechanical Linear Harmonic Oscillator

Python Quantum Mechanical Linear Harmonic Oscillator Wave Functions n0 to n10Download
Python Quantum Mechanical Linear Harmonic Oscillator WorkbookDownload

Blog Entry © Monday, May 11, 2026, by James Pate Williams, Jr., Laplace Equation in a Solid Cylinder

Blog Entry Monday May 11 2026Download
Laplace Equation in a Solid CylinderDownload

Blog Entry © Monday, April 20, 2026, by James Pate Williams, Jr., Vector Analysis Continued and Perhaps Corrected

Blog Entry Monday April 20 2026Download

Blog Entry © Sunday, April 19, 2026, by James Pate Williams, Jr., Scattering from a Spherically Symmetric Potential

Blog Entry Sunday April 19 2026Download

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);
}