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

void Psitions(double **Positions,double sigma2, int N);
void Vtesses(double **Vitesses,int N);
void allocdouble2(double ***array,int  xsize,int ysize);
void allocdouble(double **array,long  xsize);
double distance2(double **Positions,int i, int j);


void allocdouble2(double ***array,int  xsize,int ysize)

/* dynamically allocates a two dimensional (*array)[xize][ysize]
 with datasize bytes per element */
{ int i;
	*array=(double **) malloc(xsize*sizeof(double *));
	/*  if (*array ==  NULL) */
	/*   {fprintf(stderr," memory allocation problems\n"); */
	/*    exit(-1);} */
	for(i=0;i<xsize;i++){
		(*array)[i]=(double *) calloc(ysize,sizeof(double));
	}
	return;
}

void allocdouble(double **array,long  xsize)
/* dynamically allocates a one dimensional (*array)[xize]
 with datasize bytes per element */
{
	/*  if (*array != NULL)  /\* free previous allocation *\/ */
	/*    free(*array); */
	/* allocate memory for array of xsize pointers plus data */
	*array=(double *) calloc(xsize, sizeof(double));
	if (*array ==  NULL)
		{fprintf(stderr," memory allocation problems\n");
		exit(-1);}
	return;
}


// Définit la norme 2 pour les distances

double distance2(double **Positions,int i, int j)
{
	return(
		   (Positions[0][i]-Positions[0][j])*(Positions[0][i]-Positions[0][j])
		   +(Positions[1][i]-Positions[1][j])*(Positions[1][i]-Positions[1][j]));
}


/*Définit les positions*/

void Psitions(double **Positions,double sigma2, int N) {
	
	int i,j,overlap,retry;
	
	for (i=0; i<N; i++) {
		overlap=0; retry=0;
		do {
			Positions[0][i]=100*drand48();
			Positions[1][i]=100*drand48();
			/* 	printf("i%d \n",i); */
			j=0;
			overlap=0;
			while(j<i && (overlap==0)) {
				
				if (distance2(Positions,i,j)<sigma2) {
					overlap=1;
					
				}
				j++;
			}
			if (overlap==1) {
				retry=1;
			}
			else {
				retry=0;
			}
		}
		while(retry==1);
	}
}



/*Définit les vitesses (distribution gaussienne)*/

void Vtesses(double **Vitesses,int N) {
	int i;
	double vitessecentregrave[2];
	vitessecentregrave[0]=0;
	vitessecentregrave[1]=0;
	for (i=0; i<N; i++) {
		double u,v;
		u=drand48();
		v=drand48();
		Vitesses[0][i]=sqrt(-2*log(u))*cos(2*M_PI*v);
		Vitesses[1][i]=sqrt(-2*log(v))*cos(2*M_PI*u);
		vitessecentregrave[0]=vitessecentregrave[0]+Vitesses[0][i];
		vitessecentregrave[1]=vitessecentregrave[1]+Vitesses[1][i];
	}
	for (i=0; i<N; i++) {
		Vitesses[0][i]=Vitesses[0][i]-1/N*vitessecentregrave[0];
		Vitesses[1][i]=Vitesses[1][i]-1/N*vitessecentregrave[1];
	}
}


