// bulk.h  (Class Header)
// 2D Numerical Simulation Program for
// Collision of a Bulk and a Bullet
// Copyright (C) 1988,2000 Yukinori NAGATANI
//

#ifndef _bulk_h
#define _bulk_h

#include <math.h>
#include <iostream.h>
#include <stdio.h>
#include <stdlib.h>		// atof, atoi

typedef	double	Real;

void Error (char * szMessage)
{
	cerr << szMessage;
	exit (1);
}

//// Virtual Objects ////

class Object {
public:
	Real	x,	y;
	
	// New : memory allocation
	virtual Object * copyNew (void) const = 0;
	
	// get force at point (x,y) from this object
	virtual Object * getForce (Real x, Real y, Real &xF, Real &yF) const = 0;
	
	// time evolution by force and by itself.
	virtual Object * setTimeEvolution (Real xAc, Real yAc, Real dt) = 0;
	
	// print out for operator << 
	virtual ostream& print (ostream& s) const = 0;
	
	// print out for PostScript
	virtual ostream& put (ostream& s) const = 0;
	
};

ostream & operator << (ostream& s, const Object &O)
{
	O.print (s);
	return	s;
}

//// Definition of Objects ////

class PlanetObj : public Object {
public:
	Real	vx,	vy;
	Real	m;
	int	cc;
	
	PlanetObj (Real X, Real Y, Real VX, Real VY, Real M, int CC = 1)
	{
		x	= X;
		y	= Y;
		vx	= VX;
		vy	= VY;
		m	= M;
		cc	= CC;
	}
	PlanetObj (const PlanetObj &O)
	{
		x	= O.x;
		y	= O.y;
		vx	= O.vx;
		vy	= O.vy;
		m	= O.m;
		cc	= O.cc;
	}
	
	
	virtual PlanetObj * copyNew (void) const
	{
		return	new PlanetObj (*this);
	}
	virtual Object * getForce (Real x0, Real y0, Real &xF, Real &yF) const
	{
			Real	dx = x - x0,
				dy = y - y0;
			Real	d2 = dx*dx + dy*dy;
			Real	ff = m / sqrt(d2*d2*d2);
			
			xF = ff*dx;
			yF = ff*dy;
	}
	virtual Object * setTimeEvolution (Real xF, Real yF, Real dt)
	{
			vx += xF * dt;
			vy += yF * dt;
			
			x += vx * dt;
			y += vy * dt;
	}

	virtual ostream & print (ostream &s) const
	{
		s << "[";
		s << "x:"	<< x	<<  "  "
		  << "y:"	<< y	<<  "  "
		  << "vx:"	<< vx	<<  "  "
		  << "vy:"	<< vy	<<  "  "
		  << "m:"	<< m	<<  "  "
		  << "cc:"	<< cc;
		s << "]";
		
		return	s;
	}
	virtual ostream & put (ostream &s) const
	{
		s <<  ' '
		  << x	<< ' '
		  << y	<< ' '
		  << vx	<< ' '
		  << vy	<< ' '
		  << cc	<< " P";
		
		return	s;
	}
};

class FixedObj : public Object {
public:
	Real	m;
	int	cc;
	
	FixedObj (Real X, Real Y, Real M, int CC = 1)
	{
		x	= X;
		y	= Y;
		m	= M;
		cc	= CC;
	}
	FixedObj (const FixedObj &O)
	{
		x	= O.x;
		y	= O.y;
		m	= O.m;
		cc	= O.cc;
	}
	
	
	virtual FixedObj * copyNew (void) const
	{
		return	new FixedObj (*this);
	}
	virtual Object * getForce (Real x0, Real y0, Real &xF, Real &yF) const
	{
			Real	dx = x - x0,
				dy = y - y0;
			Real	d2 = dx*dx + dy*dy;
			Real	ff = m / sqrt(d2*d2*d2);
			
			xF = ff*dx;
			yF = ff*dy;
	}
	virtual Object * setTimeEvolution (Real xF, Real yF, Real dt)
	{
		// nothing to do
	}

	virtual ostream & print (ostream &s) const
	{
		s << "[";
		s << "x:"	<< x	<<  "  "
		  << "y:"	<< y	<<  "  "
		  << "m:"	<< m	<<  "  "
		  << "cc:"	<< cc;
		s << "]";
		
		return	s;
	}
	virtual ostream & put (ostream &s) const
	{
		s <<  ' '
		  << x	<< ' '
		  << y	<< ' '
		  << 0	<< ' '
		  << 0	<< ' '
		  << cc	<< " P";
		
		return	s;
	}
};

