/*----------------------------------------*/
/* 2次元ラプラス方程式の解                */
/*----------------------------------------*/
#include<stdio.h>
#include<math.h>

#define n      50        /*　一辺あたりのメッシュの数  */ 
#define max_it 10000      /*　最大計算回数　　　　　　　*/
#define eps    1.0e-6     /*　誤差の閾値                */
#define phi0   1.0e0      /*  かける電圧                */ 
#define Length 1.0e1      /*  矩形の一辺の長さ          */

int main(void)
{
  int i,j,num;
  int flag;
  double error;
  double delta_L;
  double phi[n+2][n+2],phi_old[n+2][n+2]; 
  /* ↑全部で (n+2)^2の格子点だが、端は境界条件でfixするので、*/
  /*   実際に解くのは中のn^2個であることに注意                */

  FILE *output_file;
  

  /* ---------------- 初期化 及び境界条件 ----------------- */
  /*     境界はメッシュ番号 0 とメッシュ番号 n+1 にとる     */


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

    for(j=0; j < n+2 ; j++){
      phi[i][j]=0.0;               /* 初期化及び境界の値をグラウンドに */
    }
  }

  for(i=0; i < n+2 ; i++){
    phi[i][n+1]=phi0;                /* y=L に電圧を掛ける */
    }


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

  /* -------- ガウスザイデルの繰り返しループ --------------- */    

  for(num=0; num < max_it; num++){ 
    
    flag = 0;

    for(i=1; i < n+1 ; i++){
      for(j=1; j < n+1 ; j++){
	phi_old[i][j]=phi[i][j]; /* 前のステップの値を保存 */
      }
    }
    
    for(i=1; i < n+1 ; i++){
      for(j=1; j < n+1  ; j++){
	
	error=0.0;

	phi[i][j]=0.25*(phi[i+1][j]+phi[i-1][j]+phi[i][j+1]+phi[i][j-1]);
	/* ↑ ラプラス方程式の差分方程式 */
	
	error = fabs(phi[i][j]-phi_old[i][j])/fabs(phi[i][j]+phi_old[i][j]);
	if( error > eps ) flag = 1;
	/* ↑ 各メッシュでの誤差の評価。誤差が大きければフラグを立てる */

      }
    }
    
    if( flag == 0 ) break; /* フラグが立っていなければループ離脱 */
    
  }

  /* -------- ガウスザイデルの繰り返しループ 終わり----------- */    


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

  printf("number of iteration = %d\n",num);
  
  delta_L = (float)Length/(n+1);

  output_file=fopen("output.dat","w");
  for(i=0; i < n+2 ; i++){
    for(j=0; j < n+2 ; j++){
    fprintf(output_file,"%f %f %f\n",(float)delta_L*i,(float)delta_L*j,phi[i][j]);
    }
    fprintf(output_file,"\n");
  }
  /* --------        出力ルーチン終わり              -----------*/

  return (0);
}