/*----------------------------------------*/
/* マクスウェルボルツマン分布(２次元)     */
/*----------------------------------------*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define n          1000        /*　粒子数              */ 
#define max_it     10000000    /*　最大時間ステップ数　*/
#define rmax       1e0         /*  領域の一辺の長さ    */
#define vmax       1e0         /*  初期速度            */
#define mass       1e0         /*  粒子一つの質量      */
#define pi         3.141592    /*  π                  */
#define nout       100         /*  出力コマ数 */ 
#define dene       0.1         /*  エネルギー分布の刻み幅 */
#define a          0.001       /*  衝突半径  */

#define tend       10.         /*  計算時間  */
#define dt         1.e-3       /*  時間刻幅  */

void est(FILE *output_file2,double x[][2], double v[][2]);  /* エネルギー分布の書き出し */

main(void)
{
  int i,j,k,num;
  int count;
  double x[n][2],v[n][2], v_tmp[n][2], t;
  double dx[2],dv[2];
  double dist;
  double tout,dtout;
  double scal;
  FILE *output_file,*output_file2;

  output_file=fopen("output.dat","w"); /* 出力ファイルオープン*/ 
  output_file2=fopen("output2.dat","w"); /* 出力ファイルオープン*/ 

  tout =0e0;        /* 出力時間変数 */
  
  dtout = tend /(double)nout;    /* 出力時間間隔 */

  count = 0;
  /* ---------------- 初期化 及び境界条件 ----------------- */

  t=0e0;

  /*   最初乱数を発生させて粒子の初期位置を決める    */

  srand(time(NULL)) ;       /* 乱数の初期化 */

  /* 0 < rnd < 1 の乱数の取得 */ 
  i=0;
  for(i=0; i < n; i++){
    for(k=0; k < 2 ; k++){
      x[i][k]= rmax * ((float)rand() / RAND_MAX); 
      v[i][k]= vmax; 
    }
  } 
  

  /* -------------- 初期化 及び境界条件終わり ------------- */


  /* -- 初期条件の出力 -- */

  num=0;
  printf("# output time = %f %d\n",t,num);
  fprintf(output_file,"# output time = %f %d\n",t,num);
  for(i=0; i < n ; i++){
    fprintf(output_file,"%f %f %f %f\n",(float)x[i][0],(float)x[i][1],
	    (float)v[i][0],(float)v[i][1]);
  }
  fprintf(output_file,"\n \n");
  fprintf(output_file2,"# output time = %f %d\n",t,num);
  est(output_file2,x,v);   /*  エネルギー分布の書き出し */
  fflush(output_file);
  fflush(output_file2);
  
  /* --------------------*/
  
  /* -------- 時間積分 --------- */    
  
  for(num=1; num < max_it; num++){
    
    for(i=0; i < n ; i++){
      for(k=0; k < 2 ; k++){
	x[i][k]=x[i][k]+v[i][k]*dt;
      }
    }
    
    /* 境界条件 壁で反射する */
    for(i=0; i < n ; i++){
      for(k=0; k < 2; k++){
	if((x[i][k]-rmax)*(x[i][k]-0.)>=0.){
	  v[i][k]=-1.*v[i][k];
	}   
      }   
    }

    /* ２粒子の衝突条件 */
    for(i=0; i < n ; i++){
      for(j=i+1; j < n ; j++){
	dist=0.;
	for(k=0; k < 2; k++){
	  dist+=pow((x[i][k]-x[j][k]),2.);
	}
	if(sqrt(dist)<= a ){
	  scal=0;
	  for(k=0; k < 2; k++){
	    dx[k]=x[j][k]-x[i][k];
	    dv[k]=v[j][k]-v[i][k];
	    scal+=dx[k]*dv[k];
	  }
	  for(k=0; k < 2; k++){
	    v_tmp[i][k]=dx[k]*scal/dist;
	    v[i][k]=v[i][k]+v_tmp[i][k];
	    v[j][k]=v[j][k]-v_tmp[i][k];
	  }
	}
      }
    }

    t=t+dt;

    /* --------            出力ルーチン              -----------*/

    if( t >= tout ){ 
      
      count+= 1;

      printf("number of iteration = %d %d\n",num,count);
      printf("time, tend, dt = %f %f %f\n",t, tend, dt);
      
      fprintf(output_file,"# output time = %f %d\n",t,num);
      for(i=0; i < n ; i++){
	fprintf(output_file,"%f %f %f %f\n",(float)x[i][0],(float)x[i][1]
		,(float)v[i][0],(float)v[i][1]);
      }
      fprintf(output_file,"\n \n");
      fprintf(output_file2,"# output time = %f %d\n",t,num);
      est(output_file2,x,v);   /*  エネルギー分布の書き出し */
      fflush(output_file);
      fflush(output_file2);
      tout+=dtout;
    }
    
    /* --------        出力ルーチン終わり              -----------*/

    if(t > tend) break; /* 計算終了判定 */
  }

  /* -------- 時間ステップ繰り返しループ 終わり----------- */    

}

void est(FILE *output_file2,double x[][2],double v[][2])
{
  int i,j,k;
  int number[100];
  double ene[n],enemax,energy[100];

  for(i=0; i < 100; i++){
    energy[i]=dene*((float)i);
  }

  enemax=0.;
  for(i=0; i < n ; i++){
    ene[i]=0.;
    for(k=0; k < 2; k++){
      ene[i]+=0.5*mass*pow(v[i][k],2.);
      if(enemax<ene[i])enemax=ene[i];
    }
  }
  
  for(i=0; i < 100; i++){
    number[i]=0;
  }

  for(i=0; i < n ; i++){
    for(j=1; j < 100; j++){
      if((ene[i]-energy[j-1])*(ene[i]-energy[j])<=0.){
	number[j]+=1;
	break;
      }
    }
  }

  for(j=1; j < 100; j++){
    fprintf(output_file2,"%f %d\n",0.5*(energy[j-1]+energy[j]),number[j]);
  }
  fprintf(output_file2,"\n \n");

}