/* Molecular Dynamics */
/* Reflective walls */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>
#define N 101
#define Tmax 20000
#define ListSize 150000
#define L 500
#define a 0.01
#define v0 0.0
#define vb 1.0
#define SEED 56723
#define BINS 1000
int num_hits=0;

main()
{
int m,n,i,j;
int t, which_pair;
int *wp;
float *time;
float Etime;
float List[ListSize];
float X[N],V[N],M[N];
int bins[BINS];
float *Xptr,*Vptr,*Mptr,*Lptr;
int *binptr;
void initialize(float *Xptr,float *Vptr,float *Mptr);
void sort(float *Xptr,float *Vptr);
void get_tagged(float *Xptr,float *Vptr,float *Mptr);
void update(float *Xptr,float *Vptr,float *Mptr,int *wp);
void getmintime(float *Xptr,float *Vptr,float *time,int *wp);
void move(float *Xptr,float *Vptr,float *time);
void check_pair(float *Lptr,float *Vptr,int *wp);
void elaped_time(float time,float Etime);
void data_process(float *Xptr,float *Vptr,float *Lptr,int *binptr,int num_hits);

/* Body of the program */

for(n=0;n<N;n++){
X[n]=0.0;
V[n]=0.0;
M[n]=0.0;
}
for(n=0;n<BINS;n++)
bins[n]=0;
List[0]=0.0;
Xptr=&X[0];
Vptr=&V[0];
Mptr=&M[0];
which_pair=0;
wp=&which_pair;
*time=100000.0;
srand(SEED);
initialize( Xptr,Vptr, Mptr); /* load arrays */    
sort (Xptr,Vptr);            /* sort the particles, ascending positions */ 
get_tagged(Xptr, Vptr, Mptr); /* re-initialize tagged end particles */
Lptr=&List[0];
binptr=&bins[0];

/* We are all set up, start the time loop */
for(t=1;t<Tmax;t++){
getmintime(Xptr,Vptr,time,wp); 
move(Xptr,Vptr,time);                 
update(Xptr,Vptr,Mptr,wp);         
check_pair(Lptr,Vptr,wp);
}       /* end of time loop */

/* data output routine */
data_process(Xptr,Vptr,Lptr,binptr,num_hits);

}      /* end of the program */

void initialize(float *Xptr,float *Vptr,float *Mptr)
{
int m;
float frac;
for(m=0;m<N/2;m++){
frac=(float)rand();
frac=frac/((float)RAND_MAX);
*(Xptr+m)=frac*(float)L-(float)(L/2);
frac=(float)rand();
frac=frac/((float)RAND_MAX);
*(Xptr+m+N/2)=frac*(float)L-(float)(L/2);
*(Vptr+m)=-vb;
*(Vptr+m+N/2)=vb;
*(Mptr+m)=a;
*(Mptr+m+N/2)=a;}
return;
}

void sort(float *Xptr,float *Vptr)
{
float smallest,tmp;
int n,m;
for(n=0;n<N-1;n++){
  for(m=n+1;m<N;m++){
    if(*(Xptr+m)<*(Xptr+n)){
tmp=*(Xptr+n);
*(Xptr+n)=*(Xptr+m);
*(Xptr+m)=tmp;
tmp=*(Vptr+n);
*(Vptr+n)=*(Vptr+m);
*(Vptr+m)=tmp;}
  }}
return;
}

void get_tagged(float *Xptr,float *Vptr,float *Mptr)
{

*(Vptr+N/2)=v0;
*(Mptr+N/2)=1.0;
*Xptr=-L/2.0;
*(Xptr+N-1)=L/2.0;
*Vptr=0.0;
*(Vptr+N-1)=0.0;
*Mptr=1.0E9;  /* this makes the boundaries reflective */
*(Mptr+N-1)=1.0E9;
return;
}

