/*----------------------------------------*/
/* 自己重力系のビリアル定理　　　　       */
/*----------------------------------------*/
#include<stdio.h>
#include <stdlib.h>
#include<math.h>

#define n          50          /*　粒子数              */ 
#define max_it     10000000    /*　最大時間ステップ数　*/
#define rmax       1e0         /*  初期半径            */
#define grav       1e0         /*  重力定数            */
#define mass_tot   1e0         /*  系の全質量          */
#define pi         3.141592    /*  π                  */
#define eps        rmax*1e-2   /*  重力softening length */
#define nout       100         /*  出力コマ数 */ 
 
void force(double mass, double x[][3], double v[][3], double fx[][3], double fv[][3],double pot[]);

main(void)
{
  int i,j,k,num;
  int count;
  double x[n][3],v[n][3], x_tmp[n][3],v_tmp[n][3], t, dt;
  double fx[n][3],fv[n][3],f_finx[n][3],f_finv[n][3];
  double k1x[n][3],k2x[n][3],k3x[n][3],k4x[n][3];
  double k1v[n][3],k2v[n][3],k3v[n][3],k4v[n][3];
  double pot[n];
  double radi;
  double tff,rho_ini,tend,tout,dtout,mass;

  FILE *output_file;

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

  mass = mass_tot/(double)n;    /* 一粒子質量 */

  rho_ini = mass_tot / (4e0*pi/3e0*pow(rmax,3.));  /* 初期密度 */

  tff = sqrt(3e0*pi/(32e0*grav*rho_ini));  /* 計算時間の目安 */

  tend = 10e0*tff ; /* 計算終了時刻 */

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

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

  t=0e0;
  dt=1e-5;

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

  srand(time(NULL)) ;       /* 乱数の初期化 */
  /* 0 < rnd < 1 の乱数の取得 */ 
  i=0;
  while(i < n){
    for(k=0; k < 3 ; k++){
      x[i][k]= 2.*rmax * ((float)rand() / 2147483648.-5e-1); 
    }
    radi = sqrt(x[i][0]*x[i][0] + x[i][1]*x[i][1] + x[i][2]*x[i][2]);

    v[i][0]=0e0;
    v[i][1]=0e0;
    v[i][2]=0e0;
    if( radi <= rmax ) i=i+1; 

    } 


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


  /* -- 初期条件の出力 -- */
  
  for(i=0; i < n ; i++){
    fprintf(output_file,"%f %f %f\n",(float)x[i][0],(float)x[i][1],(float)x[i][2]);
  }
  fprintf(output_file,"\n");

  
  /* --------------------*/
  
  /* -------- 4次ルンゲクッタの時間積分 --------- */    
  
  for(num=1; num < max_it; num++){
    
    if( count == 1 ){
      tout =tout + dtout;
      count=0;
    }

    force(mass,x,v,fx,fv,pot);

    
    for(i=0; i < n ; i++){
      for(k=0; k < 3 ; k++){
	k1x[i][k]=fx[i][k]*dt; 
	x_tmp[i][k]=x[i][k]+k1x[i][k]*5e-1;
	
	k1v[i][k]=fv[i][k]*dt; 
	v_tmp[i][k]=v[i][k]+k1v[i][k]*5e-1;
      }
    }
    
    force(mass,x_tmp,v_tmp,fx,fv,pot);
    
    for(i=0; i < n ; i++){
      for(k=0; k < 3 ; k++){
	
	k2x[i][k]=fx[i][k]*dt; 
	x_tmp[i][k]=x[i][k]+k2x[i][k]*5e-1;
	
	k2v[i][k]=fv[i][k]*dt; 
	v_tmp[i][k]=v[i][k]+k2v[i][k]*5e-1;
      }
    }
    
    force(mass,x_tmp,v_tmp,fx,fv,pot);
    
    for(i=0; i < n ; i++){
      for(k=0; k < 3 ; k++){
	
	k3x[i][k]=fx[i][k]*dt; 
	x_tmp[i][k]=x[i][k]+k3x[i][k];
	
	k3v[i][k]=fv[i][k]*dt; 
	v_tmp[i][k]=v[i][k]+k3v[i][k];
      }
    }
    
    force(mass,x_tmp,v_tmp,fx,fv,pot);
    
    for(i=0; i < n ; i++){
      for(k=0; k < 3 ; k++){
	
	k4x[i][k]=fx[i][k]*dt; 
	k4v[i][k]=fv[i][k]*dt; 
      }
    }
    
    for(i=0; i < n ; i++){
      for(k=0; k < 3 ; k++){
	
	x[i][k] = x[i][k] + 1e0 / 6e0 * ( k1x[i][k] + 2e0*k2x[i][k] + 2e0*k3x[i][k] + k4x[i][k] );
	f_finx[i][k] = 1e0 / 6e0 * ( k1x[i][k] + 2e0*k2x[i][k] + 2e0*k3x[i][k] + k4x[i][k] ) / dt;
	
	v[i][k] = v[i][k] + 1e0 / 6e0 * ( k1v[i][k] + 2e0*k2v[i][k] + 2e0*k3v[i][k] + k4v[i][k] );
	f_finv[i][k] = 1e0 / 6e0 * ( k1v[i][k] + 2e0*k2v[i][k] + 2e0*k3v[i][k] + k4v[i][k] ) / dt;
      }
    }

    t=t+dt;

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

    /* printf("time, tend, dt = %f %f %f\n",t, tout, dt); */

    if( t > tout ){ 
      
      count = 1;

      printf("number of iteration = %d\n",num);
      printf("time, tend, dt = %f %f %f\n",t, tend, dt);

      
      for(i=0; i < n ; i++){
	/* printf("i,fv,fx = %d %f %f \n",i,fv[i][k],fx[i][k]); */
	fprintf(output_file,"%f %f %f \n",(float)x[i][0],(float)x[i][1],(float)x[i][2]);
      }
      fprintf(output_file,"\n");
    }

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

    /* --------        time step control ---------- */
    dt = 5e-1;
      
    for(i=0; i < n ; i++){
      
      for(k=0; k < 3 ; k++){
	
	if( abs(f_finx[i][k]) > 1e-8 && fabs(x[i][k]) > 1e-8 ){
	    if( dt  > 2e-1 * fabs(eps / f_finx[i][k])){
	      dt =  2e-1* fabs(eps / f_finx[i][k]) ;
	    }
	}
	
	if( abs(f_finv[i][k]) > 1e-8 && fabs(v[i][k]) > 1e-8 ){
	  if( dt  > 2e-1 * fabs(v[i][k] / f_finv[i][k])){
	    dt =  2e-1 * fabs(v[i][k] / f_finv[i][k]) ;
	  }
	}
      }
    }
    
    /* --------        time step control 終わり ---------- */


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

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

}

void force(double mass,double x[][3], double v[][3], double fx[][3], double fv[][3],double pot[])
{
  int i,j,k;
  double r[n],r2[n];

    for(i=0; i < n ; i++){

      pot[i]=0.;

      for( k=0; k < 3 ; k++){
	fv[i][k]=0.;
      }

      for(j=0; j < n ; j++){

	r2[i]=0e0;

	for( k=0; k < 3 ; k++){
	  r2[i] = r2[i] + (x[i][k]-x[j][k])*(x[i][k]-x[j][k]);
	}

	r[i] = sqrt(r2[i]);

	if(i != j){

	  pot[i] = pot[i] - 1./r[i];

	  for( k=0; k < 3 ; k++){
	    
	    fv[i][k]= fv[i][k]+(x[j][k]-x[i][k])/(r[i]*r2[i]+eps*eps*eps);

	  }
	}
      }



      pot[i] = pot[i]*grav*mass;

      for( k=0; k < 3 ; k++){
	
	fv[i][k]=fv[i][k]*grav*mass;
	fx[i][k]=v[i][k];

	/*	printf("i,fv,fx = %d %f %f \n",i,fv[i][k],fx[i][k]); */
		  
      }
      
    }

}