Numerical Explorations

class Figure
{
public:
	static void ComputeFigure(
		bool eight,
		double& xi2, double& eta2,
		std::vector<double>& radius,
		std::vector<double>& xi,
		std::vector<double>& eta,
		std::vector<double>& vix,
		std::vector<double>& viy,
		std::vector<double>& energy1,
		double (*f)(double),
		double (*g)(double));
};

#include "Figure.h"
#include <algorithm>
#include <cmath>
#include <iterator>
#include <set>
#include <vector>

// Function to compute intersection of two sorted containers of doubles with tolerance
template <typename InputIt1, typename InputIt2, typename OutputIt>
void set_intersection_with_tolerance(InputIt1 first1, InputIt1 last1,
	InputIt2 first2, InputIt2 last2,
	OutputIt d_first, double tolerance) {
	while (first1 != last1 && first2 != last2) {
		double a = *first1;
		double b = *first2;

		if (std::fabs(a - b) <= tolerance) {
			// Values are considered equal within tolerance
			*d_first++ = /*a; // or */ (a + b) / 2.0; //if you want averaged value
			++first1;
			++first2;
		}
		else if (a < b - tolerance) {
			++first1;
		}
		else {
			++first2;
		}
	}
}

static double f(double xi)
{
	return xi * tan(xi);
}

static double g(double xi)
{
	return -xi * 1.0 / tan(xi);
}

static void Intersection(
	bool eight, std::vector<double>& intersection)
{
	double del = 0.0, radius2 = 0.0;
	double tolerance = 0.0;
	double (*h)(double);
	int ilimit = 0;
	std::set<double> eta, rad;

	if (eight)
	{
		ilimit = 10000;
		radius2 = 1.0;
		tolerance = 1.0e-1;
		h = f;
	}

	else
	{
		ilimit = 10000;
		radius2 = 4.0;
		tolerance = 1.0e-1;
		h = g;
	}

	del = radius2 / ilimit;
	intersection.clear();

	for (size_t i = 0; i < ilimit; i++)
	{
		double xi0 = i * del;

		if (eight && xi0 >= 0.65)
		{
			double et = h(xi0);
			double r2 = xi0 * xi0 + et * et;
			
			if (fabs(r2 - radius2) < tolerance)
			{
				eta.insert(et);
				rad.insert(r2);
			}
		}

		else if (!eight && xi0 >= 1.8125)
		{
			double et = h(xi0);
			double r2 = xi0 * xi0 + et * et;

			if (fabs(r2 - radius2) < tolerance)
			{
				eta.insert(et);
				rad.insert(r2);
			}
		}
	}

	if (eight)
	{
		size_t count = 0, index = 0;

		for (double val : eta)
		{
			if (index == eta.size() - 1)
			{
				intersection.push_back(val);
				break;
			}

			count++;
			index++;
		}
	}

	else
	{
		int count = 24, index = 0;

		for (double val : eta)
		{
			if (index == eta.size() - count)
			{
				intersection.push_back(val);
				break;
			}

			index++;
		}
	}

	/*set_intersection_with_tolerance(
		eta.begin(), eta.end(),
		rad.begin(), rad.end(),
		std::back_inserter(intersection),
		tolerance);*/
}

