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

/*Définit les positions*/

double distance2(double **Positions,int i, int j);
void Psitions(double **Positions,double sigma2, int N) {
	
  int i,j,overlap,retry;
  
  for (i=0; i<N; i++) {
    overlap=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;
	}
      if (overlap==1) {retry=1;} else {retry=0;}
      
    } while(retry==1);
    
    
  }
}


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

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

/* dynamically allocates a two dimensional (*array)[xize][ysize]
   with datasize bytes per element */
{ int i;
/*  if (*array != NULL)  /\* free previous allocation *\/ */
/*    free(*array); */
  /* allocate memory for array of ysize pointers plus data */
 *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 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();
  }
}
	
	

int main() {
  double **Positions;
  double **Vitesses;

	
	/*Nombre de particules*/
	int N;
	N=10;

	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*/
	int i, j, k;
	
	Psitions(Positions,sigma,N);

	/*Imprime a l'ecran les positions et les vitesses*/
		
	for (k=0; k<N; k++) {
		printf("%lf ", Positions[0][k]);
	}
	printf("\n");
	for (k=0; k<N; k++) {
		printf("%lf ", Positions[1][k]);
	}
	printf("\n ---------------- \n");

		
	for (k=0; k<N; k++) {
		printf("%lf ", Vitesses[0][k]);
	}
	printf("\n");
	for (k=0; k<N; k++) {
		printf("%lf ", Vitesses[1][k]);
	}
	
	printf("\n ---------------- \n");
	
	float t;
	t=0.0;
	

	
	/*Résoudre l'équation de degré 2 pour les collisions et trouver le temps minimal*/
	float a[N][N], b[N][N], c[N][N];
	float Delta[N][N];
	float temps[N][N][2];
	float tmin=powf(10,8);
	
	int un, deux;
	
	for (i=0; i<N; i++) {
		for (j=0; j<=i; j++) {

			a[i][j] = powf(Vitesses[0][i]-Vitesses[0][j],2)+powf(Vitesses[1][i]-Vitesses[1][j],2);
			b[i][j] = 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[i][j] = powf(Positions[0][i]-Positions[0][j],2)+powf(Positions[1][i]-Positions[1][j],2)-powf(sigma,2);
			
			Delta[i][j] = powf(b[i][j],2)-4*a[i][j]*c[i][j];
			printf("%f |", Delta[i][j]);
			
			if (Delta[i][j] >= 0) {
				temps[i][j][0]=(-b[i][j]+sqrt(Delta[i][j]))/(2*a[i][j]);
				temps[i][j][1]=(-b[i][j]+sqrt(Delta[i][j]))/(2*a[i][j]);
				
				for (k=0; k<=1; k++) {
					if (temps[i][j][k]>0 && temps[i][j][k] <= tmin) {
						tmin=temps[i][j][k];
						un=i;
						deux=j;
					}
				}
			}
		}
		
			printf("\n  %f  \n",tmin);
	}
	
	t=t+tmin;
	
	/*Définit les positions au temps t*/
	
	float Positionst[2][N];
	for (i=0; i<N; i++) {
		Positionst[0][i]=Positions[0][i]+Vitesses[0][i]*t;
	}
	
	for (i=0; i<N; i++) {
		Positionst[1][i]=Positions[1][i]+Vitesses[1][i]*t;
	}
	
	/*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]);

	
	printf("\n ---------------- \n");

	
	int y;
	for (y=0; y<N; y++) {
		printf("%f ", Positionst[0][y]);
	}
	
	printf("\n");
	for (y=0; y<N; y++) {
		printf("%f ", Positionst[1][y]);
	}
	
	
	printf("\n ---------------- \n");
	

		
	return (0);
}