/*----------------------------------------*/
/* Roe-Pike method   　　　　                 */
/*----------------------------------------*/
#include<stdio.h>
#include<math.h>

#define n      100      /*　一辺あたりのメッシュの数  */ 
#define Length 1.0e0     /*  矩形の一辺の長さ          */
#define pi     3.141592  /* π */
#define tcal   0.2       /* 計算時間 */
#define rhoL   1.e0       /* 左側の初期密度 */
#define pL     1.e0       /* 左側の初期圧力 */
#define vL     0.e0       /* 左側の初期速度 */
#define rhoR   0.125e0    /* 右側の初期密度 */
#define pR     0.1e0      /* 右側の初期圧力 */
#define vR     0.e0       /* 右側の初期速度 */
#define gamma  (4./3.)   /* 比熱比 */

main(void)
{
  int i,k,num;
  double delta_L,dt;
  double rho[n+2][2],v[n+2][2],p[n+2][2];
  double rhov[n+2][2],e[n+2][2],eth[n+2][2],cs[n+2][n];
  double rhoc[n+2],vc[n+2],hc[n+2],csc[n+2],pc[n+2];
  double flux[n+2][3],alpha[n+2][3],lambda[n+2][3];
  double pc1,pc2;
  double t,lambda1;
  double tout,dtout;

  FILE *output_file;
  
  /* 時間ステップなど */

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

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

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

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


  for(i=0; i < n+2 ; i++){
    for(k=0; k < 2 ; k++){
      if(i<=n/2){
	rho[i][k]=rhoL;
	v[i][k]=vL;
	p[i][k]=pL;
	rhov[i][k]=rho[i][k]*v[i][k]; /* 運動量 */
	eth[i][k]=p[i][k]/(gamma-1.e0)/rho[i][k]; /* 熱エネルギー */
	e[i][k]=rho[i][k]*(0.5e0*pow(v[i][k],2.)+eth[i][k]);/* 全エネルギー */
	cs[i][k]=sqrt(gamma*p[i][k]/rho[i][k]);/* 音速 */
      }else{
	rho[i][k]=rhoR;
	v[i][k]=vR;
	p[i][k]=pR;
	rhov[i][k]=rho[i][k]*v[i][k]; /* 運動量 */
	eth[i][k]=p[i][k]/(gamma-1.e0)/rho[i][k];/* 熱エネルギー */
	e[i][k]=rho[i][k]*(0.5e0*pow(v[i][k],2.)+eth[i][k]);/* 全エネルギー */
	cs[i][k]=sqrt(gamma*p[i][k]/rho[i][k]);/* 音速 */
      }
    }
  }

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


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



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

    num=1;

    while(t<=tcal){
      /* 時間刻幅の設定 */
      lambda1=0.;
      for(i=1; i<n+2; i++){
	if(lambda1<(fabs(v[i][k])+cs[i][k])){
	  lambda1=fabs(v[i][k])+cs[i][k];
	}
      }
      dt = 0.1*delta_L/lambda1;

      /* メッシュ境界での物理量の近似 */
    for(i=0; i < n+1 ; i++){
      rhoc[i]=sqrt(rho[i][0]*rho[i+1][0]);
      vc[i]=(sqrt(rho[i][0])*v[i][0]+sqrt(rho[i+1][0])*v[i+1][0])
	/(sqrt(rho[i][0])+sqrt(rho[i+1][0]));
      hc[i]=(sqrt(rho[i][0])*(e[i][0]+p[i][0])/rho[i][0]
	   +sqrt(rho[i+1][0])*(e[i+1][0]+p[i+1][0])/rho[i+1][0])
	/(sqrt(rho[i][0])+sqrt(rho[i+1][0]));
      csc[i]=sqrt((gamma-1.e0)*(hc[i]-0.5e0*pow(vc[i],2.)));
      pc[i]=(hc[i]-0.5e0*vc[i]*vc[i])/gamma*(gamma-1.0)*rhoc[i];
      alpha[i][0]=1./2./csc[i]/csc[i]*((p[i+1][0]-p[i][0])-rhoc[i]*csc[i]*(v[i+1][0]-v[i][0]));
      alpha[i][1]=(rho[i+1][0]-rho[i][0])-(p[i+1][0]-p[i][0])/csc[i]/csc[i];
      alpha[i][2]=1./2./csc[i]/csc[i]*((p[i+1][0]-p[i][0])+rhoc[i]*csc[i]*(v[i+1][0]-v[i][0]));
      lambda[i][0]=vc[i]-csc[i];
      lambda[i][1]=vc[i];
      lambda[i][2]=vc[i]+csc[i];
    }

      /* メッシュ境界でのフラックスの計算 */
    for(i=0; i < n+1 ; i++){
      flux[i][0]=0.5e0*(rho[i][0]*v[i][0]+rho[i+1][0]*v[i+1][0])-
	0.5e0*(alpha[i][0]*fabs(lambda[i][0])+alpha[i][1]*fabs(lambda[i][1])
	       +alpha[i][2]*fabs(lambda[i][2]));
      flux[i][1]=0.5e0*(rho[i][0]*v[i][0]*v[i][0]+p[i][0]
			+rho[i+1][0]*v[i+1][0]*v[i+1][0]+p[i+1][0])-
	0.5e0*(alpha[i][0]*fabs(lambda[i][0])*(vc[i]-csc[i])
	       +alpha[i][1]*fabs(lambda[i][1])*vc[i]
	       +alpha[i][2]*fabs(lambda[i][2])*(vc[i]+csc[i]));
      flux[i][2]=0.5e0*((e[i][0]+p[i][0])*v[i][0]+(e[i+1][0]+p[i+1][0])*v[i+1][0])-
	0.5e0*(alpha[i][0]*fabs(lambda[i][0])*(hc[i]-vc[i]*csc[i])
	       +alpha[i][1]*fabs(lambda[i][1])*0.5*vc[i]*vc[i]
	       +alpha[i][2]*fabs(lambda[i][2])*(hc[i]+vc[i]*csc[i]));
    }
      
      /* 物理量の更新 */
    for(i=1; i < n+1 ; i++){
      rho[i][1]=rho[i][0]-dt/delta_L*(flux[i][0]-flux[i-1][0]);
      rhov[i][1]=rhov[i][0]-dt/delta_L*(flux[i][1]-flux[i-1][1]);
      e[i][1]=e[i][0]-dt/delta_L*(flux[i][2]-flux[i-1][2]);
    }

    t=t+dt;
    num+=1;

    for(i=1; i < n+1 ; i++){  /* ステップの結果をコピー */
      rho[i][0]=rho[i][1];
      rhov[i][0]=rhov[i][1];
      e[i][0]=e[i][1];
    }

      /* 新しい時間での物理量の更新 */
    for(i=1; i < n+1 ; i++){  
      v[i][0]=rhov[i][0]/rho[i][0];
      p[i][0]=(e[i][0]-0.5e0*rhov[i][0]*v[i][0])*(gamma-1.e0);
    }
    
    for(i=1; i < n+1 ; i++){  
      for(k=0; k<2; k++){
	rhov[i][k]=rho[i][k]*v[i][k];
	eth[i][k]=p[i][k]/(gamma-1.e0)/rho[i][k];
	e[i][k]=rho[i][k]*(0.5e0*pow(v[i][k],2.)+eth[i][k]);
	cs[i][k]=sqrt(gamma*p[i][k]/rho[i][k]);
      }
    }
    
  /* --------            出力ルーチン              -----------*/

    if(t>=tout){  /* t >= tout となれば出力 */
      printf("number of iteration = %d %f %f\n",num,t,dt);
  
      fprintf(output_file,"#number of iteration = %d %f\n",num,t);
      for(i=0; i < n+2 ; i++){
	fprintf(output_file,"%f %f %f %f\n",(float)delta_L*i,rho[i][0],v[i][0],p[i][0]);
      }
      fprintf(output_file,"\n\n");
      tout+=dtout;  /* toutの更新 */
    }
  /* --------        出力ルーチン終わり              -----------*/

  }

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

}