/*----------------------------------------*/
/* 1Dポテンシャル反射・透過問題　　　　   */
/*----------------------------------------*/
#include<stdio.h>
#include <stdlib.h>
#include<math.h>

#define nin        1000          /*　grid数                 */ 
#define nbuf       nin/10       /*  片側バッファーのgrid数 */
#define n          nin+2*nbuf    /*  全grid数               */
#define max_it     10000000      /*　最大時間ステップ数　   */
#define Length     300           /*  計算領域               */
#define dL         Length/(double)(nin+1) /* mesh 間隔        */
#define pi         3.141592      /*  π                     */
#define nout       300           /*  出力コマ数             */ 
 
void bibun(double psi[][2], double fpsi[][2]);
double pot(double x);
double dump(double x);
void shoki(double x, double *psi_r, double *psi_i);

main(void)
{
  int i,j,k,num;
  int count;
  double psi[n+2][2], psi_tmp[n+2][2], fpsi[n+2][2];
  double k1[n+2][2],k2[n+2][2],k3[n+2][2],k4[n+2][2];
  double t, dt, tend, tout, dtout;
  double x, psi_r, psi_i;

  FILE *output_file;

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


  tend = 2e2 ; /* 計算終了時刻 */

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

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

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

  for(i=1; i < n+1 ; i++){
    x=(double)(i-nbuf)*dL;
    shoki(x,&psi_r,&psi_i);
    psi[i][0] = psi_r; psi[i][1] = psi_i;
  }

  psi[0][0] = 0.;
  psi[0][1] = 0.;
  psi[n+1][0] = 0.;
  psi[n+1][1] = 0.;

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


  /* -- 初期条件の出力 -- */
  
  for(i=0; i < n+2 ; i++){
    x=(double)(i-nbuf)*dL;
    fprintf(output_file,"%f %f %f %f\n",(float)x,(float)psi[i][0],(float)psi[i][1],(float)pot(x));
  }
  fprintf(output_file,"\n");
  
  /* --------------------*/
  
  /* -------- 4次ルンゲクッタの時間積分 --------- */    
  
  for(num=1; num < max_it; num++){
    
    if( count == 1 ){
      tout =tout + dtout;
      count=0;
    }

    bibun(psi,fpsi);    
    for(i=1; i < n+1 ; i++){
      for(k=0; k < 2 ; k++){
	k1[i][k]=fpsi[i][k]*dt; 
	psi_tmp[i][k]=psi[i][k]+k1[i][k]*5e-1;	
      }
    }

    bibun(psi_tmp,fpsi);    
    for(i=1; i < n+1 ; i++){
      for(k=0; k < 2 ; k++){
	k2[i][k]=fpsi[i][k]*dt; 
	psi_tmp[i][k]=psi[i][k]+k2[i][k]*5e-1;	
      }
    }

    bibun(psi_tmp,fpsi);    
    for(i=1; i < n+1 ; i++){
      for(k=0; k < 2 ; k++){
	k3[i][k]=fpsi[i][k]*dt;
	psi_tmp[i][k]=psi[i][k]+k3[i][k];
      }
    }

    bibun(psi_tmp,fpsi);    
    for(i=1; i < n+1 ; i++){
      for(k=0; k < 2 ; k++){
	k4[i][k]=fpsi[i][k]*dt; 
      }
    }
    
    
    for(i=1; i < n+1 ; i++){
      for(k=0; k < 2 ; k++){
	
	psi[i][k] = psi[i][k] + 1e0 / 6e0 * ( k1[i][k] + 2e0*k2[i][k] + 2e0*k3[i][k] + k4[i][k] );
	
      }
    }


    t=t+dt;

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

    if( t > tout ){ 
      
      count = 1;

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

      
      for(i=0; i < n+2 ; i++){
	x=(double)(i-nbuf)*dL;
	fprintf(output_file,"%f %f %f %f\n",(float)x,(float)psi[i][0],(float)psi[i][1],(float)pot(x));
      }
      fprintf(output_file,"\n");
    }

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

    /* --------        time step control ---------- */
    dt = 0.5* dL * dL;
    /* --------        time step control 終わり ---------- */


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

  }

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

}
/* ------- main 終り ------------------------------ */




/* --------微係数を与える関数--------- */
void bibun(double psi[][2], double fpsi[][2]){
  double x;
  int i;
  for(i=1; i < n+1 ; i++){
    x = (double)(i-nbuf)*dL;
    fpsi[i][0] = -(psi[i+1][1] + psi[i-1][1] - 2.*psi[i][1])/(dL*dL) + pot(x)*psi[i][1]-dump(x)*psi[i][0];
    fpsi[i][1] =  (psi[i+1][0] + psi[i-1][0] - 2.*psi[i][0])/(dL*dL) - pot(x)*psi[i][0]-dump(x)*psi[i][1];
  }
  return;
}

/* ------- Potential ------------ */
double pot(double x){
  double L1,L2, ans;
  L1=0.5;L2=0.6;
  ans=0.;
  if(x > L1*Length && x < L2*Length ){
    ans = 0.9;
  }
  return ans;
}

/*------- dumping function -----------*/
double dump(double x){
  double ans;
  ans=0.;
  if(x < 0. || x > Length ){
    ans = 0.3;
  }
  return ans;
}

/*------- initial condition ------------*/
void shoki(double x, double *psi_r, double *psi_i){
  double L1,L2;
  L1=0.;L2=L1+0.1;
  *psi_r=0.;*psi_i=0.;

  
  if(x > L1*Length && x < L2*Length ){
    *psi_r = cos(x)/sqrt((L2-L1)*Length); 
    *psi_i = sin(x)/sqrt((L2-L1)*Length);
  }
  
  return;
}