/*----------------------------------------*/
/* モンテカルロ      　　                 */
/*----------------------------------------*/
#include<stdio.h>
#include<stdlib.h>
#include<math.h>

#define n      20        /*　一辺あたりのメッシュの数  */ 
#define T0_amp 10.       /*  初期温度分布の中心温度と熱浴の温度の比  */ 
#define T_edge 1.0e1     /*  境界の熱浴の温度           */ 
#define Length 1.0e1     /*  矩形の一辺の長さ          */
#define kappa  1e0       /* κ */
#define pi     3.141592  /* π */
#define tcal   5.        /* 計算時間 */

#define nT     10      /* 温度1をエネルギー粒子何個に変換するか */
#define sigma  .3        /* 正規分布の標準偏差 */

void normrand(double *,double *); /* 正規乱数を作る関数 */

main(void)
{
  int i,j,k,num;
  int flag;
  double delta_L,dt;
  double tmp[n+2][2], t;
  double pos[nT*(n*2+(int)(T0_amp))*(int)(T_edge)][2];
  double tout,dtout;
  double dx,dy;
  int number[n+2][2],nall,nn;

  FILE *output_file;
  
  /* 乱数の発生 */
  srand(time(NULL)) ;       /* 乱数の初期化 */

  /* 時間ステップなど */

  delta_L = (float)Length/((float)n+1);
  dt = sigma*sigma/2./kappa;
  
  output_file=fopen("outmont1D.dat","w"); /* 出力ファイルオープン*/ 

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

  t=0;
  tout=0.1;     /* 計算を出力する時間 */
  dtout=0.1;    /* 計算を出力する時間刻み */

  /*     境界はメッシュ番号 0 とメッシュ番号 n+1 にとる     */

  
  for(i=0; i < n+2 ; i++){
    for(k=0; k < 2 ; k++){
      tmp[i][k]=T_edge;               /* 初期化及び境界の値をT_edgeに */
      number[i][k]=nT*(int)T_edge;    /* 各メッシュに存在する粒子数 */
    }
  }

  tmp[n/2][0]=T0_amp*T_edge;              /* 初期条件設定・真ん中一点が最初高温 */
  number[n/2][0]=nT*(int)(T0_amp*T_edge); /* 中心メッシュに存在する粒子数 */

  /* 各粒子の初期位置 */
  nall=0; 
  for(i=0; i < n+2 ; i++){
    for(nn=0; nn < number[i][0]; nn++){
      pos[nall][0]=(float)delta_L * i;
      nall+=1;
    }
  }
  
  /* -------------- 初期化 及び境界条件終わり ------------- */


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

  for(i=0; i < n+2 ; i++){
    number[i][1]=0;
  }
  
  for(nn=0; nn < nall ; nn++){
    for(i=0; i < n+2 ; i++){
      if((pos[nn][0]-delta_L*((float)i-0.5))*(pos[nn][0]-delta_L*((float)i+0.5))<=0.){
	number[i][1]+=1;
      }
    }
  } 
  
  for(i=0; i < n+2 ; i++){
    tmp[i][1]=(float)number[i][1]/(float)nT;
  }  
  
  for(i=0; i < n+2 ; i++){
    fprintf(output_file,"%f %f %d\n",(float)delta_L*i,
	    (float)tmp[i][1],number[i][1]);
  }
  fprintf(output_file,"\n");
  fprintf(output_file,"\n");

  /* --------------------*/



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

  num=1;

  while(t<=tcal){

    for(nn=0; nn < nall ; nn++){
      normrand(&dx,&dy);
      pos[nn][0]+=dx;
    }

    t=t+dt;
    num+=1;

    /* 境界条件 */ 
    for(i=0; i < n+2 ; i++){
      number[i][1]=0;
    }
    for(nn=0; nn < nall ; nn++){
      i=0;
      if((pos[nn][0]-delta_L*((float)i-0.5))*(pos[nn][0]-delta_L*((float)i+0.5))<=0.){
	number[i][1]+=1;
	if(number[i][1]>nT*(int)T_edge)pos[nn][0]=-100.;
      }
      i=n+1;
      if((pos[nn][0]-delta_L*((float)i-0.5))*(pos[nn][0]-delta_L*((float)i+0.5))<=0.){
	number[i][1]+=1;
	if(number[i][1]>nT*(int)T_edge)pos[nn][0]=-100.;
      }
    } 
    nn=0;
    i=0;
    while(number[i][1]<nT*(int)T_edge){
      nn+=1;
      number[i][1]+=1;
      pos[nall+nn][0]=(float)delta_L * i;
    }
    i=n+1;
    while(number[i][1]<nT*(int)T_edge){
      nn+=1;
      number[i][1]+=1;
      pos[nall+nn][0]=(float)delta_L * i;
    }
    nall+=nn;

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

    if(t>=tout){   /* t >= tout となれば出力 */
      printf("number of iteration = %d %f\n",num,t);

      /* どのメッシュに何個粒子が存在するか */
      for(i=0; i < n+2 ; i++){
	number[i][1]=0;
      }
      for(nn=0; nn < nall ; nn++){
	for(i=0; i < n+2 ; i++){
	  if((pos[nn][0]-delta_L*((float)i-0.5))*(pos[nn][0]-delta_L*((float)i+0.5))<=0.){
	      number[i][1]+=1;
	  }
	}
      } 

      /* 粒子数を温度に変換 */
      for(i=0; i < n+2 ; i++){
	tmp[i][1]=(float)number[i][1]/(float)nT;
      }  

      fprintf(output_file,"#number of iteration = %d %f\n",num,t);
      for(i=0; i < n+2 ; i++){
	fprintf(output_file,"%f %f %d\n",(float)delta_L*i,
		(float)tmp[i][1],number[i][1]);
      }	
      fprintf(output_file,"\n");
      fprintf(output_file,"\n");
      tout+=dtout;    /* toutの更新 */
    }

  /* --------        出力ルーチン終わり              -----------*/

  }

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

}


/* --------正規乱数を与える関数--------- */
void normrand(double *g1,double *g2){
  double a,b;

  a=((float)rand() / RAND_MAX);
  while (a ==0){
    a=((float)rand() / RAND_MAX);
  }
  b=((float)rand() / RAND_MAX);
  
  *g1=sigma*sqrt(-2.0*log(a))*cos(2.0*pi*b);
  *g2=sigma*sqrt(-2.0*log(a))*sin(2.0*pi*b);
}