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