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

#define n      50        /*　一辺あたりのメッシュの数  */ 
#define max_sum 500      /*　級数のオーダーの最大値　　*/
#define phi0   1.0e0      /*  かける電圧                */ 
#define Length 1.0e1      /*  矩形の一辺の長さ          */
#define pi     3.14159263 /* π　*/

int main(void)
{
  int i,j,k;
  double delta_L;
  double phi[n+2][n+2];
  double x,y;
  /* ↑全部で (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;               /* 初期化及び境界の値をグラウンドに */
    }
  }

 

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

  delta_L = (float)Length/(n+1);


  for(i=0; i < n+2 ; i++){
    for(j=0; j < n+2 ; j++){
      x=delta_L*i;y=delta_L*j;
      for(k=1; k < max_sum+1 ; k++){
	phi[i][j] += 2.*(1.-pow(-1,k))/(k*pi*(exp(k*pi)-exp(-k*pi)))*sin(k*pi*x/Length)*(exp(k*pi*y/Length)-exp(-k*pi*y/Length));
				       }			       
    }
  }
  
  /* --------            出力ルーチン              -----------*/


  output_file=fopen("output2.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);
}