Category: Bessel Functions of the First and Second Kind

// SecondOrderLinearODE.cpp (c) Friday March 13, 2026
// by James Pate Williams, Jr., BA, BS, MSwE, PhD
// Solves the Bessel Equation of the First Kind for Order 0
// Using the power series found in the reference below - 
// References: https://math.libretexts.org/Bookshelves/Calculus/Calculus_(OpenStax)/17%3A_Second-Order_Differential_Equations/17.04%3A_Series_Solutions_of_Differential_Equations
// https://mathworld.wolfram.com/BesselDifferentialEquation.html

#include "framework.h"
#include <Windows.h>
#include "SecondOrderLinearODE.h"
#include <string>
#include <vector>

#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 fTitle[128];
double chi[1025], eta[1025], coeff[65];

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

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

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

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

static double Factorial(int n)
{
    double factorial = 1.0;

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

    return factorial;
}

static void ComputeCoefficients()
{
    double a0 = 1.0;

    for (int k = 0; k <= 64; k++)
    {
        double kFactorial = Factorial(k);

        coeff[k] = a0 * pow(-1.0, k) /
            (pow(2.0, 2.0 * k) * pow(kFactorial, 2.0));
    }
}

static void ComputePoints()
{
    double chi0 = -15.0;
    double deltaChi = 30.0 / 1024;
 
    for (int i = 0; i <= 1024; i++)
    {
        double sum = 0.0;

        for (int k = 64; k >= 0; k--)
            sum += pow(chi0, 2 * k) * coeff[k];

        chi[i] = chi0;
        eta[i] = sum;

        chi0 += deltaChi;
    }
}

static void FindMinMax(
    double x[],
    double& min,
    double& max)
{
    min = x[0];
    max = x[0];

    for (int i = 1; i <= 1024; 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,
    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
    GetTextExtentPoint32(hdc, fTitle.c_str(), (int)fTitle.length(), &sz);
    int w = sz.cx;
    int h = sz.cy;
    TextOut(hdc, (width - w) / 2, h, fTitle.c_str(), (int)fTitle.length());

    // x-axis title: "x"
    std::wstring xTitle = L"x";
    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: "p(x)"
    std::wstring yTitle = L"J0(x)";
    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);
}

//
//  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:
        ComputeCoefficients();
        ComputePoints();
        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);

        ComputePoints();
        double xMin = 0, yMin = 0;
        double xMax = 0, yMax = 0;
        FindMinMax(chi, xMin, xMax);
        FindMinMax(eta, yMin, yMax);

        double h = 0, pi = 0, plm = 0, theta = 0;
        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 = { 0 };
        int i = 0;
        float x = (float)xMin;
        float y = (float)yMax;
        wcscpy_s(fTitle, 128, L"Bessel Function of the First Kind Order 0");
        DrawTitles(hdc, rect, fTitle, (int)sx0, (int)sx1, (int)sy0, (int)sy1);


        px = (float)chi[0];
        py = (float)eta[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 = { 0 };
            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 nPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

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

        HRGN clipRegion = CreateRectRgn(
            (int)sx0, (int)sy0,              // Left, Top
            (int)(sx1), (int)(sy1)           // Right, Bottom
        );

        // Apply clipping region

        SelectClipRgn(hdc, clipRegion);


        px = (float)chi[0];
        py = (float)eta[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 <= 1024; j++)
        {

            px = (float)chi[j];
            py = (float)eta[j];

            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        DeleteObject(nPenNew);
        nPenNew = 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;
}