Blog Entry © Sunday, April 26, 2026, by James Pate Williams, Jr. Butterfly Curve

Butterfly curve (transcendental) – Wikipedia

// ButterflyEquation.cpp : Defines the entry point for the application.
//

#include "framework.h"
#include "ButterflyEquation.h"
#include <vector>

#define MAX_LOADSTRING 100

typedef struct tagPoint3d
{
    double t, x, y;
} Point3d, * PPoint3d;

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
std::vector<Point3d> 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 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_BUTTERFLYEQUATION, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    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_BUTTERFLYEQUATION));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_BUTTERFLYEQUATION);
    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;
}

void CreateButterflyGraphPoints()
{
    double p = 4.0 * atan(1.0);
    double h = 12.0 * p / 1024.0;
    double t = 0.0;

    for (int i = 0; i <= 1024; i++)
    {
        Point3d pt = { 0 };
        double c = 2.0 * cos(4.0 * t);
        double d = pow(sin(t / 12.0), 5.0);
        double x = sin(t) * (exp(cos(t)) - c - d);
        double y = cos(t) * (exp(cos(t)) - c - d);
        pt.t = t;
        pt.x = x;
        pt.y = y;
        points.push_back(pt);
        t += h;
    }
}

static void FindMinMax(
    double& xMin, double& xMax,
    double& yMin, double& yMax)
{
    // uses global 2D double points structure

    xMin = yMin = DBL_MAX;
    xMax = yMax = DBL_MIN;

    for (size_t i = 0; i < points.size(); i++)
    {
        Point3d pt = points[i];
        double x = pt.x;
        double y = pt.y;

        if (x < xMin)
            xMin = x;
        if (x > xMax)
            xMax = x;
        if (y < yMin)
            yMin = y;
        if (y > yMax)
            yMax = y;
    }
}

static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}

//
//  FUNCTION: WndProc(HWND, UINT, WPARAM, LPARAM)
//
//  PURPOSE: Processes messages for the main window.
//
//  WM_COMMAND  - process the application menu
//  WM_PAINT    - Paint the main window
//  WM_DESTROY  - post a quit message and return
//
//
LRESULT CALLBACK WndProc(HWND hWnd, UINT message, WPARAM wParam, LPARAM lParam)
{
    switch (message)
    {
    case WM_CREATE:
        CreateButterflyGraphPoints();
        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);
            double h = 0, pi = 0, plm = 0, theta = 0;
            double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
            FindMinMax(xMin, xMax, yMin, yMax);
            float xSpan = (float)(xMax - xMin);
            float ySpan = (float)(yMax - yMin);
            RECT rect = { };
            GetClientRect(hWnd, &rect);
            float width = (float)(rect.right - rect.left + 1);
            float height = (float)(rect.bottom - rect.top - 32 + 1);
            float sx0 = 2.0f * width / 16.0f;
            float sx1 = 14.0f * width / 16.0f;
            float sy0 = 2.0f * height / 16.0f;
            float sy1 = 14.0f * height / 16.0f;
            float deltaX = xSpan / 8.0f;
            float deltaY = ySpan / 8.0f;
            float xSlope = (sx1 - sx0) / xSpan;
            float xInter = (float)(sx0 - xSlope * xMin);
            float ySlope = (sy0 - sy1) / ySpan;
            float yInter = (float)(sy0 - ySlope * yMax);
            float px = 0, py = 0, sx = 0, sy = 0;
            POINT wPt = { };
            int i = 0;
            float x = (float)xMin;
            float y = (float)yMax;
            px = x;
            py = y;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
            char buffer[128] = { };

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

                sprintf_s(buffer, "%5.4lf", x);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    buffer,
                    (int)strlen(buffer),
                    &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, buffer, 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)
                {
                    sprintf_s(buffer, "%+5.3lf", y);
                    SIZE size = { };
                    GetTextExtentPoint32A(
                        hdc,
                        buffer,
                        (int)strlen(buffer),
                        &size);
                    RECT textRect = { };
                    textRect.left = (long)(sx0 - size.cx - size.cx / 5.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, buffer, textRect);
                }

                y += deltaY;
                i++;
            }

            HGDIOBJ bPenNew = NULL;
            HGDIOBJ hPenOld = NULL;

            bPenNew = CreatePen(PS_SOLID, 2, RGB(0, 0, 255));
            hPenOld = SelectObject(hdc, bPenNew);

            px = (float)points[0].x;
            py = (float)points[0].y;
            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 < points.size(); j++)
            {
                px = (float)points[j].x;
                py = (float)points[j].y;
                sx = xSlope * px + xInter;
                sy = ySlope * py + yInter;
                LineTo(hdc, (int)sx, (int)sy);
            }

            SelectObject(hdc, hPenOld);
            DeleteObject(bPenNew);
            EndPaint(hWnd, &ps);
        }
        break;
    case WM_DESTROY:
        PostQuitMessage(0);
        break;
    default:
        return DefWindowProc(hWnd, message, wParam, lParam);
    }
    return 0;
}

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

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

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

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

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

// 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

