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

void Psitions(double **Positions,double sigma2, int N, double size);
void Vtesses(double **Vitesses,int N);
void allocdouble2(double ***array,int  xsize,int ysize);
double distance2(double **Positions,int i, int j, double size);
void allocdouble(double **array,long  xsize);
void snapshot(double **Positionst,double **Vitesses, int flag,int N, int size, int un, int deux);
void calcul_temps(double **Positions,double **Vitesses,double *tempscol,int *col,double size,double sigma, 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;
}


// Definit la norme 2 pour les distances

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


/*Definit les positions*/

void Psitions(double **Positions,double sigma2, int N, double size) {
	
	int i,j,overlap;
	
	for (i=0; i<N; i++) {
		do {
			overlap=0;
			Positions[0][i]=size*drand48();
			Positions[1][i]=size*drand48();
			j=0;
			while(j<i && (overlap==0)) {
			  if (distance2(Positions,i,j,size)<sigma2) overlap=1;
				++j;
			}
		}
		while(overlap==1);
	}
}



/*Definit les vitesses (distribution gaussienne)*/

void Vtesses(double **Vitesses,int N) {
	int i;
	double Ec=0;
	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/(double)N*vitessecentregrave[0];
		Vitesses[1][i]=Vitesses[1][i]-1/(double)N*vitessecentregrave[1];
		Ec=Ec+0.5*(Vitesses[0][i]*Vitesses[0][i]+Vitesses[1][i]*Vitesses[1][i]);
	}
	for (i=0; i<N; i++) {
		Vitesses[0][i]=Vitesses[0][i]/sqrt(Ec);
		Vitesses[1][i]=Vitesses[1][i]/sqrt(Ec);
	}
}

//Permet de creer des images

void snapshot(double **Positionst, double **Vitesses,int flag,int N, int size, int un, int deux)
{
  /* Declare the image */
    gdImagePtr im;

    /* Declare output files */
    FILE *pngout;

    /* Declare color indexes */
    int black;
    int white;
	
	//size of the picture
	int pictsize=640;
    
    /* Allocate the image: 640 pixels across by 640 pixels tall */
    im = gdImageCreate(pictsize, pictsize);
    
    /* Allocate the color black (red, green and blue all minimum).
	Since this is the first color in a new image, it will
	be the background color. */
	white = gdImageColorAllocate(im, 255, 255, 255);
  
	/* Allocate the color white (red, green and blue all maximum). */
	black = gdImageColorAllocate(im, 0, 0, 0);
	int yellow = gdImageColorAllocate(im, 255, 255, 0);
	char nnom[30];
 
	sprintf(nnom,"image%d.png",flag);
	pngout = fopen(nnom, "wb");
	
	int Posx,Posy,Vx,Vy,rad;
	rad=pictsize/size;
	
	int i;
	
	for (i=0;i<N;i++)
		{
			Posx=rint(Positionst[0][i]*pictsize/size);
			Posy=rint(Positionst[1][i]*pictsize/size);
			if (i != un && i !=deux) {
				gdImageFilledArc(im, Posx, Posy, rad, rad, 0, 360, black,gdArc);
			}
			else {
				gdImageFilledArc(im, Posx, Posy, rad, rad, 0, 360, yellow,gdArc);
			}
			Vx=Posx+rint(Vitesses[0][i]*sqrt(N)*rad);
			Vy=Posy+rint(Vitesses[1][i]*sqrt(N)*rad);
			gdImageLine(im, Posx, Posy, Vx, Vy, black);
		}
	
	gdImagePng(im, pngout);

	
	fclose(pngout);

	
	gdImageDestroy(im);
	
}

//Sous programme pour calculer les temps de collision

void calcul_temps(double **Positionst,double **Vitesses,double *tempscol,int *col,double size,double sigma, int i, int j)
{
	
	double temps=1e8,a,b,c,Delta;
	a = pow(Vitesses[0][i]-Vitesses[0][j],2.0)+pow(Vitesses[1][i]-Vitesses[1][j],2.0);
	
	double Posimodij[2];
	
	Posimodij[0]=Positionst[0][i]-Positionst[0][j]-size*rint((Positionst[0][i]-Positionst[0][j])/size);
	Posimodij[1]=Positionst[1][i]-Positionst[1][j]-size*rint((Positionst[1][i]-Positionst[1][j])/size);
	
	
	b = 2*(Posimodij[0]*(Vitesses[0][i]-Vitesses[0][j])+Posimodij[1]*(Vitesses[1][i]-Vitesses[1][j]));
	
	c = pow(Posimodij[0],2.0)+pow(Posimodij[1],2.0)-sigma*sigma;
	
	Delta = b*b-4*a*c;

	if (Delta >= 0) {
		double t1,t2,ttmin,ttmax;
		t1=(-b-sqrt(Delta))/(2*a);
		
		t2=(-b+sqrt(Delta))/(2*a);
		if (t1<t2) {ttmin=t1;ttmax=t2;}else{ttmin=t2;ttmax=t1;}
		if (ttmin>0) {temps=ttmin;}
		else { if (ttmax>0) {temps=ttmax;}}
		
	}
	
	if (temps<tempscol[i]) {
		tempscol[i]=temps;
		col[i]=j;
	}
	if (temps<tempscol[j]) {
		tempscol[j]=temps;
		col[j]=i;
	}
}

