//---------------------------------------------------------------------------
// File: universe.cxx
// Written by: Michael Main and Allison Brown
// This program illustrates a first use of structs
// for an animation of elementary particles.
#include <cctype>
#include <cmath>
#include <iostream>
#include <ctime>
#include <vector>
#include "graphics.h"
using namespace std;
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// Global constants regarding graphics and animation
const int S = 500;                // Size of a graphics window
const int FRAMES_PER_SECOND = 30; // For the animation
// Physical constants
const double WMAX = 5.0e-2;             // Maximum world coordinate
const double WMIN = -WMAX;              // Minimum world coordinate
const double SIMRATE = 1.0e-5;          // Seconds time in each frame (s)
const double E_CHARGE = -1.602e-19;     // Charge of electron (Coulombs)
const double P_CHARGE = +1.602e-19;     // Charge of proton (Coulombs)
const double E_MASS = 9.109e-31;        // Mass of electron (kg)
const double P_MASS = 1.673e-27;        // Mass of proton (kg)
const double K = 8.988e9;               // Coulomb force constant (Nm^2/C^2)
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
class particle
{
public:
    // Get and set functions are implemented inline
    double get_x() const { return x; }  // gets the value of the x coordinate
    double get_y() const {return y; }   // gets the value of the y coordinate
    void set_properties(
	double new_x, double new_y,
	double new_vx, double new_vy,
	double new_m, double new_c
	)
	{
	    x = new_x;
	    y = new_y;
	    vx = new_vx;
	    vy = new_vy;
	    mass = new_m;
	    charge = new_c;
	}
    
    // Compute the acceleration that this particle feels from some other
    // particle in the x- or y- dimension.  The answer, which takes into
    // account only the electrical force (not gravity), is in m/s^2.
    double accx(particle p2);
    double accy(particle p2);

    // Draw this particle in an already open square graphics screen:
    void draw(int pixel_size, double wmin, double wmax);

    // Simulate the change in this particle over t seconds, given that
    // the vector contains all the other particles in the universe:
    void simtime(double t, const vector<particle>& universe);
private:
    double x, y;   // position in world coordinates (in meters)
    double vx, vy; // how fast particle is going (in m/s)
    double mass;   // in kg
    double charge; // in Coulombs 
};
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
// Function prototypes
// Interact with mouse clicks to add particles to the universe:
void add_particles(vector<particle>& universe);

// Interact with keystrokes to change the current seconds_per_frame:
void change_rate(double& seconds_per_frame);

// pixel(v, v0, v1, pmax) converts v from an interval that ranges from v0 to v1
// into an integer interval that ranges from 0 to pmax. Note that when v = v0,
// the answer is always 0, and when v = v1, the answer is always pmax.
// Example 1: Suppose an x-coordinate called x ranges from -2 (on the left
// side of the screen) to +3 (on the right side of the screen), and
// we want to figure out the corresponding pixel x coordinate for a screen
// that is 400 pixel wide. Then call pixel(x, -2, 3, 400).
// Example 2: Suppose a y-coordinate called y ranges from +4 (on the top
// of the screen) to -2 (on the bottom of the screen), and
// we want to figure out the corresponding pixel y coordinate for a screen
// that is 300 pixels tall. Then call pixel(y, 4, -2, 300). Note that in this
// case v0 > v1 since we want bigger coordinates at the top of the screen.
int pixel(double v, double v0, double v1, int pmax);

// world(p, v0, v1, pmax) converts p from an interval that ranges from 0 to
// pmax into a double number interval that ranges from v0 to v1.
double world(int p, double v0, double v1, int pmax);
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
int main( )
{
    vector<particle> universe;
    vector<particle>::size_type i;
    double seconds_per_frame = SIMRATE;
    
    initwindow(S,S,"Universe",0,0,true);

    while (true)
    {
	// 1. Calculations
	change_rate(seconds_per_frame);
	add_particles(universe);
	for (i = 0; i < universe.size( ); ++i)
	{
	    universe[i].simtime(SIMRATE, universe);
	}
	
	// 2. Draw the next frame of animation
	clearviewport( );
	for (i = 0; i < universe.size( ); ++i)
	{
	    universe[i].draw(S, WMIN, WMAX);
	}
	
	// 3. Swap and delay
	swapbuffers( );
	delay(1000/FRAMES_PER_SECOND);
    }
}
//---------------------------------------------------------------------------

