Blog Entry © Thursday, October 30, 2025, by James Pate Williams, Jr. Quantum Mechanical Harmonic Oscillator See Introduction to Quantum Mechanics with Applications to Chemistry by Linus Pauling and E. Bright Wilson, Jr. Chapter III

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

#include "pch.h"
#include "framework.h"
#include "QMHarmonicOscillator.h"
#include "Integration.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 buffer[128];
RECT rect;
bool draw = false;
int n;
std::wstring fTitle, yTitle;
std::vector<double> xi, fx;

// 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_QMHARMONICOSCILLATOR, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    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 = { };

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

#define IDC_COMPUTE_BUTTON  1010
#define IDC_CANCEL_BUTTON   1020
#define IDC_CLEAR_BUTTON    1030

static double Factorial(int n)
{
    // n! = n ... 2 1

    double factorial = 1.0;

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

    return factorial;
}

static double N(int n)
{
    // normilization factor

    double pi = 4.0 * atan(1.0);
    double factor1 = sqrt(1.0 / pi);
    double factor2 = 1.0 / (pow(2.0, n) * Factorial(n));
    return sqrt(factor1 * factor2);
}

static double H(double xi, int n)
{
    // Hermite polynomial

    if (n == 1)
    {
        return 2.0 * xi;
    }

    else if (n == 2)
    {
        return 4.0 * xi * xi - 2.0;
    }

    else if (n == 3)
    {
        return xi * (8.0 * xi * xi - 12.0);
    }

    else if (n == 4)
    {
        return 16.0 * pow(xi, 4.0) - 48.0 * pow(xi, 2.0) + 12.0;
    }

    else if (n == 5)
    {
        return 32.0 * pow(xi, 5.0) - 160.0 * pow(xi, 3.0) +
            120.0 * xi;
    }

    else if (n == 6)
    {
        return 64.0 * pow(xi, 6.0) - 480.0 * pow(xi, 4.0) +
            720.0 * xi * xi - 120.0;
    }

    else if (n == 7)
    {
        return 128.0 * pow(xi, 7.0) - 1344.0 * pow(xi, 5.0) +
            3360.0 * xi * xi * xi - 1680.0 * xi;
    }

    else if (n == 8)
    {
        return 256.0 * pow(xi, 8.0) - 3584.0 * pow(xi, 6.0) +
            13440.0 * pow(xi, 4.0) - 13440.0 * xi * xi + 1680.0;
    }

    else if (n == 9)
    {
        return 512.0 * pow(xi, 9.0) - 9216.0 * pow(xi, 7.0) +
            48384.0 * pow(xi, 5.0) - 80640.0 * xi * xi * xi + 30240.0 * xi;
    }

    else if (n == 10)
    {
        return 1024.0 * pow(xi, 10.0) - 23040.0 * pow(xi, 8.0) +
            161280.0 * pow(xi, 6.0) - 403200.0 * pow(xi, 4.0) +
            302400.0 * xi * xi - 30240.0;
    }

    else
    {
        return 0.0;
    }
}

static double I(double xi, int n)
{
    // first derivative of H(x) is H'(x)

    if (n == 1)
    {
        return 2.0;
    }

    else if (n == 2)
    {
        return 8.0 * xi;
    }

    else if (n == 3)
    {
        return 24.0 * xi * xi - 12.0;
    }

    else if (n == 4)
    {
        return 64.0 * pow(xi, 3.0) - 96.0 * xi;
    }

    else if (n == 5)
    {
        return 32.0 * 5.0 * pow(xi, 4.0) - 3.0 * 160.0 * pow(xi, 2.0) +
            120.0;
    }

    else if (n == 6)
    {
        return 6.0 * 64.0 * pow(xi, 5.0) - 4.0 * 480.0 * pow(xi, 3.0) +
            2.0 * 720.0 * xi;
    }

    else if (n == 7)
    {
        return 7.0 * 128.0 * pow(xi, 6.0) - 5.0 * 1344.0 * pow(xi, 4.0) +
            3.0 * 3360.0 * xi * xi - 1680.0;
    }

    else if (n == 8)
    {
        return 8.0 * 256.0 * pow(xi, 7.0) - 6.0 * 13584.0 * pow(xi, 5.0) +
            4.0 * 13440.0 * pow(xi, 3.0) - 2.0 * 13440.0 * xi;
    }

    else if (n == 9)
    {
        return 9.0 * 512.0 * pow(xi, 8.0) - 7.0 * 9216.0 * pow(xi, 6.0) +
            5.0 * 48384.0 * pow(xi, 4.0) - 3.0 * 80640.0 * xi * xi + 30240.0;
    }

    else if (n == 10)
    {
        return 10.0 * 1024.0 * pow(xi, 9.0) - 8.0 * 23040.0 * pow(xi, 7.0) +
            6.0 * 161280.0 * pow(xi, 5.0) - 4.0 * 403200.0 * pow(xi, 3.0) +
            2.0 * 302400.0 * xi;
    }

    else
    {
        return 0;
    }
}