void update(float *Xptr,float *Vptr,float *Mptr,int *wp)
{
float tmpx,tmpy,m1,m2,alpha,beta,gamma,den;
tmpx=*(Vptr+*wp);
tmpy=*(Vptr+*wp+1);
m1=*(Mptr+*wp);
m2=*(Mptr+*wp+1);
den=m1+m2;
alpha=(m1-m2)/den;
beta=2.0*m1/den;
gamma=2.0*m2/den;
*(Vptr+*wp)=alpha*tmpx+gamma*tmpy;
*(Vptr+*wp+1)=beta*tmpx-alpha*tmpy;
return;
}

void getmintime(float *Xptr,float *Vptr,float *time,int *wp)
{
int n,m;
float A,B;
if((*Vptr-*(Vptr+1))>0.0)
B=(*(Xptr+1)-*Xptr)/(*Vptr-*(Vptr+1));
else
B=1.0E5;
for(m=0;m<N-1;m++){
if((*(Vptr+m)-*(Vptr+m+1))>0.0)
A=(*(Xptr+m+1)-*(Xptr+m))/(*(Vptr+m)-*(Vptr+m+1));
else
A=1.0E5;
if(A<=B){
B=A;
*wp=m;}
}
*time=B;
return;
}

void move(float *Xptr,float *Vptr,float *time)
  {
int m;
for(m=1;m<N-1;m++){
*(Xptr+m)=*(Xptr+m)+*time*(*(Vptr+m));
}
return;
  }

void check_pair(float *Lptr,float *Vptr,int *wp)
  {
    if(*wp==(N/2-1) || *wp==N/2){
*(Lptr+num_hits)=*(Vptr+N/2);
/*printf("%f\n",*(Lptr+num_hits));*/
num_hits=num_hits+1;}
return;
  }

void elaped_time(float time,float Etime)
  {
Etime=Etime+time;
return;
  }

void data_process(float *Xptr,float *Vptr,float *Lptr,int *binptr,int num_hits)
{
float smallest;
float  largest;
float div;
int j,k;
int b_value;
float v_value;
float tmp;
float bound1, bound2;
float normalization;
float normal_val;
float vsum;
float v2sum;
float sigma;
smallest = *Lptr;
largest =*Lptr;
 
 for (j = 0; j < num_hits; j++)
    {
      if (*(Lptr+j) > largest)
        largest = *(Lptr+j);
      if (*(Lptr+j) < smallest)
        smallest = *(Lptr+j);
    }
                    
 /* create BINS bins */
  div = (largest - smallest) / (float)BINS;
 /* compute the normalization constant */
  normalization =div*(float)num_hits;
for( j=0;j<BINS;j++){
bound1=smallest+div*((float)j);
bound2=bound1+div;
  for (k = 0; k < num_hits; k++){
    if(*(Lptr+k)>bound1 && *(Lptr+k)<=bound2)
    *(binptr+j)=*(binptr+j)+1;
  }}
    
/* lines specific to GNU Plotutils graphing program */
printf("# The mass ratio, bath speed and tagged initial speed are \
%f\t%f\t%f\n",a,vb,v0);
printf("# There were %d  collisions involving tagged particle in time \
%ld\n",num_hits,Tmax); 
printf("# velocity limits are %f\t%f , there are %d velocty bins\n", \
smallest, largest, BINS); 
printf("# there are %d particles, system size is %d\n", N,L); 
printf("# Division size is %f The normalization constant is %f\n", \
div,normalization); 
vsum=0.0;
v2sum=0.0;
for(j=0;j<BINS;j++){
v_value=smallest+(float)j*div+div/2.0;
normal_val=(float)(*(binptr+j))/normalization;
vsum=vsum+v_value*normal_val*div;
v2sum=v2sum+v_value*v_value*normal_val*div;
printf("%f\t%f\n", v_value,normal_val);
}
sigma=sqrt(v2sum-vsum*vsum);
printf("# Standard deviation is %f\n", sigma);
}