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

// c++ TColl.cc ; ./a.out
// Writes in the file res.dat :
// first line : N,theta,dt,eps,deltaplot,thresh
// other lines : time log_{10}(Z[2]) log_{10}(Z[3]) ... log_{10}(Z[N])
// where Z[i]=min_{k=1,...,N} R_{k,i}
// where R_{k,i}=||X_k-Y||^2, where Y is the (i-1)-th neighbor of X_i,
// i.e. if N=4 and ||X_2-X_3||<||X_2-X_1||<||X_2-X_4||, then 
// R_{2,2}=||X_2-X_3||^2, R_{2,3}=||X_2-X_1||^2 and R_{2,4}=||X_2-X_4||^2.

// In gnuplot :
// set title "aaa"; set key outside; set key box
// plot for [i=2:N] 'res.dat' using 1:(column(i)+2*i) w l title sprintf("%d",i)

const int N=23; // Number of particles
const double T=1; // Maximum final time
const double theta=2.7; // Parameter of interaction
const double dt=1e-9; // Time step (OK with 1e-9 or 1e-10)
const double eps=1e-8; // Smoothing parameter (X/|X|^2 replaced by X/(|X|^2+eps) (OK with 1e-8)
const long deltaplot=1/dt/1e6; // Saves the data each "deltaplot" time steps
const double thresh=1e-6; // Stops at time T or when Z[N]<thresh (OK with 1e-5)
const double pi=3.14159265359;

double Normale()
{
  double u = ((double) rand()+1)/((double) RAND_MAX+1);
  double v = ((double) rand())/((double) RAND_MAX);
  double w = sqrt(-2.0*log(u))*cos(2*pi*v);
  return (w);
}

int comp (const void * a, const void * b)
{
    if (*(double*)a > *(double*)b) return 1;
    else if (*(double*)a < *(double*)b) return -1;
    else return 0;
}

int main()
{
  double ect,dtdrift,temps=0,tac,trucx,trucy;
  double X[N],Y[N],XX[N],YY[N],A[N][N],Z[N+1];
  long i=0;
  int j,k=0,l,m,kz,k1,k2;
  bool boo=0;
  FILE* pf;
  
  pf=fopen("res.dat","w");  
  srand((unsigned)(clock()+time(NULL)));  
  fprintf(pf,"N=%d   theta=%f   dt=%e   eps=%e   deltaplot=%ld   thresh=%f \n", N,theta,dt,eps,deltaplot,thresh);
  ect=sqrt(dt);
  dtdrift=dt*theta/(double)N; 
  kz=ceil(2*N/theta);
  k1=kz-1;
  k2=kz-2;  
  if ((k2-1)*(2-k2*theta/(double)N)>2) {k2=kz-1;};
  printf("N=%d, theta=%f, dt=%e, eps=%e, deltaplot=%ld. We have k0=%d, k1=%d, k2=%d. \n", N,theta,dt,eps,deltaplot,kz,k1,k2);

  for (j=0;j<N;j++)
    {
       X[j]=0.2*Normale();
       Y[j]=0.2*Normale();
    };

  do
    { 
      temps=temps+dt;
      i=i+1;
      for (j=0;j<N;j++)
        {
          XX[j]=X[j]+ect*Normale();
          YY[j]=Y[j]+ect*Normale();
        };
      for (j=0;j<N;j++)
        {
          trucx=0;
          trucy=0;
          for (l=0;l<N;l++)
            { 
              tac= (XX[j]-XX[l])*(XX[j]-XX[l])+(YY[j]-YY[l])*(YY[j]-YY[l])+eps;
              trucx=trucx+(XX[j]-XX[l])*dtdrift/tac;
              trucy=trucy+(YY[j]-YY[l])*dtdrift/tac;
            };
          X[j]=XX[j]-trucx;
          Y[j]=YY[j]-trucy;
        };  
      if (i>deltaplot)
        {
          i=0;
          for (j=0;j<N;j++)
            {
              for (l=0;l<N;l++)
               { 
                 A[j][l]=(X[j]-X[l])*(X[j]-X[l])+(Y[j]-Y[l])*(Y[j]-Y[l]);
               };
            };  
          for (j=0;j<N;j++)
            {
              qsort(A[j],N,sizeof(double),comp);
            };
          for (l=2;l<N+1;l++)
            {
              tac=A[0][l-1];
              for (m=1;m<N;m++)
                {
                  if (A[m][l-1]<tac) {tac=A[m][l-1];};
                };
              Z[l]=tac;
            };
          printf("Temps= %f, Zstop= %f \r",temps,Z[N]);
          fprintf(pf,"%f ", temps);
          for (m=2;m<N+1;m++)
            {
              fprintf(pf,"%f ",log(Z[m])/log(10));
            };
          fprintf(pf,"\n");
          if (Z[N]<thresh) {boo=1;}
        };        
    }
  while ((temps<T)&(boo==0));
  printf("\n");
  return 0;
  fclose(pf);
}