//-----------------------------------------------------------------------------
// START of the particle member functions
double particle::accx(particle p2)
{
    double dx = (x - p2.x);
    double dy = (y - p2.y);
    double denominator = pow(dx*dx + dy*dy, 1.5)*mass;

    if (denominator == 0)
    {
	return 0;
    }
    else
    {
	return K*charge*p2.charge*dx/denominator;
    }
}

double particle::accy(particle p2)
{
    double dx = (x - p2.x);
    double dy = (y - p2.y);
    double denominator = pow(dx*dx + dy*dy, 1.5)*mass;

    if (denominator == 0)
    {
	return 0;
    }
    else
    {
	return K*charge*p2.charge*dy/denominator;
    }
}

void particle::draw(int pixel_size, double wmin, double wmax)
{
    int radius = int(2.5 + 4e9*pow(mass, 0.33));
    int px = pixel(x, wmin, wmax, pixel_size);
    int py = pixel(y, wmax, wmin, pixel_size);

    if (charge < 0)
    {
	setfillstyle(SOLID_FILL, RED);
    }
    else if (charge > 0)
    {
	setfillstyle(SOLID_FILL, BLUE);
    }
    else
    {
	setfillstyle(SOLID_FILL, GREEN);
    }
    fillellipse(px, py, radius, radius);
}

void particle::simtime(double t, const vector<particle>& universe)
{
    vector<particle>::size_type i;

    // Update the velocity:
    for (i = 0; i < universe.size( ); ++i)
    {
	// Change the velocity of the particle that activated this method
	// by using the force that's felt from universe[i].
	vx += accx(universe[i])*t;
	vy += accy(universe[i])*t;
    }

    // Update the position:
    x += vx*t;
    y += vy*t;
}
// END of the particle member functions
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
void add_particles(vector<particle>& universe)
{
    int px, py;
    double wx, wy;
    particle p;
    int down_time, hold_time;
    double vx;
    
    if (ismouseclick(WM_LBUTTONDOWN))
    {   // Make a new electron
	down_time = time(NULL);
	getmouseclick(WM_LBUTTONDOWN, px, py);
	wx = world(px, WMIN, WMAX, S);
	wy = world(py, WMAX, WMIN, S);
	while (!ismouseclick(WM_LBUTTONUP))
	{
	    delay(10);
	}
	hold_time = time(NULL) - down_time;
	clearmouseclick(WM_LBUTTONUP);
	vx = (hold_time/1000.0)*(WMAX-WMIN)/SIMRATE;
	p.set_properties(wx, wy, vx, 0, E_MASS, E_CHARGE);
	universe.push_back(p);
    }

    if (ismouseclick(WM_RBUTTONDOWN))
    {   // Make a new electron
	down_time = time(NULL);
	getmouseclick(WM_RBUTTONDOWN, px, py);
	wx = world(px, WMIN, WMAX, S);
	wy = world(py, WMAX, WMIN, S);
	while (!ismouseclick(WM_RBUTTONUP))
	{
	    delay(10);
	}
	hold_time = time(NULL) - down_time;
	clearmouseclick(WM_RBUTTONUP);
	vx = (hold_time/1000.0)*(WMAX-WMIN)/SIMRATE;
	p.set_properties(wx, wy, vx, 0, P_MASS, P_CHARGE);
	universe.push_back(p);
    }
}
//---------------------------------------------------------------------------

//---------------------------------------------------------------------------
void change_rate(double& seconds_per_frame)
{
    static double old_rate = SIMRATE;
    
    if (kbhit( ))
    {
	switch (toupper(getch( )))
	{
	case '-':
	    seconds_per_frame /= 1.2;
	    break;
	case '+':
	    seconds_per_frame *= 1.2;
	    break;
	case 'P':
	    if (seconds_per_frame != 0)
	    {
		old_rate = seconds_per_frame;
		seconds_per_frame = 0;
	    }
	    break;
	case 'R':
	    if (seconds_per_frame == 0)
	    {
		seconds_per_frame = old_rate;
	    }
	    break;
	}
    }
}
//---------------------------------------------------------------------------

//-----------------------------------------------------------------------------
int pixel(double v, double v0, double v1, int pmax)
{
    return int(((v - v0)/(v1 - v0)) * pmax);
}
//-----------------------------------------------------------------------------

//----------------------------------------------------------------------------
double world(int p, double v0, double v1, int pmax)
{
    double fraction = double(p)/double(pmax);
    double distance_from_v0 = fraction*(v1 - v0);
    
    return v0 + distance_from_v0;
}
//----------------------------------------------------------------------------