/* movie.c		inertia, gravity, bouncing */

#include "rtd.h"
#include <math.h>

double dot();
float tx43;
#define cube(a) ((tx43=a)*(tx43)*(tx43))
#define density 1.0

movie(bl, nob, step, Step, gravity, bounce, attraction)
struct ball *bl;
int nob, step, Step, bounce;
float gravity, attraction;
{	register int i, j;
	register struct ball *bli, *blj;
	struct vector f;
	for (step = 0; step<Step; step++)
	{	for (i=0, bli=bl; i<nob; i++, bli++)
		{	if (bounce) for (j=i+1, blj=bli+1; j<nob; j++, blj++)
			{	register float 
				      massi, massj, /* mass of balls (i&j) */
				      fi, fj,	    /* mass fraction (i&j) */
				      tmp;
				struct vector cmv,  /* center mass velocity*/
				       rvi, rvj,    /* relative vel (i&j)  */
				       d;	    /* distance between i&j*/
				sv(&d, &bli->s.cent, &blj->s.cent);
				vecl(&d);
				if (d.l <= bli->s.rad + blj->s.rad)
				{	tmp = bli->s.rad; massi = tmp*tmp*tmp;
					tmp = blj->s.rad; massj = tmp*tmp*tmp;
					fi = massi/(massi+massj);
					fj = massj/(massi+massj);
					cmv.x = fi*bli->v.x + fj*blj->v.x;
					cmv.y = fi*bli->v.y + fj*blj->v.y;
					cmv.z = fi*bli->v.z + fj*blj->v.z;

					  /*    rvi = bli->v - cmv      */
					sv(&rvi, &bli->v, &cmv);
					tmp = dot(&d,&rvi);
					  /*    d = 2*|d.rvi|/|d.d|*d   */	
					scamult(2.0*tmp/(d.l*d.l),&d);
					  /*    bli->v = d - rvi + cmv  */	
					sv(&rvi, &rvi,&d);
					av(&bli->v, &cmv, &rvi);

					vecl(&d);
					sv(&rvj, &blj->v, &cmv);
					tmp = dot(&d,&rvj);
					scamult(2.0*tmp/(d.l*d.l),&d);
					sv(&rvj, &rvj,&d);
					av(&blj->v, &cmv, &rvj);
			}	}
			if (attraction != 0.0) 
			{	register float dl, tmp;
				struct vector d;
				f.x = f.y = f.z = 0.0;
				for (j=0, blj=bl; j<nob; j++, blj++) if (i != j)
				{
					sv(&d, &bli->s.cent, &blj->s.cent);
					dl = sqrt(d.x*d.x + d.y*d.y + d.z*d.z);
					tmp = blj->s.rad/dl; tmp=tmp*tmp*tmp;
					scamult(tmp, &d); /* mass/(r*r) */
					av(&f, &f, &d);
				}
				scamult(attraction/(density*Step), &f);
				sv(&bli->v,  &bli->v, &f);
		}	}
		for (i=0, bli=bl; i<nob; i++, bli++)
		{	bli->s.cent.x += bli->v.x/Step;

			bli->v.y -= gravity;
			bli->s.cent.y += bli->v.y/Step;
			if (bli->s.cent.y < bli->s.rad && bli->v.y < 0)
				bli->v.y = - bli->v.y;

			bli->s.cent.z += bli->v.z/Step;
			if (bli->s.cent.z < bli->s.rad && bli->v.z < 0)
				bli->v.z = - bli->v.z;
}	}	}