Tag: fresnel-integrals

// FresnelFunctions.cpp : Defines the entry point for the application.
// (c) Thursday, January 23, 2025, by James Pate Williams, Jr.
// Tables of Fresnel Integrals C(x) and S(x)
// References: "A Numerical Library in C for Scientists and Engineers"
// (c) 1995 by H. T. Lau also the National Bureau of Standards and
// https://nvlpubs.nist.gov/nistpubs/jres/102/3/j23mie.pdf
// https://digital.library.unt.edu/ark:/67531/metadc40301/m2/1/high_res_d/applmathser55_w.pdf

#include "framework.h"
#include "FresnelFunctions.h"
#include <chrono>
#include <cwchar>
#include <vector>

#define MAX_LOADSTRING 100

// resource definitions

#define IDC_EDIT1 201
#define IDC_EDIT2 202
#define IDC_EDIT3 203
#define IDC_EDIT4 204
#define IDC_EDIT5 205

#define IDC_STATIC1 301
#define IDC_STATIC2 302
#define IDC_STATIC3 303
#define IDC_STATIC4 304

#define IDC_BUTTON_COMPUTE 401
#define IDC_BUTTON_CLEAR 402

// 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[16834], line[128];

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

    // Perform application initialization:

    if (!InitInstance (hInstance, nCmdShow))
    {
        return FALSE;
    }

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

    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_FRESNELFUNCTIONS));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_FRESNELFUNCTIONS);
    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 fg_lau(double, double*, double*);

void fresnel_lau(double x, double* c, double* s)
{
	double absx, x3, x4, a, p, q, f, g, c1, s1;

	absx = fabs(x);
	if (absx <= 1.2) {
		a = x * x;
		x3 = a * x;
		x4 = a * a;
		p = (((5.47711385682687e-6 * x4 - 5.28079651372623e-4) * x4 +
			1.76193952543491e-2) * x4 - 1.99460898826184e-1) * x4 + 1.0;
		q = (((1.18938901422876e-7 * x4 + 1.55237885276994e-5) * x4 +
			1.09957215025642e-3) * x4 + 4.72792112010453e-2) * x4 + 1.0;
		*c = x * p / q;
		p = (((6.71748466625141e-7 * x4 - 8.45557284352777e-5) * x4 +
			3.87782123463683e-3) * x4 - 7.07489915144523e-2) * x4 +
			5.23598775598299e-1;
		q = (((5.95281227678410e-8 * x4 + 9.62690875939034e-6) * x4 +
			8.17091942152134e-4) * x4 + 4.11223151142384e-2) * x4 + 1.0;
		*s = x3 * p / q;
	}
	else if (absx <= 1.6) {
		a = x * x;
		x3 = a * x;
		x4 = a * a;
		p = ((((-5.68293310121871e-8 * x4 + 1.02365435056106e-5) * x4 -
			6.71376034694922e-4) * x4 + 1.91870279431747e-2) * x4 -
			2.07073360335324e-1) * x4 + 1.00000000000111e0;
		q = ((((4.41701374065010e-10 * x4 + 8.77945377892369e-8) * x4 +
			1.01344630866749e-5) * x4 + 7.88905245052360e-4) * x4 +
			3.96667496952323e-2) * x4 + 1.0;
		*c = x * p / q;
		p = ((((-5.76765815593089e-9 * x4 + 1.28531043742725e-6) * x4 -
			1.09540023911435e-4) * x4 + 4.30730526504367e-3) * x4 -
			7.37766914010191e-2) * x4 + 5.23598775598344e-1;
		q = ((((2.05539124458580e-10 * x4 + 5.03090581246612e-8) * x4 +
			6.87086265718620e-6) * x4 + 6.18224620195473e-4) * x4 +
			3.53398342767472e-2) * x4 + 1.0;
		*s = x3 * p / q;
	}
	else if (absx < 1.0e15) {
		fg_lau(x, &f, &g);
		a = x * x;
		a = (a - floor(a / 4.0) * 4.0) * 1.57079632679490;
		c1 = cos(a);
		s1 = sin(a);
		a = (x < 0.0) ? -0.5 : 0.5;
		*c = f * s1 - g * c1 + a;
		*s = -f * c1 - g * s1 + a;
	}
	else
		*c = *s = ((x > 0.0) ? 0.5 : -0.5);
}

