#include <math.h>
#include <stdlib.h>
#include <stdio.h>

#define N 20								// Nombre de grains, pas d'alloc. dynamique

extern void gnuplot( int, int, char c[50]);	// fct externe qui sert a tracer des graphes avec gnuplot


int main(void){

double  dx, d=0.1, e=.5;			// e = coeff de restitution, elastique si e=1
						// d = rayon des grains
double tps, time=0, tmin, t0=1E8;		// t0 grand, pour la recherche de collision suivante
double V[N+1], X[N+1], Vmoy, vi, vj;		// liste de positions et vitesse des grains
int n, i, j, iold=0, bcl, bclmax=5000;		// bclmax, inutile (il vaudrait mieux un while)
char nom[50];					// chaine de caracteres pour fichier de sortie
FILE *fout, *fVit;				// fichiers trajectoires & vitesses finales


// init
for(n=0;n<=N;n++)
{V[n]=0; X[n]=n;}				// un grain "tous les 1" avec V=0
V[1]=1.3;					// 1er grain
V[N]=-1.;					// dernier grain

// on écrit les positions initiales dans le fichier fout
sprintf(nom,"Trajec.xls");			// nom du fichier des trajectoires
fout=fopen(nom,"w");				// on ouvre le fichier
fprintf(fout, "%f\t",time);			// 1ere colonne : time
for(n=1;n<=N;n++)				
fprintf(fout, "\t%.2f", X[n]);			// on ecrit les positions a t=0
fprintf(fout, "\n");				// on ferme le fichier
//******************


for(bcl=0;bcl<=bclmax;bcl++)	// boucle principale (sur les collisions)
{

tmin=t0;							// attention : t remis a zero a chaque coll.
for(n=1;n<=N-1;n++)						// on calcul pour chaque grain n l'instant de l'impact avec n+1
{
	dx = X[n]-X[n+1]+2*d;					// distance entre les 2 grains
	if(V[n]!=V[n+1])	{ tps = dx/(V[n+1]-V[n]);}	// instant de l'impact entre i et j
	else				{ tps =t0;}		// par convention t0 si jms de contact

	// on cherche la prochaine collision
	if(tps>0 && tps<tmin && fabs(dx)>1E-11)	// il faut que tps>0 et dx>... pour le bruit numerique
	{
		tmin=tps;					// la prochaine collision aura lieu a t=tim
		i=n;						// entre les grains i
		j=n+1;						// et j
	}	
}	// end for n


if(tmin==t0)							// si on n'a pas trouve de collision
{
	printf("\n\n fini : t= %f \n\n", time);			// affiche un message à l'écran
	fVit=fopen("Vit.xls","w");				// on va ecrire les vitesses finales dans un fichier
	for(n=1;n<=N;n++)					// ...
	fprintf(fVit, "%d\t%f\n", n, V[n]);			// ...
	fclose(fVit);

	gnuplot(N, (int) time +1, nom);				// on ecrit un fichier bash pour gnuplot
	system("bash 'plots.bash'\n");				// lance le script bash

	break;							// on arrete le prgm maintenant qu'il n'y a plus de collision
}

time+=tmin;						// le temps mesuré depuis le début du code
printf("coll entre %d et %d, \t dans tmin=%f\t a time=%f\n", i,j,tmin, time);
// On connait maintenant le tps de la prochaine collision : tmin entre i et j

// tout le monde bouge
for(n=1;n<=N;n++) X[n]+=V[n]*tmin;

// Nouvelles vitesses 
Vmoy=0.5* (V[i]+V[j]);			// dans le referentiel barycentrique
vi = V[i] - Vmoy;		V[i] = -e*vi + Vmoy;
vj = V[j] - Vmoy;		V[j] = -e*vj + Vmoy;

fprintf(fout, "%f", time);		// fichier trajectoires
for(n=1;n<=N;n++)
fprintf(fout, "\t%.2f", X[n]);
fprintf(fout, "\n");

}	// fin bcl
fclose(fout);

}	// fin prgm



// On ecrit un script bash pour gnuplot
// commande : bash 'plots.bash'
void gnuplot(int a, int b, char c[50])
{
FILE *fgnu;
int nn;

fgnu = fopen("plots.bash", "w");		// ouverture du fichier

fprintf(fgnu, "#!/bin/bash\n");
fprintf(fgnu, "gnuplot <<FIN\n");
fprintf(fgnu, "set terminal postscript eps\n");
fprintf(fgnu, "set output 'Traj.eps'\n");
fprintf(fgnu, "set multiplot\n");

for(nn=1;nn<=a;nn++)
fprintf(fgnu, "plot [0:%d] [0:%d] '%s' u 1:%d w l\n", b, a, c, nn+1);
//		    rangex rangey  nom	    with line
fprintf(fgnu, "FIN\n");
fprintf(fgnu, "gv Traj.eps &\n");
fclose(fgnu);

return;
}