#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);
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;
}

// 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();
			j=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*/

void Vtesses(double **Vitesses,int N) {
		int i;
		for (i=0; i<N; i++) {
			Vitesses[0][i]=10*drand48();
			Vitesses[1][i]=10*drand48();
		}
}

// Calcule le temps de collision
/*
void calc_colis(double **Vitesses, double **Positions, int N, double sigma, double tmin) {
	double a, b, c;
	double Delta;
	int un,deux;
	int i,j;
 
	for (i=0; i<N; i++) {
		for (j=0; j<N; j++) {
			
 
			a = powf(Vitesses[0][i]-Vitesses[0][j],2)+powf(Vitesses[1][i]-Vitesses[1][j],2);
			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]));
			c = powf(Positions[0][i]-Positions[0][j],2)+powf(Positions[1][i]-Positions[1][j],2)-sigma*sigma;
 
			Delta = b*b-4*a*c;
 // printf("%lf |", Delta);
			double temps=powf(10,8);
			if (Delta >= 0) {
				if (a>0) temps=(-b-sqrt(Delta))/(2*a); else temps=(-b+sqrt(Delta))/(2*a);
			}
 
			if (temps<tmin && temps>0) {
				tmin=temps;
				un=i;
				deux=j;
			}
			printf("%f \n",tmin);
		}
	}
}
*/

int main() {
	double **Positions;
	double **Vitesses;
	
	double tmin;
	
	/*Nombre de particules*/
	int N;
	N=100;

	allocdouble2(&Positions,2,N);
	allocdouble2(&Vitesses,2,N);
 
	printf("totototo\n");

  
	
	/*Diamètre de collision*/
	double sigma;
	sigma=1.0;
	
	/*Coefficient élastique*/
	double epsilon;
	epsilon=1.0;
	
	
	/*Pour changer de nombre aléatoire à chaque fois*/
	//srand48(time(NULL));
	
	/*Définit les compteurs*/
	
	Psitions(Positions,sigma,N);
	
	Vtesses(Vitesses, N);
	
	int i, j;
	
	/*
	for (i=0; i<2; i++) {
		for (j=0; j<N; j++) {
			printf("%lf |", Positions[i][j]);
		}
		printf("\n");
	}
	
	for (i=0; i<2; i++) {
		for (j=0; j<N; j++) {
			printf("%lf |", Vitesses[i][j]);
		}
		printf("\n");
	}
	*/


	
	//calc_colis(Vitesses,Positions,N,sigma,tmin);
	double t=0.0;
	double tmax=10.0;
	
	double Positionst[2][N];
	
	for (i=0; i<N; i++) {
		Positionst[0][i]=Positions[0][i];
		Positionst[1][i]=Positions[1][i];
	}		
	
	
	while (t<tmax) {
		double Ec=0;
		tmin=powf(10,8);
		
	// Calcule le temps de collision
	
	double a, b, c;
	double Delta;
	int un,deux;
	
	
	for (i=0; i<N; i++) {
		for (j=0; j<N; j++) {
			
			a = powf(Vitesses[0][i]-Vitesses[0][j],2)+powf(Vitesses[1][i]-Vitesses[1][j],2);
			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]));
			c = powf(Positions[0][i]-Positions[0][j],2)+powf(Positions[1][i]-Positions[1][j],2)-sigma*sigma;
			
			Delta = b*b-4*a*c;
			// printf("%lf |", Delta);
			double temps=powf(10,8);
			if (Delta >= 0) {
				if (a>0) temps=(-b-sqrt(Delta))/(2*a); else temps=(-b+sqrt(Delta))/(2*a);
			}
			
			if (temps<tmin && temps>0) {
				tmin=temps;
				un=i;
				deux=j;
			}
		}
	}
		
		//calcule l'énergie cinétique
		
		for (i=0; i<N; i++) {
			Ec=Ec+0.5*(Vitesses[0][i]*Vitesses[0][i]+Vitesses[1][i]*Vitesses[1][i]);
		}
		
	//	incrémente le temps
	
		t=t+tmin;
	
		printf("%f | %f \n",t,Ec);

	
	/*Définit les positions au temps t*/
	
		
		
	for (i=0; i<N; i++) {
		Positionst[0][i]=Positionst[0][i]+Vitesses[0][i]*tmin;
		Positionst[1][i]=Positionst[1][i]+Vitesses[1][i]*tmin;
	}
	
	
	/*Vitesses après le choc*/
	
	Vitesses[0][un]=(Vitesses[0][un]+epsilon*((Vitesses[0][deux]-Vitesses[0][un])*(Positions[0][un]-Positions[0][deux])+(Vitesses[1][deux]-Vitesses[1][un])*(Positions[1][un]-Positions[1][deux]))/(powf(Positions[0][un]-Positions[0][deux],2)+powf(Positions[1][un]-Positions[1][deux],2))*(Positions[0][un]-Positions[0][deux]));
	
	Vitesses[1][un]=(Vitesses[1][un]+epsilon*((Vitesses[0][deux]-Vitesses[0][un])*(Positions[0][un]-Positions[0][deux])+(Vitesses[1][deux]-Vitesses[1][un])*(Positions[1][un]-Positions[1][deux]))/(powf(Positions[0][un]-Positions[0][deux],2)+powf(Positions[1][un]-Positions[1][deux],2))*(Positions[1][un]-Positions[1][deux]));
	
	Vitesses[0][deux]=(Vitesses[0][deux]+epsilon*((Vitesses[0][un]-Vitesses[0][deux])*(Positions[0][un]-Positions[0][deux])+(Vitesses[1][un]-Vitesses[1][deux])*(Positions[1][un]-Positions[1][deux]))/(powf(Positions[0][un]-Positions[0][deux],2)+powf(Positions[1][un]-Positions[1][deux],2))*(Positions[0][un]-Positions[0][deux]));
	
	Vitesses[1][deux]=(Vitesses[1][deux]+epsilon*((Vitesses[0][un]-Vitesses[0][deux])*(Positions[0][un]-Positions[0][deux])+(Vitesses[1][un]-Vitesses[1][deux])*(Positions[1][un]-Positions[1][deux]))/(powf(Positions[0][un]-Positions[0][deux],2)+powf(Positions[1][un]-Positions[1][deux],2))*(Positions[0][un]-Positions[0][deux]));

	}

		
	return (0);
}