class AccelObj : public Object {
public:
	Real	vx,	vy;
	Real	m;
	int	cc;
	Real	Accel;
	
	AccelObj (Real X, Real Y, Real VX, Real VY, Real M, Real ACCEL, int CC = 1)
	{
		x	= X;
		y	= Y;
		vx	= VX;
		vy	= VY;
		m	= M;
		cc	= CC;
		Accel	= ACCEL;
	}
	AccelObj (const AccelObj &O)
	{
		x	= O.x;
		y	= O.y;
		vx	= O.vx;
		vy	= O.vy;
		m	= O.m;
		cc	= O.cc;
		Accel	= O.Accel;
	}
	
	
	virtual AccelObj * copyNew (void) const
	{
		return	new AccelObj (*this);
	}
	virtual Object * getForce (Real x0, Real y0, Real &xF, Real &yF) const
	{
		Real	dx = x - x0,
			dy = y - y0;
		Real	d2 = dx*dx + dy*dy;
		Real	ff = m / sqrt(d2*d2*d2);
		
		xF = ff*dx;
		yF = ff*dy;
	}
	virtual Object * setTimeEvolution (Real xF, Real yF, Real dt)
	{
		Real	velocity = sqrt (vx*vx + vy*vy),
			xTan = vx / velocity,
			yTan = vy / velocity;
		
		vx += (xF + Accel*xTan) * dt;
		vy += (yF + Accel*yTan) * dt;
		
		x += vx * dt;
		y += vy * dt;
	}

	virtual ostream & print (ostream &s) const
	{
		s << "[";
		s << "x:"	<< x		<<  "  "
		  << "y:"	<< y		<<  "  "
		  << "vx:"	<< vx		<<  "  "
		  << "vy:"	<< vy		<<  "  "
		  << "m:"	<< m		<<  "  "
		  << "Accel:"	<< Accel	<<  "  "
		  << "cc:"	<< cc;
		s << "]";
		
		return	s;
	}
	virtual ostream & put (ostream &s) const
	{
		s <<  ' '
		  << x	<< ' '
		  << y	<< ' '
		  << vx	<< ' '
		  << vy	<< ' '
		  << cc	<< " P";
		
		return	s;
	}
};

class BindObj : public Object {
public:
	Real	vx,	vy;
	Real	m;
	int	cc;
	
	BindObj (Real X, Real Y, Real VX, Real VY, Real M, int CC = 1)
	{
		x	= X;
		y	= Y;
		vx	= VX;
		vy	= VY;
		m	= M;
		cc	= CC;
	}
	BindObj (const BindObj &O)
	{
		x	= O.x;
		y	= O.y;
		vx	= O.vx;
		vy	= O.vy;
		m	= O.m;
		cc	= O.cc;
	}
	
	
	virtual BindObj * copyNew (void) const
	{
		return	new BindObj (*this);
	}
	virtual Object * getForce (Real x0, Real y0, Real &xF, Real &yF) const
	{
			Real	dx = x - x0,
				dy = y - y0;
			Real	rr = dx*dx + dy*dy;
			if (rr > 4) {
				xF = 0.0;
				yF = 0.0;
			} else {
				Real	ff = + (rr - 1) * exp(-4*rr);
				xF = ff*dx;
				yF = ff*dy;
			}
	}
	virtual Object * setTimeEvolution (Real xF, Real yF, Real dt)
	{
			vx += xF * dt / m;
			vy += yF * dt / m;
			
			x += vx * dt;
			y += vy * dt;
	}

	virtual ostream & print (ostream &s) const
	{
		s << "[";
		s << "x:"	<< x	<<  "  "
		  << "y:"	<< y	<<  "  "
		  << "vx:"	<< vx	<<  "  "
		  << "vy:"	<< vy	<<  "  "
		  << "m:"	<< m	<<  "  "
		  << "cc:"	<< cc;
		s << "]";
		
		return	s;
	}
	virtual ostream & put (ostream &s) const
	{
		s <<  ' '
		  << x	<< ' '
		  << y	<< ' '
		  << vx	<< ' '
		  << vy	<< ' '
		  << cc	<< " P";
		
		return	s;
	}
};

//// Node and List