void Figure::ComputeFigure(
	bool eight,
	double& xi2, double& eta2,
	std::vector<double>& radius,
	std::vector<double>& xi,
	std::vector<double>& eta,
	std::vector<double>& vix,
	std::vector<double>& viy,
	std::vector<double>& intersection,
	double (*f)(double),
	double (*g)(double))
{
	double xi0 = 0.0;
	double xi1 = 3.5;
	double eta0 = 0.0, eta1 = 0.0;
	double (*h)(double);

	xi.clear();
	eta.clear();

	if (eight)
		h = f;
	else
	{
		xi0 = 1.5;
		h = g;
	}

	eta0 = h(xi0);
	eta1 = h(xi1);

	double deltaXi = (xi1 - xi0) / 10000.0;
	int count = 0;

	vix.clear();
	viy.clear();

	for (int j = 0; j <= 10000; j++)
	{
		double x = j * deltaXi;
		double hx = h(x);
		double vx = x * x + hx * hx;

		if (eight && x >= xi0 && x <= xi1)
		{
			if (count == 0 && x >= xi0)
			{
				xi.push_back(x);
				eta.push_back(hx);
				count = 1;
			}

			else
			{
				xi.push_back(x);
				eta.push_back(hx);
			}

			for (int k = 0; k < 2; k++)
			{
				double r = radius[k];

				vix.push_back(xi[xi.size() - 1]);
				viy.push_back(r * r);
			}
		}

		else if (!eight && x >= xi0 && x <= xi1)
		{
			if (count == 0 && x >= xi0)
			{
				xi.push_back(x);
				eta.push_back(hx);
				count = 1;
			}

			else if (count == 1)
			{
				xi.push_back(x);
				eta.push_back(hx);
			}

			if (xi.size() > 0)
			{
				for (int k = 0; k < 1; k++)
				{
					double r = radius[k];

					vix.push_back(xi[xi.size() - 1]);
					viy.push_back(r * r);
				}
			}
		}

		if (eight && xi[j] >= 0.0 && eta[j] >= 3.5)
			break;
		
		else if (!eight &&
			xi.size() > 0 &&
			xi[xi.size() - 1] > 1.5 &&
			eta.size() > 0 &&
			eta[eta.size() - 1] >= 3.5)
			break;
	}
	
	Intersection(eight, intersection);
}

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

#include "framework.h"
#include "SchiffChapter2Fig8and9.h"
#include <time.h>
#include <algorithm>
#include <vector>
#include "Figure.h"

#define MAX_LOADSTRING 100

typedef struct tagPoint2d
{
    double x, y;
} Point2d, * PPoint2d;

// Global Variables:
HINSTANCE hInst;                                // current instance
WCHAR szTitle[MAX_LOADSTRING];                  // The title bar text
WCHAR szWindowClass[MAX_LOADSTRING];            // the main window class name
WCHAR line[8192], text[8192];                   // wide character buffers
WCHAR title[65536];                             // window title
bool eight;                                     // true = plot figure 8
double xi2, eta2;                               // energy cordinates
std::vector<Point2d> points;                    // plotting 2d points
std::vector<double> radius;                     // radius vector { 1.0, 2.0 }
std::vector<double> xi;                         // Greek x-coordinate
std::vector<double> eta;                        // Greek y-coordinate
std::vector<double> V0;                         // potential of the well
std::vector<double> vix;                        // potential x-coordinate
std::vector<double> viy;                        // potential y-coordinate
std::vector<double> energy1;                    // energy eigenvalues

// 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_PTR CALLBACK    DrawEtaDialog(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_SCHIFFCHAPTER2FIG8AND9, szWindowClass, MAX_LOADSTRING);
    MyRegisterClass(hInstance);

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

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

    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_SCHIFFCHAPTER2FIG8AND9));
    wcex.hCursor        = LoadCursor(nullptr, IDC_ARROW);
    wcex.hbrBackground  = (HBRUSH)(COLOR_WINDOW+1);
    wcex.lpszMenuName   = MAKEINTRESOURCEW(IDC_SCHIFFCHAPTER2FIG8AND9);
    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 fx(double x)
{
    return x * tan(x);
}