int main(int argc, char *argv[]) {
	
	
	double **Positions,**Positionst;
	double **Vitesses;
	
	//initialise les temps
	
	double tmin=1e8;
	
	double t=0.0;
	int N;
	double alpha;
	/* quelques lignes pour se rappeler comment ca marche */
	if (argc != 2)
		{
			printf("usage: a2  <fichier de donnees>\n");
			exit (1);
		}
	
	/*Nombre de particules  et alpha sont les parametres qui varient entre deux simulations*/
	/*  on les mets dans un fichier pour éviter de recompiler le programme*/
	
	{
		FILE *debut;
		if ((debut = fopen (argv[1],"r")) == NULL)
			{ fprintf(stderr,"Probleme d'ouverture de %s\n",argv[1]);
				exit (1);
			}
		fscanf(debut,"%d%*[^\n]", &N);
		fscanf(debut,"%le%*[^\n]", &alpha);
		
	}
	
	
	
	
	
	allocdouble2(&Positions,2,N);
	allocdouble2(&Positionst,2,N);
	allocdouble2(&Vitesses,2,N);
	
	/*Diamètre de collision*/
	const double sigma=1;
	
	
	
	
	
	/* on stoppe la simulation quand on quitte le free cooling state */
	double tmax=1/(1-alpha*alpha);
	
	/*Pour changer de nombre aléatoire à chaque fois*/
	//srand48(time(NULL));
	
	
	//Positions et vitesses pour t=0
	
	printf("initialisation!\n");
	Psitions(Positions,sigma,N);
	Vtesses(Vitesses,N);
	
	
	printf("on démarre!\n");
	
	/*Définit les compteurs*/
	int i, j;  
	for (i=0; i<N; i++) {
		Positionst[0][i]=Positions[0][i];
		Positionst[1][i]=Positions[1][i];
	}		
	
	double tempscol[N];
	for (i=0; i<N; i++){
		tempscol[i]=1e8;
	}
	
	double col[N];
	for (i=0; i<N; i++) {
		col[i]=-1;
	}
	double Ec=0;
	
	for (i=0; i<N; i++) {
		Ec=Ec+0.5*(Vitesses[0][i]*Vitesses[0][i]+Vitesses[1][i]*Vitesses[1][i]);
	}
	
	
	// Calcule le temps de la premiere collision
	
	
	double Delta;
	int un,deux;
	double temps;
	
	
	for (i=0; i<N; i++) {
		for (j=0; j<N; j++) {
			// attention il faut écarter les termes self       
			if (i != j) {
				temps=1e8;
				double a = powf(Vitesses[0][i]-Vitesses[0][j],2)+powf(Vitesses[1][i]-Vitesses[1][j],2);
				double b = 2*((Positions[0][i]-Positions[0][j])*(Vitesses[0][i]-Vitesses[0][j])+(Positions[1][i]-Positions[1][j])*(Vitesses[1][i]-Vitesses[1][j]));
				double c = powf(Positions[0][i]-Positions[0][j],2)+powf(Positions[1][i]-Positions[1][j],2)-sigma*sigma;
				
				double Delta = b*b-4*a*c;
				
				if (Delta >= 0) {
					if (a>0) {
						temps=(-b-sqrt(Delta))/(2*a);
					}
					else {
						temps=(-b+sqrt(Delta))/(2*a);
					}
				}
				
				if (temps<tempscol[i] && temps>0) {
					tempscol[i]=temps;
					col[i]=j;
				}
				
				if (temps<tempscol[j] && temps>0) {
					tempscol[j]=temps;
					col[j]=i;
				}
				
			}
		}
	}
	
	for (i=0; i<N; i++) {
		if (tempscol[i]<tmin) {
			tmin=tempscol[i];
			un=i;
			deux=col[i];
		}
	}
	
	
	
	
	//	incrémente le temps
	
	t+=tmin;
	
	printf("%f %f \n",t,Ec);	
	
	
	/*Définit les positions au temps tmin*/
	
	
	
	for (i=0; i<N; i++) {
		Positionst[0][i]=fmod(Positionst[0][i]+Vitesses[0][i]*tmin,100);
		Positionst[1][i]=fmod(Positionst[1][i]+Vitesses[1][i]*tmin,100);
	}
	
	
	/*Vitesses après le choc*/
	
	double DeltaVDeltaRsurR2;
	DeltaVDeltaRsurR2 = ((Vitesses[0][deux]-Vitesses[0][un])*(Positionst[0][un]-Positionst[0][deux])+(Vitesses[1][deux]-Vitesses[1][un])*(Positionst[1][un]-Positionst[1][deux]))/(powf(Positionst[0][un]-Positionst[0][deux],2)+powf(Positionst[1][un]-Positionst[1][deux],2));
	
	Vitesses[0][un]=Vitesses[0][un]+(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[0][un]-Positionst[0][deux]);
	
	Vitesses[1][un]=Vitesses[1][un]+(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[1][un]-Positionst[1][deux]);
	
	Vitesses[0][deux]=Vitesses[0][deux]-(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[0][un]-Positionst[0][deux]);
	
	Vitesses[1][deux]=Vitesses[1][deux]-(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[1][un]-Positionst[1][deux]);
	
	
	
	
	while (t<tmax) {
		
		for (i=0; i<N; i++) {
			tempscol[i]=tempscol[i]-tmin;
			if (tempscol[i]<=0) {
				tempscol[i]=1e8;
			}
		}	
		
		tmin=1e8;
		
		
		
		for (i=0; i<N; i++) {
			if (i != deux) {
				temps=1e8;
				
				double a = powf(Vitesses[0][i]-Vitesses[0][un],2)+powf(Vitesses[1][i]-Vitesses[1][un],2);
				double b = 2*((Positionst[0][i]-Positionst[0][un])*(Vitesses[0][i]-Vitesses[0][un])+(Positionst[1][i]-Positionst[1][un])*(Vitesses[1][i]-Vitesses[1][un]));
				double c = powf(Positionst[0][i]-Positionst[0][un],2)+powf(Positionst[1][i]-Positionst[1][un],2)-sigma*sigma;
				
				Delta = b*b-4*a*c;
				
				if (Delta >= 0) {
					if (a>0) temps=(-b-sqrt(Delta))/(2*a); else temps=(-b+sqrt(Delta))/(2*a);
				}
				
				if (temps<tempscol[i] && temps>0) {
					tempscol[i]=temps;
					col[i]=un;
				}
				
				if (temps<tempscol[un] && temps>0) {
					tempscol[un]=temps;
					col[un]=i;
				}
			}
			
			if (i != un) {
				
				temps=1e8;
				
				double a = powf(Vitesses[0][i]-Vitesses[0][deux],2)+powf(Vitesses[1][i]-Vitesses[1][deux],2);
				double b = 2*((Positionst[0][i]-Positionst[0][deux])*(Vitesses[0][i]-Vitesses[0][deux])+(Positionst[1][i]-Positionst[1][deux])*(Vitesses[1][i]-Vitesses[1][deux]));
				double c = powf(Positionst[0][i]-Positionst[0][deux],2)+powf(Positionst[1][i]-Positionst[1][deux],2)-sigma*sigma;
				
				Delta = b*b-4*a*c;
				
				if (Delta >= 0) {
					if (a>0) temps=(-b-sqrt(Delta))/(2*a); else temps=(-b+sqrt(Delta))/(2*a);
				}
				
				if (temps<tempscol[i] && temps>0) {
					tempscol[i]=temps;
					col[i]=un;
				}
				
				if (temps<tempscol[deux] && temps>0) {
					tempscol[deux]=temps;
					col[deux]=i;
				}
			}
			
		}
		
		
		for (i=0; i<N; i++) {
			if (tempscol[i]<tmin) {
				tmin=tempscol[i];
				un=i;
				deux=col[i];
			}
		}
		
		
		//calcule l'énergie cinétique
		Ec=0;
		for (i=0; i<N; i++) {
			Ec=Ec+0.5*(Vitesses[0][i]*Vitesses[0][i]+Vitesses[1][i]*Vitesses[1][i]);
		}
		
		
		
		printf("%f %f \n",t,Ec);
		//	printf("\n ---- %f ---- \n",tmin);
		
		//Définit les positions au temps t
		
		for (i=0; i<N; i++) {
			Positionst[0][i]=fmod(Positionst[0][i]+Vitesses[0][i]*tmin,100);
			Positionst[1][i]=fmod(Positionst[1][i]+Vitesses[1][i]*tmin,100);
		}
		
		
		//Vitesses après le choc
		
		
		DeltaVDeltaRsurR2 = ((Vitesses[0][deux]-Vitesses[0][un])*(Positionst[0][un]-Positionst[0][deux])+(Vitesses[1][deux]-Vitesses[1][un])*(Positionst[1][un]-Positionst[1][deux]))/(powf(Positionst[0][un]-Positionst[0][deux],2)+powf(Positionst[1][un]-Positionst[1][deux],2));
		
		Vitesses[0][un]=Vitesses[0][un]+(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[0][un]-Positionst[0][deux]);
		
		Vitesses[1][un]=Vitesses[1][un]+(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[1][un]-Positionst[1][deux]);
		
		Vitesses[0][deux]=Vitesses[0][deux]-(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[0][un]-Positionst[0][deux]);
		
		Vitesses[1][deux]=Vitesses[1][deux]-(1+alpha)/2*DeltaVDeltaRsurR2*(Positionst[1][un]-Positionst[1][deux]);
		
		//	incrémente le temps
		t=t+tmin;
		
		
	}
	
	
	return (0);
}