void fg_lau(double x, double* f, double* g)
{
	double absx, c, s, c1, s1, a, xinv, x3inv, c4, p, q;

	absx = fabs(x);
	if (absx <= 1.6) {
		fresnel_lau(x, &c, &s);
		a = x * x * 1.57079632679490;
		c1 = cos(a);
		s1 = sin(a);
		a = (x < 0.0) ? -0.5 : 0.5;
		p = a - c;
		q = a - s;
		*f = q * c1 - p * s1;
		*g = p * c1 + q * s1;
	}
	else if (absx <= 1.9) {
		xinv = 1.0 / x;
		a = xinv * xinv;
		x3inv = a * xinv;
		c4 = a * a;
		p = (((1.35304235540388e1 * c4 + 6.98534261601021e1) * c4 +
			4.80340655577925e1) * c4 + 8.03588122803942e0) * c4 +
			3.18309268504906e-1;
		q = (((6.55630640083916e1 * c4 + 2.49561993805172e2) * c4 +
			1.57611005580123e2) * c4 + 2.55491618435795e1) * c4 + 1.0;
		*f = xinv * p / q;
		p = ((((2.05421432498501e1 * c4 + 1.96232037971663e2) * c4 +
			1.99182818678903e2) * c4 + 5.31122813480989e1) * c4 +
			4.44533827550512e0) * c4 + 1.01320618810275e-1;
		q = ((((1.01379483396003e3 * c4 + 3.48112147856545e3) * c4 +
			2.54473133181822e3) * c4 + 5.83590575716429e2) * c4 +
			4.53925019673689e1) * c4 + 1.0;
		*g = x3inv * p / q;
	}
	else if (absx <= 2.4) {
		xinv = 1.0 / x;
		a = xinv * xinv;
		x3inv = a * xinv;
		c4 = a * a;
		p = ((((7.17703249365140e2 * c4 + 3.09145161574430e3) * c4 +
			1.93007640786716e3) * c4 + 3.39837134926984e2) * c4 +
			1.95883941021969e1) * c4 + 3.18309881822017e-1;
		q = ((((3.36121699180551e3 * c4 + 1.09334248988809e4) * c4 +
			6.33747155851144e3) * c4 + 1.08535067500650e3) * c4 +
			6.18427138172887e1) * c4 + 1.0;
		*f = xinv * p / q;
		p = ((((3.13330163068756e2 * c4 + 1.59268006085354e3) * c4 +
			9.08311749529594e2) * c4 + 1.40959617911316e2) * c4 +
			7.11205001789783e0) * c4 + 1.01321161761805e-1;
		q = ((((1.15149832376261e4 * c4 + 2.41315567213370e4) * c4 +
			1.06729678030581e4) * c4 + 1.49051922797329e3) * c4 +
			7.17128596939302e1) * c4 + 1.0;
		*g = x3inv * p / q;
	}
	else {
		xinv = 1.0 / x;
		a = xinv * xinv;
		x3inv = a * xinv;
		c4 = a * a;
		p = ((((2.61294753225142e4 * c4 + 6.13547113614700e4) * c4 +
			1.34922028171857e4) * c4 + 8.16343401784375e2) * c4 +
			1.64797712841246e1) * c4 + 9.67546032967090e-2;
		q = ((((1.37012364817226e6 * c4 + 1.00105478900791e6) * c4 +
			1.65946462621853e5) * c4 + 9.01827596231524e3) * c4 +
			1.73871690673649e2) * c4 + 1.0;
		*f = (c4 * (-p) / q + 0.318309886183791) * xinv;
		p = (((((1.72590224654837e6 * c4 + 6.66907061668636e6) * c4 +
			1.77758950838030e6) * c4 + 1.35678867813756e5) * c4 +
			3.87754141746378e3) * c4 + 4.31710157823358e1) * c4 +
			1.53989733819769e-1;
		q = (((((1.40622441123580e8 * c4 + 9.38695862531635e7) * c4 +
			1.62095600500232e7) * c4 + 1.02878693056688e6) * c4 +
			2.69183180396243e4) * c4 + 2.86733194975899e2) * c4 + 1.0;
		*g = (c4 * (-p) / q + 0.101321183642338) * x3inv;
	}
}

