#include <stdio.h>
#include <string.h>
#include <math.h>

/* var */
double temp[101];	/* 結果の温度格納変数	*/

/* all */
double f;
double dtime;		/* 変化時間  Δt	*/
int ni;			/* 分割数		*/
int nim1;		/* 分割数ー１		*/

/* geom */
double x[101];		/* 計算点間隔格納変数	*/
double dxep[101];	/* E-P間の距離		*/
double dxpw[101];	/* P-W間の距離		*/
double sew[101];	/* 計算領域サイズ	*/

/* coef */
double ap[101];		/* 係数 AP	*/
double ae[101];		/* 係数 AE	*/
double aw[101];		/* 係数 AW	*/
double su[101];		/* 係数 	*/
double sp[101];		/* 係数 	*/

/* TDMA(Tri-Diagonal Matrix Algorithm) */
double tdma_temp[101];	/* TDMA受け渡し変数 温度配列	*/
int tdma_ni;		/* TDMA受け渡し変数 分割数	*/

main(argc,argv)
int argc;
char **argv;
{
  FILE *fw;
  double dx;		/* 計算対象領域のサイズ  Δx	*/
  double xmax;		/* 計算対象領域の全サイズ	*/
  double tmax;		/* 計算最大時間			*/
  double r;		/* フーリエ数			*/
  double time;		/* 計算中の時間			*/
  int i;		/* 変数				*/
  int flag;		/* 画面表示ON/OFFフラグ		*/
  int niter;		/* 計算反復回数			*/

  /* 画面表示 ON/OFF + fourier数の設定 */
  flag = 1 ;
  if ( argc >= 1 ) {
    sscanf(argv[1],"%lf",&r);
    flag = 0;
  } else {
    printf("Local fourier number [r] = ") ;
    scanf("%lf",&r) ;
  }

  /* 出力ファイルのオープン */
  if ( (fw=fopen("heat.csv","w")) == NULL ) {
    printf("file not found.[heat.dat]");
    exit(1);
  }

  /* 初期設定 */
  ni	= 21 ;			/* 21分割		*/
  nim1	= ni-1 ;
  xmax	= 1.0 ;			/* X最大寸法 1.0	*/
  tmax	= 0.5 ;			/* t最大値 0.5		*/
  dx	= xmax/(double)nim1 ;

  dtime	= r*dx*dx ;		/* Δt 			*/
/* f 	= 0.0 ;			/* f=0  :陽解法 */
  f 	= 0.5 ;			/* f=0.5:クランク・ニコルソンの陰解法 */
/*  f 	= 1.0 ;			/* f:1.0:完全陰解法 */

  /* 変数のイニシャライズ */
  init();

  /* 計算条件の表示 */
  fprintf(fw,"number of grid ,=, %d\n",ni);
  fprintf(fw,"Xmax	,=,  %f\n",xmax);
  fprintf(fw,"Tmax	,=,  %f\n",tmax);
  fprintf(fw,"Dtime	,=,  %f\n",dtime);
  for(i=1;i<=ni;i++) fprintf(fw,"%f,",x[i]);
  fprintf(fw,"\n");
  
  /* 反復回数と計算時間の初期化 */
  niter	= 0 ;
  time	= 0.0 ;
  
  /* 結果タイトルの表示＆書き込み */
  if ( flag == 1 ) {
    for(i=1;i<=ni;i++) printf("point%02d,",i);
    printf("\n");
  }

  /* ループ計算開始 */
  while ( time < tmax ) {
    niter	= niter + 1 ;
    time	= time + dtime ;

    /* calcoe 呼び出し */
    calcoe();
        tdma_ni = ni ;
        for(i=1;i<=ni;i++) tdma_temp[i] = temp[i] ;

    /* TDMA 処理 */
    tdma();
        ni = tdma_ni ;
        for(i=1;i<=ni;i++) temp[i] = tdma_temp[i] ;
    if ( flag == 1 ) {
      printf("niter	= % 4d    time	= %10.4f\n",niter,time);
      printf("%lf - %lf --- %lf - %lf \n",temp[1],temp[2],temp[ni-1],temp[ni]);
    }
    if ( flag == 1 ) prin();

    fprintf(fw,"%f,",time);
    for(i=1;i<=ni;i++) {
      fprintf(fw,"%f,",temp[i]);
    }
    fprintf(fw,"\n");
  }
  printf("calculation end.\n");
  exit(0);
}