class Node {
	friend class List;
	friend ostream & operator << (ostream& s, const List &L);
private:
	Node	*next;
	Object	*O;
public:
	Node (const Object &toNew, Node *node = NULL)
	{
		next  = node;
		O = toNew.copyNew ();
	}
	void free (void)
	{
		delete	O;
	}
};

class List {
private:
	Node		*root;
public:
	List (void)
	{
		root = NULL;
	}
	List (const List &L)
	{
		root	= NULL;
		
		for (Node *node = L.root ; node ; node = node->next)
			add (*(node->O));
	}
	List & add (const Object &toNew)
	{
		Node	*NewNode = new Node (toNew, root);
		if (NewNode == NULL)
			Error ("Memery Shotage in Add List.\n");
		
		root = NewNode;
		return	*this;
	}
	List & add (const List &L)
	{
		for (Node *node = L.root ; node ; node = node->next)
			add (*(node->O));
		
		return	*this;
	}

	List & addTail (const Object &toNew)
	{
		Node	*NewNode = new Node (toNew, NULL);
		if (NewNode == NULL)
			Error ("Memery Shotage in Add List to Tail.\n");
		
		if (root != NULL) {
			Node	*node = root;
			for (; node->next ; node = node->next);
			node->next = NewNode;
		} else {
			root = NewNode;
		}
		return	*this;
	}
	List & addTail (const List &L)
	{
		for (Node *node = L.root ; node ; node = node->next)
			addTail (*(node->O));
		
		return	*this;
	}
	List & freeAll (void)
	{
		Node	*node,	*nodeNext;
		for (node = root ; node ; node = nodeNext) {
			nodeNext = node->next;
			node->free ();
		}
		root = NULL;
		return	*this;
	}

	int nSize (void) const
	{
		int	n = 0;
		for (Node *node = root ; node ; node = node->next)
			n++;
		return	n;
	}
	friend ostream & operator << (ostream& s, const List &L)
	{
		s << '[';
		for (const Node *node = L.root ; node ; node = node->next)
			s << *(node->O) << ',' << '\n';
		s << ']';
		return	s;
	}

	const List & put (ostream &s = cout) const
	{
		for (const Node *node = root ; node ; node = node->next)
			node->O->put (s);
		s << endl;
		
		return	*this;
	}
	
	const List & getForce (const Node *node, Real &xF, Real &yF) const
	{
		xF = 0;
		yF = 0;
		for (const Node *no = root ; no ; no = no->next) {
			
        		if (node == no)
				continue;
			
			Real	xF1,	yF1;
			no->O->getForce (node->O->x, node->O->y, xF1, yF1);
			xF += xF1;
			yF += yF1;
		}
		return	*this;
	}
	List & step (Real dt, Real rG = 1)
	{
		for (Node *node = root ; node ; node = node->next) {
			
			Real	xF = 0,
				yF = 0;
			getForce (node, xF, yF);
			xF *=  rG;
			yF *=  rG;
			
			node->O->setTimeEvolution (xF, yF, dt);
        	}
	}
};


//// Put PostScript ////

