/*
 *	(c) 2001, Michel Beaudouin-Lafon, mbl@lri.fr
 *
 *	Simple particle system
 */

#include <stdlib.h>
#include <stdio.h>
#include <GL/glut.h>
#include <math.h>

inline float random (float mean, float variance)
{
  double res;
  res = 2.0 * (double) rand() / (double) RAND_MAX - 0.5;     // (-1 <= res <= 1)
  return res*variance + mean;
}

// -------- Vector --------

class Vector {
public:
	float x, y, z;
	
	Vector ();
	Vector (float vx, float vy, float vz);
	
	Vector& operator += (const Vector& v);
	
	float Norm ();
	void Normalize ();
};

Vector::Vector ()
{
	x = y = z = 0.0;
}

Vector::Vector (float vx, float vy, float vz)
{
	x = vx;
	y = vy;
	z = vz;
}

Vector&
Vector::operator += (const Vector& v)
{
	x += v.x;
	y += v.y;
	z += v.z;
	return *this;
}

float
Vector::Norm ()
{
	return sqrtf (x*x + y*y + z*z);
}

void
Vector::Normalize ()
{
	float norm = Norm ();
	if (norm != 0.0) {
		x /= norm;
		y /= norm;
		z /= norm;
	}
}

Vector
operator - (const Vector& v2, const Vector& v1)
{
	return Vector (v2.x - v1.x, v2.y - v1.y, v2.z - v1.z);
}

Vector
operator * (float a, const Vector& v)
{
	return Vector (a * v.x, a * v.y, a * v.z);
}


// -------- Particle --------

class Particle {
protected:
	float mass;
	Vector position;
	Vector speed;
	Vector force;
	int fixed;
public:
	Particle (float px, float py, bool pfixed = false);
	virtual ~Particle ();
	
	const Vector& Position ()	{ return position; }
	const Vector& Speed ()		{ return speed; }
	const Vector& Force ()		{ return force; }
	
	void ResetForce ();
	void AddForce (const Vector& v);
	
	virtual void Step (float dt);
	virtual void Render ();
	
	virtual void Dump ();
};


Particle::Particle (float px, float py, bool pfixed)
: position(), speed(), force()
{
	mass = 1.0;
	position.x = px;
	position.y = py;
	fixed = pfixed;
}

Particle::~Particle ()
{
	// rien
}

void
Particle::ResetForce ()
{
	force = Vector ();
}

void
Particle::AddForce (const Vector& v)
{
	force += v;
}

void
Particle::Step (float dt)
{
	if (! fixed) {
		speed += (dt / mass) * force;
		position += dt * speed;
	}
}

void
Particle::Render ()
{
	glVertex3f (position.x, position.y, position.z);
}

void
Particle::Dump ()
{
	printf ("%f %f %f\n", position.x, position.y, position.z);
}


// -------- System --------

#define MAX_PARTICLES 10

class System {
protected:
	Particle* particles[MAX_PARTICLES];
	int maxp;
	
public:
	System ();
	~System ();
	
	void AddParticle (Particle* p);
	void Spring (Particle* p1, Particle* p2, float k, float d0);
	
	virtual void Simulate (float dt);
	virtual void Render ();
	
	virtual void Dump ();
};

System::System ()
{
	maxp = MAX_PARTICLES;
	for (int p = 0; p < maxp; p++)
		particles[p] = new Particle (200.0 - p * 2., 200.00, p == 0);
}

System::~System ()
{
	for (int p = 0; p < maxp; p++)
		delete particles[p];
}

void
System::Simulate (float dt)
{
	int p;
	for (p = 0; p < maxp; p++) {
		particles[p]->ResetForce();
		/* gravity */
		particles[p]->AddForce (Vector (0.0, -9.81, 0.0));
		/* spring between each pair of successive particles */
		if (p > 0)
			Spring (particles[p-1], particles[p], 10.0, 2.0);
	}
	for (p = 0; p < maxp; p++)
		particles[p]->Step(dt);

}

void
System::Spring (Particle* p1, Particle* p2, float k, float d0)
{
	Vector u = p2->Position () - p1->Position ();
	float d = u.Norm ();
	u.Normalize ();
	
	p1->AddForce (( k*(d-d0)) * u);
	p2->AddForce ((-k*(d-d0)) * u);
}

void
System::Render ()
{
	glBegin (GL_LINE_STRIP);
	for (int p = 0; p < maxp; p++) {
		particles[p]->Render();
	}
	glEnd ();
}

void
System::Dump ()
{
	printf ("%d particles\n", maxp);
	for (int p = 0; p < maxp; p++) {
		particles[p]->Dump();
		if (p > 10) {
			printf ("...\n");
			break;
		}
	}
}


// -------- Main Program --------

float dt = 0.05;
System* psystem = 0;

void
init ()
{
	if (psystem)
		delete psystem;
	
	psystem = new System;
}

void
idle ()
{
	psystem->Simulate (dt);
	glutPostRedisplay ();
}

void
display ()
{
	glClear (GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);
	psystem->Render ();
	glutSwapBuffers ();
}

void
keyboard (unsigned char key, int x, int y)
{
	switch (key) {
	case 27:
	case 'q':
		exit(0);
		break;
	case ' ':
		init();
		glutPostRedisplay();
		break;
	case 'd':
		psystem->Dump ();
	}
}

void 
reshape(int w, int h)
{
	glViewport(0, 0, w, h);
	glMatrixMode(GL_PROJECTION);
	glLoadIdentity();
	glOrtho(0.0, (GLdouble) w, 0.0, (GLdouble) h, -1.0, 1.0);
	glMatrixMode(GL_MODELVIEW);
	glLoadIdentity();
}

int 
main(int argc, char **argv)
{
	glutInit (&argc, argv);
	glutInitWindowSize (512, 256);
	glutInitDisplayMode (GLUT_DOUBLE | GLUT_RGB | GLUT_DEPTH);
	glutCreateWindow("particle system");
	
	init();
	
	glutReshapeFunc(reshape);
	glutDisplayFunc(display);
	glutKeyboardFunc(keyboard);
	glutIdleFunc (idle);
	
	glutMainLoop();
	return 0;
}