//-------------------------------------------------------------
// Title:		Nonlinear Pendulum Auxiliary Functions
// Author:		Tony Kim 
// Date:		2006 Nov 10 
// Description: The workings of NPendulum
//
//		The heart of this program is the following
//		differential equation:
//
//		d2x/dt2 + (1/q)*dx/dt + sin(x) = gcos(w_d*t)
//
//		with the parameters:
//			q:		measure of damp-freeness (quality)
//			g:		forcing amplitude
//			w_d:	forcing frequency
//
//-------------------------------------------------------------
#pragma once

#include <vector>
#include <math.h>
#include "FieldEngine.h"
#include "MyCoord.h"

//--------------------------------------------------
// Display parameters
//--------------------------------------------------
#define RADIUS	600
#define MARGIN	15
std::vector<MyCoord> g_coord;

//--------------------------------------------------
// Math constants
//--------------------------------------------------
#define M_PI	3.141592654

//--------------------------------------------------
// Simulation parameters
//--------------------------------------------------
#define q 2
#define g 6
#define w_d 2
#define dt 0.1

//--------------------------------------------------
// Calculate pendulum state
//--------------------------------------------------
struct Pendulum
{
	double theta;
	double theta_dot;
} g_pendulum;

std::vector<Pendulum> g_timeseries;

//--------------------------------------------------
// Update the simulation parameters
//--------------------------------------------------
void Update( int MyTime )
{
	g_pendulum.theta += dt*g_pendulum.theta_dot;
	g_pendulum.theta_dot += dt*(-(double)1/(double)q*g_pendulum.theta_dot - sin( g_pendulum.theta ) - (double)g*cos(w_d*dt*MyTime));
	if( g_pendulum.theta > 2*M_PI )
	{
		g_pendulum.theta -= 4*M_PI;
	}
	else if( g_pendulum.theta < -2*M_PI )
	{
		g_pendulum.theta += 4*M_PI;
	}
	g_timeseries.push_back( g_pendulum );
}

//--------------------------------------------------
// Draw pendulum
//--------------------------------------------------
void DrawPendulum( HDC hDC )
{
	HPEN hPen = CreatePen( PS_SOLID, 2, RGB(255,255,255) );
	HBRUSH hBrush = CreateSolidBrush( RGB(0,0,0) );

	double theta = g_pendulum.theta;
	SelectObject( hDC, hPen );

	MoveToEx( hDC, wx(0,g_coord[0]), wy(0,g_coord[0]), NULL );
	LineTo( hDC, wx(RADIUS*sin(theta),g_coord[0]), wy(-RADIUS*cos(theta),g_coord[0]));

	int size = 80;
	SelectObject( hDC, hBrush );
	Ellipse( hDC, wx(RADIUS*sin(theta)-size,g_coord[0]), wy(-RADIUS*cos(theta)-size,g_coord[0]), 
			 wx(RADIUS*sin(theta)+size,g_coord[0]), wy(-RADIUS*cos(theta)+size,g_coord[0]) );
	
	DeleteObject( hBrush );
	DeleteObject( hPen );
}	

//--------------------------------------------------
// Draw phase portrait
//--------------------------------------------------
void DrawPortrait( HDC hDC )
{
	HPEN hPen = CreatePen( PS_SOLID, 1, RGB(255,255,255) );
	HPEN hGreyPen = CreatePen( PS_SOLID, 1, RGB(200,200,200) );
	SelectObject( hDC, hPen );

	MoveToEx( hDC, wx(g_coord[1].GetMinX(),g_coord[1]),wy(0,g_coord[1]), NULL );
	LineTo( hDC, wx(g_coord[1].GetMaxX(),g_coord[1]),wy(0,g_coord[1]));

	MoveToEx( hDC, wx(0,g_coord[1]),wy(g_coord[1].GetMaxY(),g_coord[1]), NULL );
	LineTo( hDC, wx(0,g_coord[1]),wy(g_coord[1].GetMinY(),g_coord[1]));

	double size = 0.1;
	Ellipse( hDC, wx(g_pendulum.theta-size,g_coord[1]),wy(g_pendulum.theta_dot-size,g_coord[1]),
		wx(g_pendulum.theta+size,g_coord[1]),wy(g_pendulum.theta_dot+size,g_coord[1]));

	SelectObject( hDC, hGreyPen );
	MoveToEx( hDC, wx(g_timeseries[0].theta,g_coord[1]),wy(g_timeseries[0].theta_dot,g_coord[1]), NULL );
	for(int i = 1; i < (int)g_timeseries.size(); i++)
	{
		if( ((g_timeseries[i].theta > 0) && (g_timeseries[i-1].theta > 0)) || 
			((g_timeseries[i].theta < 0) && (g_timeseries[i-1].theta < 0)) || 
			(abs(g_timeseries[i].theta) < M_PI) )
			LineTo( hDC, wx(g_timeseries[i].theta,g_coord[1]),wy(g_timeseries[i].theta_dot,g_coord[1]));
		else
			MoveToEx( hDC, wx(g_timeseries[i].theta,g_coord[1]),wy(g_timeseries[i].theta_dot,g_coord[1]), NULL );
	}

	DeleteObject( hGreyPen );
	DeleteObject( hPen );
}

//--------------------------------------------------
// Draw poincare section
//--------------------------------------------------
void DrawPoincare( HDC hDC )
{
	HPEN hPen = CreatePen( PS_SOLID, 1, RGB(255,255,255) );
	HPEN hGreyPen = CreatePen( PS_SOLID, 1, RGB(200,200,200) );
	SelectObject( hDC, hPen );

	MoveToEx( hDC, wx(g_coord[2].GetMinX(),g_coord[2]),wy(0,g_coord[2]), NULL );
	LineTo( hDC, wx(g_coord[2].GetMaxX(),g_coord[2]),wy(0,g_coord[2]));

	MoveToEx( hDC, wx(0,g_coord[2]),wy(g_coord[2].GetMaxY(),g_coord[2]), NULL );
	LineTo( hDC, wx(0,g_coord[2]),wy(g_coord[2].GetMinY(),g_coord[2]));

	SelectObject( hDC, hGreyPen );

	double size = 0.1;
	for(int i = 0; i < (int)g_timeseries.size(); i++)
	{
		if( i % (int)(2*(double)M_PI/((double)w_d*(double)dt)) == 0 )
		{
			Ellipse( hDC, wx(g_timeseries[i].theta-size,g_coord[2]),wy(g_timeseries[i].theta_dot-size,g_coord[2]),
						  wx(g_timeseries[i].theta+size,g_coord[2]),wy(g_timeseries[i].theta_dot+size,g_coord[2]));
		}
	}

	DeleteObject( hGreyPen );
	DeleteObject( hPen );
}

//--------------------------------------------------
// Return which coordinate was selected
//	-1 if none
//--------------------------------------------------
int ClickedCoord( int wx, int wy )
{
	int select = -1;
	for(int i=0; i < (int)g_coord.size(); i++)
	{
		// Check if it is in the xbound
		if( (wx >= g_coord[i].GetAnchX()) && (wx <= g_coord[i].GetAnchX() + g_coord[i].GetWidth()) )
		{
			if( (wy >= g_coord[i].GetAnchY()) && (wy <= g_coord[i].GetAnchY() + g_coord[i].GetHeight()) )
			{
				select = i;
			}
		}
	}
	return select;
}