void PutPrologue (ostream &s = cout)
{
	s <<	"%!PS-Adobe-2.0\n"
		"%%Creator: bulk (C) 1998,2000 Y.Nagatani\n"
		"%%Title: Collision of a Bulk and a Bullet\n"
		"%%Pages: (attend)\n"
		"%%DocumentPaperSizes: A4\n"
		"%%BoundingBox: 0 0 600 600\n"
		"%%Orientation: Portrait\n"
		"%%EndComments\n"
		"\n"
		"/mm {72 25.4 div mul} def\n"
		"/xPaperSize 210 mm def		% A3 size\n"
		"/yPaperSize 297 mm def\n"
		"\n"
		"/DefaultState save def\n"
		"20 dict begin\n"
		"\n"
		"/xMax 30 def	% Matrix Size\n"
		"/yMax 30 def\n"
		"\n"
		"/xSize 500 def	% Region Size (pt)\n"
		"/ySize 500 def\n"
		"\n"
		"/HELV {/Helvetica findfont exch scalefont setfont} def\n"
		"\n"
		"/T { %% x y point Text  T  -\n"
		"	4 2 roll	%% point Text x y\n"
		"	moveto		%% point Text\n"
		"	exch		%% Text point\n"
		"	Scale div HELV	%% Text\n"
		"	show\n"
		"} bind def\n"
		"\n"
		"/L { %% x0 y0 x1 y1  L  -\n"
		"	4 2 roll\n"
		"	moveto lineto stroke\n"
		"} bind def\n"
		"\n"
		"/Init {			% initialize func.\n"
		"	0 setlinewidth\n"
		"	0 0 0 setrgbcolor\n"
		"	0 0 moveto\n"
		"	xPaperSize 0 lineto\n"
		"	xPaperSize yPaperSize lineto\n"
		"	0 yPaperSize lineto\n"
		"	closepath fill stroke\n"
		"	\n"
		"	xPaperSize 2 div\n"
		"	yPaperSize 2 div translate\n"
		"	\n"
		"	1 1 1 setrgbcolor\n"
		"	-100 100 moveto 60 HELV\n"
		"	cvi 30 string cvs dup stringwidth pop neg 0 rmoveto\n"
		"	show\n"
		"	\n"
		"	/xScale xSize xMax div def \n"
		"	/yScale ySize yMax div def \n"
		"	xScale yScale scale\n"
		"	\n"
		"} bind def\n"
		"\n"
		"/COL { %% strength  COL  -\n"
		"	/v exch def\n"
		"	v 0 lt {/v 0 def} if\n"
		"	v 1 gt {/v 1 def} if\n"
		"	1 v neg add 0.7 mul 1 1 sethsbcolor\n"
		"} bind def\n"
		"\n"
		"/P { %% x y  vx vy  type  P  -\n"
		"	3 1 roll\n"
		"	dup mul 2 1 roll\n"
		"	dup mul add sqrt 20 mul COL\n"
		"	-1 add 0.1 mul 1 add 2 div\n"
		"	newpath 0 360 arc fill stroke\n"
		"} bind def\n"
		"\n"
		"%%EndProlog\n";
}

void PutPageHead (int nT, ostream &s = cout)
{
	s <<	"%%Page: " << nT << endl
	  <<	"gsave "   << nT << " Init" << endl;
}

void PutPageTail (ostream &s = cout)
{
	s <<	"grestore showpage" << endl;
}

void PutEpilogue (ostream &s = cout)
{
	s <<	"%%Trailer\n"
		"%%EOF" << endl;
}

//// Test Codes ////

void Test1 (void)
{
	cout << "-------------------------------- Start of Test1" << endl;
	
	List	LS;
	cout << "LS:"		<< LS		<< endl;
	cout << "LSmenbers:"	<< LS.nSize ()	<< endl;
	
	LS.add (PlanetObj (1,2, 3,4, 5));
	LS.add (PlanetObj (6,7, 8,9, 10));
	LS.add (PlanetObj (11,12, 13,14, 15));
	LS.add (FixedObj (10,20, 3));
	LS.add (FixedObj (40,50, 6));
	LS.add (AccelObj (-1,-2, 3,4, 5, 6));
	LS.addTail (PlanetObj (101,102, 103,104, 105));
	
	cout << "LS:\n"		<< LS		<< endl;
	cout << "LSmenbers:"	<< LS.nSize ()	<< endl;
	
	cout << "-------------------------------- 1 Step Time Evoluted" << endl;
	LS.step (0.001);
	
	cout << "LS:\n"		<< LS		<< endl;
	cout << "LSmenbers:"	<< LS.nSize ()	<< endl;
	
	cout << "-------------------------------- End " << endl;
}

void Test2 (void)
{
	List	LS;
	LS.add (PlanetObj (-5,-5,	 1,-1,	5,	2));
	LS.add (PlanetObj ( 5 ,5,	-1, 1,	5,	2));
	LS.add (PlanetObj ( 0, 0,	 2, 0,	0,	1));
		
	Real	dt = 0.001,
		rG = 10.0;
	int	nMax = 100000;
	
	PutPrologue ();
	for (int n = 0 ; n < nMax ; n++) {
		LS.step (dt, rG);
		
		if (n % 100 != 0)
			continue;
		
		LS.put ();
	}
	PutEpilogue ();
}

void Test3 (void) // constant accelerating satellite
{
	List	LS;
	LS.add (FixedObj  ( 0 ,0,			1,		5));
	LS.add (AccelObj  ( 1, 0,	 0,     1,	1,	0.01,	1));
	
	Real	dt = 0.0001,
		rG = 1.0;
	int	nMax = 100000;
	
	PutPrologue ();
	for (int n = 0 ; n < nMax ; n++) {
		LS.step (dt, rG);
		
		if (n % 100 != 0)
			continue;
		
		LS.put ();
	}
	PutEpilogue ();
}

#endif	// _bulk_h