// 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 © Saturday April 4, 2026, by James Pate Williams, Jr., Hydrogen-like Radial Electron Distribution Functions in CPP

#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 Entry © Thursday, April 2, 2026, by James Pate Williams, Jr., Matrix Inverses that Exist and Their Characteristic Polynomials

#pragma once
#include <vector>

class Inverse
{
public:
	static bool Compute(
		int n,
		std::vector<std::vector<double>>& M,
		std::vector<std::vector<double>>& X);
	static void CharacteristicPolynomial(
		int n, std::vector<double>& a,
		std::vector<std::vector<double>>& C,
		std::vector<std::vector<double>>& M,
		std::vector<std::vector<double>>& Madj);
};

// Reference: 5.2: The Characteristic Polynomial - Mathematics LibreTexts
#include "Inverse.h"
#include <vector>

bool Inverse::Compute(
	int n,
	std::vector<std::vector<double>>& M,
	std::vector<std::vector<double>>& X) {
	double d = 0.0;
	int i = 0, j = 0, k = 0, l = 0;
	std::vector<double> C(n);
	std::vector<std::vector<double>> B;
	B.resize(n);
	X.resize(n);
	for (i = 0; i < n; i++)
	{
		B[i].resize(n);
		X[i].resize(n);
		for (j = 0; j < n; j++)
		{
			B[i][j] = 0.0;
			X[i][j] = 0.0;
		}
		B[i][i] = 1.0;
	}
	j = -1;
Step2:
	j++;
	if (j == n)
		goto Step6;
	for (i = j; i < n; i++) {
		if (M[i][j] != 0)
			break;
	}
	if (i == n)
		return false;
	if (i > j) {
		for (l = j; l < n; l++) {
			double temp = M[i][l];
			M[i][l] = M[j][l];
			M[j][l] = temp;
		}
		for (l = 0; l < n; l++) {
			double temp = B[i][l];
			B[i][l] = B[j][l];
			B[j][l] = temp;
		}
	}
	d = 1.0 / M[j][j];
	for (k = j + 1; k < n; k++) {
		C[k] = d * M[k][j];
	}
	for (k = j + 1; k < n; k++) {
		for (l = j + 1; l < n; l++) {
			M[k][l] -= C[k] * M[j][l];
		}
	}
	for (k = j + 1; k < n; k++) {
		for (l = 0; l < n; l++) {
			B[k][l] -= C[k] * B[j][l];
		}
	}
	goto Step2;
Step6:
	for (i = n - 1; i >= 0; i--) {
		std::vector<double> sum(n, 0);
		for (j = i + 1; j < n; j++) {
			for (k = 0; k < n; k++) {
				sum[k] += M[i][j] * X[j][k];
			}
		}
		for (j = 0; j < n; j++)
			X[i][j] = (B[i][j] - sum[j]) / M[i][i];
	}
	return true;
}

double delta(int i, int j) {
	return i == j ? 1.0 : 0.0;
}

void Inverse::CharacteristicPolynomial(
	int n, std::vector<double>& a,
	std::vector<std::vector<double>>& C,
	std::vector<std::vector<double>>& M,
	std::vector<std::vector<double>>& Madj) {
	int i = 0, j = 0, k = 0, l = 0;
	std::vector<double> X(n);
	std::vector<std::vector<double>> MC;
	a.resize(n + 1LL);
	C.resize(n);
	M.resize(n);
	MC.resize(n);
	Madj.resize(n);
	for (i = 0; i < n; i++) {
		C[i].resize(n);
		M[i].resize(n);
		MC[i].resize(n);
		Madj[i].resize(n);
		for (j = 0; j < n; j++) {
			C[i][j] = delta(i, j);
		}
	}
	a[0] = 1.0;
	i = 0;
Step2:
	i++;
	if (i == n) {
		for (j = 0; j < n; j++) {
			for (k = 0; k < n; k++) {
				double sum = 0.0;
				for (l = 0; l < n; l++) {
					sum += M[j][l] * C[l][k];
				}
				MC[j][k] = sum;
			}
		}
		double trace = 0.0;
		for (j = 0; j < n; j++) {
			trace += MC[j][j];
		}
		a[n] = -trace / n;
		for (j = 0; j < n; j++) {
			for (k = 0; k < n; k++) {
				Madj[j][k] = pow(-1, n - 1) * C[j][k];
			}
		}
		return;
	}
	for (j = 0; j < n; j++) {
		for (k = 0; k < n; k++) {
			double sum = 0.0;
			for (l = 0; l < n; l++) {
				sum += M[j][l] * C[l][k];
			}
			MC[j][k] = sum;
		}
	}
	for (j = 0; j < n; j++) {
		for (k = 0; k < n; k++) {
			C[j][k] = MC[j][k];
		}
	}
	double tr = 0.0;
	for (j = 0; j < n; j++) {
		tr += C[j][j];
	}
	a[i] = -tr / i;
	for (j = 0; j < n; j++) {
		for (k = 0; k < n; k++) {
			C[j][k] += a[i] * delta(j, k);
		}
	}
	goto Step2;
}