double S(double z)
{
	int n = 768;
	double a = 0.0;
	double b = z;
	double h = (b - a) / n;
	double h2 = 2.0 * h;
	double s = 0.0;
	double t = 0.0;
	double x = a + h;
	double pi = 4.0 * atan(1.0);
	
	for (int i = 1; i < n; i += 2)
	{
		s += sin(pi * x * x / 2.0);
		x += h2;
	}

	x = a + h2;

	for (int i = 2; i < n; i += 2)
	{
		t += sin(pi * x * x / 2.0);
		x += h2;
	}

	double endA = 0.0;
	double endB = sin(pi * z * z / 2.0);
	return h * (endA + 4 * s + 2 * t + endB) / 3.0;
}

void fresnel_nbs(double x, double* c, double* s)
{
	double pi = 4.0 * atan(1.0);

	*c = x - pi * pi * pow(x, 5.0) / 40.0 -
		pi * pi * pi * pi * pow(x, 9.0) / 3456.0;
	// the following computation is in error
	*s = -pi / 6.0 + pi * pi * pi * pow(x, 7.0) / 336.0 -
		pow(pi, 5.0) * pow(x, 11.0) / 42240.0;
	// resort to Simpson's Rule with n = 16384
	*s = S(x);
} 

//
//  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)
{
	double lt = 0.0, rt = 0.0;
	int cors = 1, nPts = 0;
	WCHAR* endPt;
	static HWND hwndEdit1 = 0, hwndEdit2 = 0;
	static HWND hwndEdit3 = 0, hwndEdit4 = 0, hwndEdit5 = 0;

	switch (message)
	{
	case WM_CREATE:
	{
		CreateWindowEx(0, L"STATIC", L"Lt Point:", WS_CHILD | WS_VISIBLE,
			10, 10, 80, 20, hWnd, (HMENU)IDC_STATIC1, GetModuleHandle(NULL), NULL);
		CreateWindowEx(0, L"STATIC", L"Rt Point:", WS_CHILD | WS_VISIBLE,
			10, 40, 80, 20, hWnd, (HMENU)IDC_STATIC2, GetModuleHandle(NULL), NULL);
		CreateWindowEx(0, L"STATIC", L"n Points:", WS_CHILD | WS_VISIBLE,
			10, 70, 80, 20, hWnd, (HMENU)IDC_STATIC3, GetModuleHandle(NULL), NULL);
		CreateWindowEx(0, L"STATIC", L"C 1, S, 2:", WS_CHILD | WS_VISIBLE,
			10, 100, 80, 20, hWnd, (HMENU)IDC_STATIC4, GetModuleHandle(NULL), NULL);

		hwndEdit1 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
			100, 10, 200, 20, hWnd, (HMENU)IDC_EDIT1, GetModuleHandle(NULL), NULL);
		hwndEdit2 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
			100, 40, 200, 20, hWnd, (HMENU)IDC_EDIT2, GetModuleHandle(NULL), NULL);
		hwndEdit3 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
			100, 70, 200, 20, hWnd, (HMENU)IDC_EDIT3, GetModuleHandle(NULL), NULL);
		hwndEdit4 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER,
			100, 100, 200, 20, hWnd, (HMENU)IDC_EDIT4, GetModuleHandle(NULL), NULL);
		hwndEdit5 = CreateWindowEx(0, L"EDIT", NULL, WS_CHILD | WS_VISIBLE | WS_BORDER
			| ES_MULTILINE, 100, 170, 400, 400, hWnd, (HMENU)IDC_EDIT5,
			GetModuleHandle(NULL), NULL);

		CreateWindowEx(0, L"BUTTON", L"Compute", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
			10, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_COMPUTE, GetModuleHandle(NULL), NULL);
		CreateWindowEx(0, L"BUTTON", L"Clear", WS_CHILD | WS_VISIBLE | BS_PUSHBUTTON,
			220, 130, 80, 30, hWnd, (HMENU)IDC_BUTTON_CLEAR, GetModuleHandle(NULL), NULL);

		SetDlgItemText(hWnd, IDC_EDIT1, L"-0.62");
		SetDlgItemText(hWnd, IDC_EDIT2, L"+0.62");
		SetDlgItemText(hWnd, IDC_EDIT3, L"25");
		SetDlgItemText(hWnd, IDC_EDIT4, L"1");
	}
    case WM_COMMAND:
        {
            int wmId = LOWORD(wParam);

            // Parse the menu selections:
			
			switch (wmId)
			{
			case IDC_BUTTON_COMPUTE:
			{
				buffer[0] = { '/0' };
				GetWindowText(hwndEdit1, line, 127);
				lt = wcstod(line, &endPt);

				if (!(lt >= -0.62 && lt <= 0.50))
				{
					MessageBox(hWnd, L"Must be -0.62 <= lt <= 0.50", L"Warning",
						MB_OK | MB_ICONWARNING);
					return 0;
				}

				GetWindowText(hwndEdit2, line, 127);
				rt = wcstod(line, &endPt);

				if (!(rt >= 0.51 && rt <= 0.62))
				{
					MessageBox(hWnd, L"Must be 0.51 <= rt <= 0.62", L"Warning",
						MB_OK | MB_ICONWARNING);
					return 0;
				}

				nPts = GetDlgItemInt(hWnd, IDC_EDIT3, FALSE, FALSE);

				if (nPts < 1 || nPts > 25)
				{
					MessageBox(hWnd, L"n-points must be >= 1 & <= 25", L"Warning",
						MB_OK | MB_ICONWARNING);
					return 0;
				}

				cors = GetDlgItemInt(hWnd, IDC_EDIT4, FALSE, FALSE);

				if (!(cors == 1 || cors == 2))
				{
					MessageBox(hWnd, L"C(x) = 1 or S(x) = 2", L"Warning",
						MB_OK | MB_ICONWARNING);
					return 0;
				}

				if (cors == 1)
				{
					wsprintf(line, L"Fresnel Cosine Integral\r\n\r\n");
					wcscpy_s(buffer, 16383, line);
					wsprintf(line, L"x\t\tC(x) Lau\t\tC(x) NBS\t\t%% Difference\r\n");
					wcscat_s(buffer, 16383, line);
				}

				else if (cors == 2)
				{
					wsprintf(line, L"Fresnel Sine Integral\r\n\r\n");
					wcscpy_s(buffer, 16383, line);
					wsprintf(line, L"x\t\tS(x) Lau\t\tS(x) Integral\t%% Difference\r\n");
					wcscat_s(buffer, 16383, line);
				}

				double deltaX = (rt - lt) / nPts, x = lt;

				for (int i = 0; i <= nPts; i++)
				{
					double c1 = 0.0, c2 = 0.0, s1 = 0.0, s2 = 0.0;
					double percentDifference = 0.0;
					double cx = x, sx = x;

					fresnel_lau(cx, &c1, &s1);
					fresnel_nbs(sx, &c2, &s2);

					if (cors == 1)
					{
						percentDifference = fabs(c1 - c2) / fabs((c1 + c2) / 2.0);
						swprintf_s(
							line, L"%13.6lf\t%13.10lf\t%13.10lf\t%13.10lf\r\n",
							cx, c1, c2, 100.0 * percentDifference);
						wcscat_s(buffer, 16383, line);
					}

					else if (cors == 2)
					{
						percentDifference = fabs(s1 - s2) / fabs((s1 + s2) / 2.0);
						swprintf_s(
							line, L"%13.10lf\t%13.10lf\t%13.10lf\t%13.10lf\r\n",
							sx, s1, s2, 100.0 * percentDifference);
						wcscat_s(buffer, 16383, line);
					}

					x += deltaX;
				}

				SetDlgItemText(hWnd, IDC_EDIT5, buffer);
			}
			break;
			case IDC_BUTTON_CLEAR:
			{
				wcscpy_s(buffer, 16383, L"\0");
				SetDlgItemText(hWnd, IDC_EDIT5, buffer);
			}
			break;
            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);
            // TODO: Add any drawing code that uses hdc here...
            EndPaint(hWnd, &ps);
        }
        break;
	case WM_SIZE:
		if (hwndEdit5)
		{
			RECT rcClient;
			GetClientRect(hWnd, &rcClient);
			SetWindowPos(hwndEdit5, NULL, 10, 170, rcClient.right - 20,
				rcClient.bottom - 20, SWP_NOZORDER);
		}
		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;
}