static double J(double xi, int n)
{
    // second derivative of H(x) is H''(x)

    if (n == 2)
    {
        return 8.0;
    }

    else if (n == 3)
    {
        return 48.0 * xi * xi;
    }

    else if (n == 4)
    {
        return 3.0 * 64.0 * pow(xi, 2.0) - 96.0;
    }

    else if (n == 5)
    {
        return 4.0 * 32.0 * 5.0 * pow(xi, 3.0) - 2.0 * 3.0 * 160.0 * xi;
    }

    else if (n == 6)
    {
        return 5.0 * 6.0 * 64.0 * pow(xi, 4.0) - 3.0 * 4.0 * 480.0 * pow(xi, 2.0) +
            2.0 * 720.0;
    }

    else if (n == 7)
    {
        return 6.0 * 7.0 * 128.0 * pow(xi, 5.0) - 4.0 * 5.0 * 1344.0 * pow(xi, 3.0) +
            2.0 * 3.0 * 3360.0 * xi;
    }

    else if (n == 8)
    {
        return 7.0 * 8.0 * 256.0 * pow(xi, 6.0) - 5.0 * 6.0 * 13584.0 * pow(xi, 4.0) +
            3.0 * 4.0 * 13440.0 * pow(xi, 2.0) - 2.0 * 13440.0;
    }

    else if (n == 9)
    {
        return 8.0 * 9.0 * 512.0 * pow(xi, 7.0) - 6.0 * 7.0 * 9216.0 * pow(xi, 5.0) +
            4.0 * 5.0 * 48384.0 * pow(xi, 3.0) - 2.0 * 3.0 * 80640.0 * xi;
    }

    else if (n == 10)
    {
        return 9.0 * 10.0 * 1024.0 * pow(xi, 8.0) - 7.0 * 8.0 * 23040.0 * pow(xi, 6.0) +
            5.0 * 6.0 * 161280.0 * pow(xi, 4.0) - 3.0 * 4.0 * 403200.0 * pow(xi, 2.0) +
            2.0 * 302400.0;
    }

    else
    {
        return 0;
    }
}

static double Psi(double xi, int n)
{
    return N(n) * exp(-xi * xi / 2.0) * H(xi, n);
}

static void FindMinMax(
    std::vector<double> x,
    double& min,
    double& max)
{
    min = x[0];
    max = x[0];

    for (int i = 1; i < x.size(); i++)
    {
        if (x[i] < min)
            min = x[i];
        if (x[i] > max)
            max = x[i];
    }
}

static void DrawTitles(
    HDC hdc, RECT clientRect, const std::wstring& fTitle,
    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());

    // x-axis title: "x"

    std::wstring xTitle = L"xi";
    GetTextExtentPoint32W(hdc, xTitle.c_str(), (int)xTitle.length(), &sz);
    w = sz.cx;
    TextOutW(hdc, sx0 + (sx1 - sx0 - w) / 2, sy1 + 2 * h, xTitle.c_str(), (int)xTitle.length());

    // y-axis title: "?"
    
    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 double Fx1(double xi)
{
    double psi = Psi(xi, n);
    return xi * psi * psi;
}