/*
 *  変数の初期化
 */
int init()
{
  double dx ;
  int i;

  dx	= 1.0 / (double) (nim1) ;	/* 分割サイズ計算 */
  x[1]	= 0.0 ;				/* 点位置計算 */
  for(i=2;i<=ni;i++) {			/* 　｜　　　 */
    x[i]=x[i-1]+dx;			/*   ↓　　　 */
  }					/* 点位置計算 */

  dxpw[1]=0.0 ;				/* 境界と端のサイズ */
  dxep[ni]=0.0 ;			/* 境界と端のサイズ */
  for ( i=1;i<=nim1;i++) {		/* 隣接点間の距離 */
    dxep[i] = x[i+1]-x[i];		/*  　｜　　　　  */
    dxpw[i+1] = dxep[i] ;		/*    ↓　　　　  */
  }					/* 隣接点間の距離 */

  sew[1] = 0.0 ;			/* */
  sew[ni] = 0.0 ;			/* */
  for ( i=2;i<=nim1;i++ ) {		/*  */
    sew[i]=0.5 * ( dxep[i] + dxpw[i] );	/*     ↓　　　　 */
  }					/*  */

  for(i=1;i<=ni;i++ ) {			/* 初期温度 */
    temp[i]=1 ;				/*    ↓　　*/
  }					/* 初期温度 */
  return(0);
}

/*
 *  TDMA計算用係数の算出
 */
int calcoe()
{
  double f1,aeor,awor,apo ;
  int i ;

  /* 境界条件の初期化 */
  temp[1] = 0.0 ;
  temp[ni] = 0.0 ;
  ap[1] = 1.0 ;
  ae[1] = 1.0 ;
  su[1] = temp[1] ;
  ap[ni] = 1.0 ;
  aw[ni] = 0.0 ;
  su[ni] = temp[ni];

  f1 = 1.0-f ;  			/* (1-f) */

  for( i=2;i<=nim1;i++ ) {
    aeor  = 1.0 / dxep[i] ;		/* */
    awor  = 1.0 / dxpw[i] ;		/* */

    /* P13 (2.5.7) */
    apo   = sew[i] / dtime ; 		/* */
    ap[i] = f*aeor+f*awor+apo;		/* */
    ae[i] = aeor*f ;			/* */
    aw[i] = awor*f ;			/* */
    su[i] =   aeor * f1 * temp[i+1]
            + awor * f1 * temp[i-1]
            + ( apo - f1 * aeor - f1 * awor ) * temp[i];
    sp[i] = 0 ;
  }
  return(0);
}

/*
 *  TDMA 処理
 */
int tdma()
{
  double p[21],q[21],term;
  int i;

  nim1 = tdma_ni-1 ;
  p[1] = ae[1] / ap[1] ;
  q[1] = su[1] / ap[1] ;

  for(i=2;i<=tdma_ni;i++) {
    term = 1.0 / ( ap[i]-aw[i]*p[i-1] ) ;
    p[i] = ae[i] * term ; 
    q[i] = ( su[i] + aw[i] * q[i-1] ) * term ;
  }
  tdma_temp[tdma_ni] = q[tdma_ni];
  for(i=nim1;i>=2;i--) {
    tdma_temp[i] = p[i] * tdma_temp[i+1] + q[i] ;
  }
  return(0);
}

/*
 *  画面表示
 */
prin()
{
  int i;
  for(i=1;i<=ni;i++) {
    printf("%8.3f,",temp[i]);
  }
  printf("\n");
}