/*###################################################################
 ###################################################################*/


int main(int argc, char *argv[]) {

	/*Definit les compteurs*/
	int i, j;
		
	double freq2[20];
	for (i=0;i<20;i++) {
		freq2[i]=0;
	}
	
	//Determine le nombre de fois ou la simulation est repetee pour realiser les moyennes
	int param,iter;
	for (param=0;param<iter;param++) {
	
	double **Positions, **Positionst;
	double **Vitesses;
	
	//initialise les temps
	
	double tmin=1e8;
	double size=20.0;

	/*Nombre de particules*/
	int N;
	
	/*Diametre de collision*/
	double sigma;
	sigma=1.0;
	
	/*Coefficient elastique*/
	double alpha;
		
	double t=0.0;
	int flag=0;
	//double tmax=1e6/(1-alpha*alpha);
	double tmax;

	double tmes;
	/* 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 eviter 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);
		fscanf(debut,"%d%*[^\n]", &iter);
		fscanf(debut,"%le%*[^\n]", &tmes);

			fscanf(debut,"%le%*[^\n]", &tmax);
		
	}


	allocdouble2(&Positions,2,N);
	allocdouble2(&Positionst,2,N);
	allocdouble2(&Vitesses,2,N);
			
	//Positions et vitesses pour t=0
	Psitions(Positions,sigma,N, size);
	Vtesses(Vitesses,N);
	
	//Initialise les positions dependants du temps
	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;
	}

	int col[N];
	for (i=0; i<N; i++) {
	  col[i]=-1;
	}

	// Calcul le temps de la premiere collision
	int un,deux;
	for (i=0; i<N; i++) {
		for (j=0; j<N; j++) {
			if (j != i) {
		  calcul_temps(Positions,Vitesses,tempscol,col,size,sigma,i,j);
			}
		}
	}

	for (i=0; i<N; i++) {
	  if (tempscol[i]<tmin) {
	    tmin=tempscol[i];
	    un=i;
	    deux=col[i];
		}
	}

		//Calcul de l'energie cinetique
		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]);
		}
	
		
	//	incremente le temps
	t=t+tmin;
	flag=flag+1;
	
	
	/*Definit les positions au temps tmin*/
	for (i=0; i<N; i++) {
	  Positionst[0][i]=Positionst[0][i]+Vitesses[0][i]*tmin-size*floor((Positionst[0][i]+Vitesses[0][i]*tmin)/size);
	  Positionst[1][i]=Positionst[1][i]+Vitesses[1][i]*tmin-size*floor((Positionst[1][i]+Vitesses[1][i]*tmin)/size);
	}

	double Posimodundeuxt[2];
	Posimodundeuxt[0]=Positionst[0][un]-Positionst[0][deux]-size*rint((Positionst[0][un]-Positionst[0][deux])/size);
	Posimodundeuxt[1]=Positionst[1][un]-Positionst[1][deux]-size*rint((Positionst[1][un]-Positionst[1][deux])/size);
	
	/*Vitesses apres le choc*/
	double DeltaVDeltaRsurR2;
	DeltaVDeltaRsurR2 = ((Vitesses[0][deux]-Vitesses[0][un])*Posimodundeuxt[0]+(Vitesses[1][deux]-Vitesses[1][un])*Posimodundeuxt[1])/(pow(Posimodundeuxt[0],2)+pow(Posimodundeuxt[1],2));
	
	Vitesses[0][un]=Vitesses[0][un]+(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[0];
	Vitesses[1][un]=Vitesses[1][un]+(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[1];
	
	Vitesses[0][deux]=Vitesses[0][deux]-(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[0];
	Vitesses[1][deux]=Vitesses[1][deux]-(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[1];
	
	

/*###################################################################
 ###################################################################*/

#ifdef ANAL2
	//Definit le temps de mesure de la distribution des vitesses
	double p=tmes;
#endif
		
	// Boucle avec le temps qui passe
	while (t<tmax) {
		
		//Met a jour les temps de collisions des particules dans la boite
		for (i=0; i<N; i++) {
			tempscol[i]=tempscol[i]-tmin;
			if (i==un || i==deux || tempscol[i]<=0) {
				tempscol[i]=1e8;
				col[i]=-1;
			}
		}
		
		tmin=1e8;

		//Calcul les collisions des particules qui devaient rencontrer une particule ayant deja collisionnee
		for (i=0; i<N; i++) {
		  if (i != un && i != deux) {
		    if (col[i]==un || col[i]==deux) {
		      tempscol[i]=1e8;
		      col[i]=-1;
		      for (j=0; j<N; j++) {
				  if (j != i) {
					  calcul_temps(Positionst,Vitesses,tempscol,col,size,sigma,i,j);
				  }
		      }
		    }
		  }
		}


		//Calcul les collisions des deux particules venant de collisionner
		for (i=0; i<N; i++) {
			if (i != deux && i !=un) {
			  calcul_temps(Positionst,Vitesses,tempscol,col,size,sigma,i,un);
			}
		  
			if (i != un && i != deux) {
			  calcul_temps(Positionst,Vitesses,tempscol,col,size,sigma,i,deux);
			}
		}
		
		for (i=0; i<N; i++) {
			if (tempscol[i]<tmin) {
				tmin=tempscol[i];
				un=i;
				deux=col[i];
			}
		}
		
		
		//Calcul l'energie cinetique
		Ec=0;
		for (i=0; i<N; i++) {
		  Ec=Ec+0.5*(Vitesses[0][i]*Vitesses[0][i]+Vitesses[1][i]*Vitesses[1][i]);
		}
		
#ifdef ANAL
	// Permet de suivre l'evolution de l'energie cinetique en fonction du temps
	printf("%f %f %f %f \n",t,Ec,log(t),log(Ec));
#endif

#ifdef ANAL2
	// Nous donne la distribution des vitesses lors d'une simulation
	if (t>p) {
		int freq[20];
		int ibis=0;
		for (ibis=0;ibis<20;ibis++) {
			freq[ibis]=0;
		}
		
		for (i=0;i<N;i++) {
			for (ibis=0;ibis<20;ibis++) {
				if (Vitesses[0][i]>(-0.2+ibis*0.02) && Vitesses[0][i]<=(-0.2+(ibis+1)*0.02)) {
					freq[ibis]=freq[ibis]+1;
				}
			}
		}
		for (ibis=0;ibis<20;ibis++) {
			freq2[ibis]=freq2[ibis]+freq[ibis];
		}
		p=tmax+1;
	}
#endif

		
#ifdef DEBUG
		for (i=0; i<N; i++) {
		    for (j=0; j<i; j++) {
		      if (distance2(Positionst,i,j,size)<(sigma-1e-8)) {
			printf("overlapa %i %i %e %e paire %e %e %d \n",i,j,t,distance2(Positionst,i,j,size)-sigma,Vitesses[0][i]-Vitesses[0][j],Vitesses[1][i]-Vitesses[1][j],flag);
			printf("********************\n");
			//exit(1);
		      }
		    }
		    //	  }
		}
#endif
		
		//Definit les positions au temps t
		for (i=0; i<N; i++) {
		  Positionst[0][i]=Positionst[0][i]+Vitesses[0][i]*tmin-size*floor((Positionst[0][i]+Vitesses[0][i]*tmin)/size);
		  Positionst[1][i]=Positionst[1][i]+Vitesses[1][i]*tmin-size*floor((Positionst[1][i]+Vitesses[1][i]*tmin)/size);
		}
		

#ifdef PICT
		// Cree une image a chaque choc ce qui peut permettre de realiser un film
		snapshot(Positionst,Vitesses,flag,N,size,un,deux);
#endif

		
#ifdef DEBUG
		for (i=0; i<N; i++) {
		    for (j=0; j<i; j++) {
		      if (distance2(Positionst,i,j,size)<(sigma-1e-8)) {
			printf("overlap %i %i %e %e paire %i %e %e  %d \n",i,j,t,distance2(Positionst,i,j,size)-sigma,col[i],tempscol[i],tmin,flag);
			printf("----------------------\n");
			//exit(1);
		      }
		    }
		    //	  }
		}
#endif
		

		/*Vitesses apres le choc*/
		double Posimodundeuxt[2];
		Posimodundeuxt[0]=Positionst[0][un]-Positionst[0][deux]-size*rint((Positionst[0][un]-Positionst[0][deux])/size);
		Posimodundeuxt[1]=Positionst[1][un]-Positionst[1][deux]-size*rint((Positionst[1][un]-Positionst[1][deux])/size);
		
		double DeltaVDeltaRsurR2;
		DeltaVDeltaRsurR2 = ((Vitesses[0][deux]-Vitesses[0][un])*Posimodundeuxt[0]+(Vitesses[1][deux]-Vitesses[1][un])*Posimodundeuxt[1])/(pow(Posimodundeuxt[0],2)+pow(Posimodundeuxt[1],2));
		
		Vitesses[0][un]=Vitesses[0][un]+(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[0];
		Vitesses[1][un]=Vitesses[1][un]+(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[1];
		
		Vitesses[0][deux]=Vitesses[0][deux]-(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[0];
		Vitesses[1][deux]=Vitesses[1][deux]-(1+alpha)/2.0*DeltaVDeltaRsurR2*Posimodundeuxt[1];
		

		//	incremente le temps
		t=t+tmin;
		flag=flag+1;
		
	
	}
	}
	
#ifdef ANAL2
	//Moyenne la distribution des vitesses sur l'ensemble des simulations
	for (i=0;i<20;i++) {
		freq2[i]=freq2[i]/(param+1);
		printf("%f %f \n",freq2[i],-0.2+(i+0.5)*0.02);
	}
#endif
	
	return (0);
}