Category: Elementary Physics

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