static double Fx2(double xi)
{
    double psi = Psi(xi, n);
    return xi * xi * psi * psi;
}

static double Fx3(double xi)
{
    double psi = Psi(xi, n);
    return xi * xi * xi * psi * psi;
}

static double Fx4(double xi)
{
    double psi = Psi(xi, n);
    return pow(xi, 4.0) * psi * psi;
}

//  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)
{
    static HFONT hFont = nullptr;
    static HWND hCombo1 = nullptr;
    static HWND hCombo2 = nullptr;
    static HWND hCombo3 = nullptr;

    switch (message)
    {
    case WM_CREATE:
        CreateWindowEx(0, L"STATIC", L"draw:", WS_CHILD | WS_VISIBLE,
            10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
        hCombo1 = CreateWindowEx(0, L"COMBOBOX", nullptr,
            WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
            ES_AUTOVSCROLL, 120, 10, 120, 100, hWnd, nullptr, hInst, nullptr);
        CreateWindowEx(0, L"STATIC", L"average:", WS_CHILD | WS_VISIBLE,
            10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
        hCombo2 = CreateWindowEx(0, L"COMBOBOX", nullptr,
            WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
            ES_AUTOVSCROLL, 120, 40, 120, 150, hWnd, nullptr, hInst, nullptr);
        CreateWindowEx(0, L"STATIC", L"n:", WS_CHILD | WS_VISIBLE,
            10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC, hInst, NULL);
        hCombo3 = CreateWindowEx(0, L"COMBOBOX", nullptr,
            WS_CHILD | WS_VISIBLE | CBS_DROPDOWNLIST | WS_VSCROLL |
            ES_AUTOVSCROLL, 120, 70, 120, 150, hWnd, nullptr, hInst, nullptr);
        CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            10, 100, 80, 30, hWnd, (HMENU)IDC_COMPUTE_BUTTON, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Cancel", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            10, 140, 80, 30, hWnd, (HMENU)IDC_CANCEL_BUTTON, hInst, NULL);
        CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
            10, 180, 80, 30, hWnd, (HMENU)IDC_CLEAR_BUTTON, hInst, NULL);
        hFont = CreateFont(16, 0, 0, 0, FW_BOLD, FALSE, FALSE, FALSE,
            UNICODE, OUT_DEFAULT_PRECIS, CLIP_DEFAULT_PRECIS,
            DEFAULT_QUALITY, DEFAULT_PITCH | FF_SWISS, L"Courier New");
        SendMessage(hCombo1, WM_SETFONT, (WPARAM)hFont, 0);
        SendMessage(hCombo2, WM_SETFONT, (WPARAM)hFont, 0);
        SendMessage(hCombo3, WM_SETFONT, (WPARAM)hFont, 0);
        SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Hermite");
        SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"Psi");
        SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"|Psi^2|");
        SendMessage(hCombo1, CB_ADDSTRING, 0, (LPARAM)L"zeros");
        SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x");
        SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x^2");
        SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x^3");
        SendMessage(hCombo2, CB_ADDSTRING, 0, (LPARAM)L"x^4");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"1");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"2");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"3");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"4");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"5");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"6");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"7");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"8");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"9");
        SendMessage(hCombo3, CB_ADDSTRING, 0, (LPARAM)L"10");
        SendMessage(hCombo1, CB_SETCURSEL, 0, 0);
        SendMessage(hCombo2, CB_SETCURSEL, 0, 0);
        SendMessage(hCombo3, CB_SETCURSEL, 0, 0);
        ShowWindow(hWnd, SW_SHOWMAXIMIZED);
        buffer[0] = L'\0';
        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 IDC_COMPUTE_BUTTON:
        {
            if (hCombo3)
            {
                LRESULT result1 = SendMessage(hCombo1, CB_GETCURSEL, 0, 0);
                LRESULT result2 = SendMessage(hCombo2, CB_GETCURSEL, 0, 0);
                GetWindowText(hCombo3, buffer, 128);
                std::wstring nStr(buffer);
                n = stoi(nStr);
                double x0 = -10.0, x1 = 10.0, h = (x1 - x0) / 256;
                double xValue = -10.0, yValue = 0;

                xi.clear();
                fx.clear();

                if (result1 == 0)
                {
                    wcscpy_s(buffer, L"xi Versus Hermite Polynomial H(xi)");
                    fTitle = std::wstring(buffer);
                    wcscpy_s(buffer, L"H(xi, n)");
                    yTitle = std::wstring(buffer);

                    while (xValue <= x1)
                    {
                        xi.push_back(xValue);
                        yValue = H(xValue, n);
                        fx.push_back(yValue);
                        xValue += h;
                    }
                }

                else if (result1 == 1)
                {
                    wcscpy_s(buffer, L"xi Versus Wave Function Psi(xi, n)");
                    fTitle = std::wstring(buffer);
                    wcscpy_s(buffer, L"Psi(xi, n)");
                    yTitle = std::wstring(buffer);

                    while (xValue <= x1)
                    {
                        xi.push_back(xValue);
                        yValue = Psi(xValue, n);
                        fx.push_back(yValue);
                        xValue += h;
                    }
                }

                else if (result1 == 2)
                {
                    wcscpy_s(buffer, L"xi Versus |Psi^2|");
                    fTitle = std::wstring(buffer);
                    wcscpy_s(buffer, L"|Psi(xi, n)^2|");
                    yTitle = std::wstring(buffer);

                    while (xValue <= x1)
                    {
                        double psi = Psi(xValue, n);
                        xi.push_back(xValue);
                        yValue = psi * psi;
                        fx.push_back(yValue);
                        xValue += h;
                    }
                }

                GetClientRect(hWnd, &rect);
                rect.left = 300;
                draw = true;
                InvalidateRect(hWnd, &rect, TRUE);

                double integ = 0.0;

                if (result2 == 0)
                {
                    integ = Integration::SimpsonsRule(
                        1024, -10.0, 10.0, Fx1);
                    swprintf_s(buffer, L"xi * Psi^2(xi, n) = %lf", integ);
                }

                if (result2 == 1)
                {
                    integ = Integration::SimpsonsRule(
                        1024, -10.0, 10.0, Fx2);
                    swprintf_s(buffer, L"xi^2 * Psi^2(xi, n) = %lf", integ);
                    double integ1

                }

                if (result2 == 2)
                {
                    integ = Integration::SimpsonsRule(
                        1024, -10.0, 10.0, Fx3);
                    swprintf_s(buffer, L"xi^3 * Psi^2(xi, n) = %lf", integ);
                }

                if (result2 == 3)
                {
                    integ = Integration::SimpsonsRule(
                        1024, -10.0, 10.0, Fx4);
                    swprintf_s(buffer, L"xi^4 * Psi^2(xi, n) = %lf", integ);
                }

                MessageBox(hWnd, buffer, L"Information",
                    MB_OK | MB_ICONINFORMATION);
            }

            break;
        }
        case IDC_CLEAR_BUTTON:
            buffer[0] = '\0';
            draw = false;
            break;
        case IDC_CANCEL_BUTTON:
        case IDM_EXIT:
            DestroyWindow(hWnd);
            break;
        default:
            return DefWindowProc(hWnd, message, wParam, lParam);
        }
    }
    case WM_PAINT:
    {
        if (draw)
        {
            PAINTSTRUCT ps;
            POINT wPt = { };
            HDC hdc = BeginPaint(hWnd, &ps);
            double xMin = 0, yMin = 0;
            double xMax = 0, yMax = 0;
            FindMinMax(xi, xMin, xMax);
            FindMinMax(fx, yMin, yMax);
            double h = 0, pi = 0, plm = 0, theta = 0;
            float xSpan = (float)(xMax - xMin);
            float ySpan = (float)(yMax - yMin);
            float width = (float)(rect.right - rect.left + 1);
            float height = (float)(rect.bottom - rect.top - 32 + 1);
            float sx0 = 2.0f * width / 16.0f + 300;
            float sx1 = 14.0f * width / 16.0f + 300;
            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;
            int i = 0;
            float x = (float)xMin;
            float y = (float)yMax;
            DrawTitles(
                hdc, rect, fTitle, yTitle,
                (int)sx0, (int)sx1, (int)sy0, (int)sy1);
            px = (float)xi[0];
            py = (float)fx[0];
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            MoveToEx(hdc, (int)sx, (int)sy0, &wPt);

            while (i <= 8)
            {
                sx = xSlope * x + xInter;
                wPt.x = wPt.y = 0;
                MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
                LineTo(hdc, (int)sx, (int)sy1);
                char numberStr[128] = "";
                sprintf_s(numberStr, 128, "%5.4lf", x);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    numberStr,
                    (int)strlen(numberStr),
                    &size);
                RECT textRect = { };
                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)xi[0];
            py = (float)fx[0];
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            wPt.x = wPt.y = 0;
            MoveToEx(hdc, (int)sx, (int)sy, &wPt);

            for (size_t j = 1; j < (int)xi.size(); j++)
            {
                px = (float)xi[j];
                py = (float)fx[j];
                sx = xSlope * px + xInter;
                sy = ySlope * py + yInter;
                LineTo(hdc, (int)sx, (int)sy);
            }

            DeleteObject(bPenNew);
            bPenNew = NULL;

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

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

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

#pragma once
class Integration
{
public:
	static double TrapezoidalRule(
		int n, double a, double b, double(*fx)(double));
	static double SimpsonsRule(
		int n, double a, double b, double(*fx)(double));
	static void GenerateGaussianLegendreAbscissasAndWeights(
		int n, double x[], double w[]);
	static double GLIntegrate11(
		int n, double x[], double w[], double(*fx)(double));
	static double GLIntegrateAB(int n, double a, double b,
		double x[], double w[], double(*fx)(double));
};

#include "pch.h"
#include "Integration.h"
#include "OrthogonalPolynomials.h"
#include "Zeros.h"

double Integration::TrapezoidalRule(
    int n, double a, double b, double(*fx)(double))
{
    double h = (b - a) / n;
    double s = 0.5 * (fx(a) + fx(b));
    double x = a + h;

    for (int i = 1; i < n; i++)
    {
        s += fx(x);
        x += h;
    }

    return h * s;
}

double Integration::SimpsonsRule(
    int n, double a, double b, double(*fx)(double))
{
    double h = (b - a) / n;
    double h2 = 2.0 * h;
    double s = 0.0;
    double t = 0.0;
    double x = a + h;

    for (int i = 1; i < n; i += 2)
    {
        s += fx(x);
        x += h2;
    }

    x = a + h2;

    for (int i = 2; i < n; i += 2)
    {
        t += fx(x);
        x += h2;
    }

    return h * (fx(a) + 4 * s + 2 * t + fx(b)) / 3.0;
}

void Integration::GenerateGaussianLegendreAbscissasAndWeights(
    int n, double x[], double w[])
{
    char name[32] = { 0 };
    OrthogonalPolynomials op;
    std::vector<double> roots;
    Zeros::Zero("Legendre", 0.0, 0.0, n, roots);

    strcpy_s(name, sizeof(name), "Legendre");

    for (int i = 0; i < n; i++)
    {
        double xi = roots[i];
        double fd = op.g(name, xi, n);
        double x2 = 1.0 - xi * xi;

        x[i] = xi;
        w[i] = 2.0 / (x2 * fd * fd);
    }
}

double Integration::GLIntegrate11(
    int n, double x[], double w[], double(*fx)(double))
{
    double sum = 0.0;

    for (int i = 0; i < n; i++)
        sum += w[i] * fx(x[i]);

    return sum;
}

double Integration::GLIntegrateAB(int n, double a, double b,
    double x[], double w[], double(*fx)(double))
{
    double c = (b - a) / 2.0;
    double d = (b + a) / 2.0;
    double sum = 0.0;

    for (int i = 0; i < n; i++)
        sum += w[i] * fx(c * x[i] + d);

    return c * sum;
}