/* gcc -O6 poincare.c -lm */
#include<stdio.h>
#include<stdlib.h>
#include<math.h>


#define MAX 1000

float x[MAX][2],x0[MAX][2];
float v[MAX][2],v0[MAX][2],mu[MAX];
float vmax,vcom_x, vcom_y, dt,t,test;
char sbuff[7];
int N,i,j,n, pl_handle;
long seed=1003;

main(int argc, char *argv[]){
    srand((unsigned int)getpid());
    vcom_x=0.0;
    vcom_y=0.0;
    dt=0.2;
    if(argc != 2){
	printf("./ideal N\n");
        exit(0);
    }
    vmax=2.0;
    N=atoi(argv[1]);     /* number of particles */
    /* randomly initialize the positions and velocities */
    /* and store the initial setup                      */
    for(n=0;n<N;n++){
	x[n][0]=(2.0)*(2.0*(float)rand()/(float)RAND_MAX-1.0);
	x0[n][0]=x[n][0];
	x[n][1]=(2.0)*(2.0*(float)rand()/(float)RAND_MAX-1.0);
	x0[n][1]=x[n][1];
	v[n][0]=vmax*(2.0*(float)rand()/(float)RAND_MAX-1.0);
	v[n][1]=vmax*(2.0*(float)rand()/(float)RAND_MAX-1.0);
	vcom_x=vcom_x+v[n][0];
	vcom_y=vcom_y+v[n][1];
    }
    /* subtract COM motion */
   vcom_x=vcom_x/(float)N;
   vcom_y=vcom_y/(float)N;

    for(n=0;n<N;n++){
	v[n][0]=v[n][0]-vcom_x;
	v[n][1]=v[n][1]-vcom_y;
	v0[n][0]=v[n][0];
	v0[n][1]=v[n][1];
    }
  

  for(t=0;t<2000;t++){
   
      for(n=0;n<N;n++){
	  mu[n]=0.0;  
	  x[n][0]=x[n][0]+v[n][0]*dt;
	  x[n][1]=x[n][1]+v[n][1]*dt;
	  if(x[n][0] <= -2.0 || x[n][0]>=2.0) {
	      v[n][0]=-v[n][0];  }
	  if(x[n][1] <= -2.0 || x[n][1]>=2.0) {
	      v[n][1]=-v[n][1];  }
	  /* how far are we from initial state? */
	  mu[n]=mu[n]+(x[n][0]-x0[n][0])*(x[n][0]-x0[n][0])+ \
	      (x[n][1]-x0[n][1])*(x[n][1]-x0[n][1]);
	  mu[n]=mu[n]+(v[n][0]-v0[n][0])*(v[n][0]-v0[n][0])+ \
	      (v[n][1]-v0[n][1])*(v[n][1]-v0[n][1]);
      }
      /* test for closeness, if all pass, announce it */
      test=0.0;
      for(n=0;n<N;n++)
	  test=test+mu[n];
      test=test/(float)N;
	  
      printf("%f\t%f\n", t,test);
    
  }

     
      exit(0);   
}