static double gx(double x)
{
    return -x * 1.0 / tan(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)
{
    switch (message)
    {
    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_FIGURE8:
            {
                eight = true;
                text[0] = L'\0';

                clock_t clock0 = clock();
                radius.push_back(1.0);
                radius.push_back(2.0);
                Figure::ComputeFigure(
                    eight, xi2, eta2, radius, xi, eta, vix, viy, energy1, fx, gx);
                clock_t clock1 = clock() - clock0;
                double runtime = (double)clock1 / CLOCKS_PER_SEC;
                swprintf_s(line, 8192, L"Runtime in Seconds = %lf\r\n", runtime);
                wcscat_s(text, 8192, line);

                points.clear();
               
                for (size_t j = 0; j < xi.size(); j++)
                {
                    Point2d pt = { 0 };
                    pt.x = xi[j];
                    pt.y = eta[j];
                    points.push_back(pt);
                }
                
                for (size_t i = 0; i < energy1.size(); i++)
                {
                    swprintf_s(line, 8192, 
                        L"E[%zu] = %lf\r\n", i, energy1[i]);
                    wcscat_s(text, 8192, line);
                }

                MessageBox(hWnd, text, L"Energy Information",
                    MB_OK | MB_ICONINFORMATION);

                DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_ETA_DIALOG), hWnd, DrawEtaDialog);
                break;
            }
            case IDM_FIGURE9:
            {
                eight = false;
                text[0] = L'\0';

                clock_t clock0 = clock();
                radius.push_back(2.0);
                Figure::ComputeFigure(
                    eight, xi2, eta2, radius, xi, eta, vix, viy, energy1, fx, gx);
                clock_t clock1 = clock() - clock0;
                double runtime = (double)clock1 / CLOCKS_PER_SEC;
                swprintf_s(line, 8192, L"Runtime in Seconds = %lf\r\n", runtime);
                wcscat_s(text, 8192, line);

                points.clear();
                
                for (size_t j = 0; j < xi.size(); j++)
                {
                    Point2d pt = { 0 };
                    pt.x = xi[j];
                    pt.y = eta[j];
                    points.push_back(pt);
                }

                for (size_t i = 0; i < energy1.size(); i++)
                {
                    swprintf_s(line, 8192,
                        L"E[%zu] = %lf\r\n", i, energy1[i]);
                    wcscat_s(text, 8192, line);
                }

                MessageBox(hWnd, text, L"Energy Information",
                    MB_OK | MB_ICONINFORMATION);

                DialogBox(hInst, MAKEINTRESOURCE(IDD_DRAW_ETA_DIALOG), hWnd, DrawEtaDialog);
                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_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;
}

static void FindMinMax(
    std::vector<Point2d>& points,
    double& xMin, double& xMax,
    double& yMin, double& yMax)
{
    // uses global 2D double points structure

    xMin = yMin = DBL_MAX;
    xMax = yMax = DBL_MIN;

    for (size_t i = 0; i < points.size(); i++)
    {
        Point2d pt = points[i];
        double x = pt.x;
        double y = pt.y;

        if (x < xMin)
            xMin = x;
        if (x > xMax)
            xMax = x;
        if (y < yMin)
            yMin = y;
        if (y > yMax)
            yMax = y;
    }
}

static void DrawFormattedText(HDC hdc, char text[], RECT rect)
{
    // Draw the text with formatting options
    DrawTextA(hdc, text, -1, &rect, DT_SINGLELINE | DT_NOCLIP);
}

static void DrawQuarterCircleArc(
    HDC hdc,
    float xSlope, float ySlope,
    float xInter, float yInter,
    float radius, bool topToRight)
{
    auto mapX = [&](float x)
        {
            return (int)lroundf(xSlope * x + xInter);
        };

    auto mapY = [&](float y)
        {
            return (int)lroundf(ySlope * y + yInter);
        };

    int x1 = mapX(-radius);
    int y1 = mapY(+radius);
    int x2 = mapX(+radius);
    int y2 = mapY(-radius);

    int left = min(x1, x2);
    int right = max(x1, x2);
    int top = min(y1, y2);
    int bottom = max(y1, y2);

    int xTop = mapX(0.0f);
    int yTop = mapY(radius);

    int xRight = mapX(radius);
    int yRight = mapY(0.0f);

    if (topToRight)
    {
        Arc(hdc, left, top, right, bottom,
            xTop, yTop, xRight, yRight);
    }

    else
    {
        Arc(hdc, left, top, right, bottom,
            xRight, yRight, xTop, yTop);
    }
}

INT_PTR CALLBACK DrawEtaDialog(
    HWND hDlg, UINT message, WPARAM wParam, LPARAM lParam)
{
    UNREFERENCED_PARAMETER(lParam);
    switch (message)
    {
    case WM_INITDIALOG:
        SetWindowText(hDlg, title);
        return (INT_PTR)TRUE;
    case WM_COMMAND:
        if (LOWORD(wParam) == IDOK || LOWORD(wParam) == IDCANCEL)
        {
            EndDialog(hDlg, LOWORD(wParam));
            return (INT_PTR)TRUE;
        }
        break;
    case WM_PAINT:
        double h = 0, pi = 0, plm = 0, theta = 0;
        double xMax = 0, xMin = 0, yMax = 0, yMin = 0;
        FindMinMax(points, xMin, xMax, yMin, yMax);
        float xSpan = (float)(xMax - xMin);
        float ySpan = (float)(yMax - yMin);
        RECT rect = { };
        GetClientRect(hDlg, &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, xInter, ySlope, yInter;
        xSlope = (sx1 - sx0) / xSpan;
        xInter = (float)(sx0 - xSlope * xMin);
        ySlope = (sy0 - sy1) / ySpan;
        yInter = (float)(sy0 - ySlope * yMax);
        float px = 0, py = 0, sx = 0, sy = 0;
        float vx = 0, vy = 0;
        PAINTSTRUCT ps;
        POINT wPt = { };
        HDC hdc = BeginPaint(hDlg, &ps);
        int i = 0;
        float x = (float)xMin;
        float y = (float)yMax;
        px = x;
        py = y;
        sx = xSlope * px + xInter;
        sy = ySlope * py + yInter;
        MoveToEx(hdc, (int)sx, (int)sy0, &wPt);
        char buffer[128] = { };

        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);

            sprintf_s(buffer, "%5.4lf", x);
            SIZE size = { };
            GetTextExtentPoint32A(
                hdc,
                buffer,
                (int)strlen(buffer),
                &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, buffer, 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)
            {
                sprintf_s(buffer, "%5.3lf", y);
                SIZE size = { };
                GetTextExtentPoint32A(
                    hdc,
                    buffer,
                    (int)strlen(buffer),
                    &size);
                RECT textRect = { };
                textRect.left = (long)(sx0 - size.cx - size.cx / 5.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, buffer, textRect);
            }

            y += deltaY;
            i++;
        }

        HGDIOBJ bPenNew = NULL;
        HGDIOBJ hPenOld = NULL;

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

        HRGN hrgn = CreateRectRgn((int)sx0, (int)sy0, (int)sx1, (int)sy1);
        
        // Select the clipping region into the DC
        if (SelectClipRgn(hdc, hrgn) == ERROR) {
            MessageBox(hDlg, L"Failed to select clip region", 
                L"Error", MB_ICONERROR);
            return (INT_PTR)FALSE;
        }

        SelectClipRgn(hdc, hrgn);

        px = (float)points[0].x;
        py = (float)points[0].y;
        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 < points.size(); j++)
        {
            px = (float)points[j].x;
            py = (float)points[j].y;
            sx = xSlope * px + xInter;
            sy = ySlope * py + yInter;
            LineTo(hdc, (int)sx, (int)sy);
        }

        float radius = 0.0f;

        if (eight)
            radius = 1.0f;
        else
            radius = 2.0f;

        DrawQuarterCircleArc(
            hdc, xSlope, ySlope,
            xInter, yInter, radius, false);

        if (eight)
        {
            radius = 2.0f;

            DrawQuarterCircleArc(
                hdc, xSlope, ySlope,
                xInter, yInter, radius, false);
        }

        SelectObject(hdc, hPenOld);
        DeleteObject(bPenNew);

        return (INT_PTR)FALSE;
    }

    return (INT_PTR)FALSE;
}