/* 
************************************************************************
*  heat.c:   差分法を用いた熱伝導方程式の解法                          *
************************************************************************
*/

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

#define size 11			/* 格子数 */

main()
{ 
   int i,j,max;
   double pi, c, r, lambda, x, dx, dt, u[size][2];     
   
   FILE *output;			/* heat.datにデータを保存 */
   output = fopen("heat.dat","w");
   
   pi=3.1415926;

   printf("Input: r\n");
   scanf("%lf", &r);
   printf("r=%g\n", r);
   printf("Input: max\n");
   scanf("%d", &max);
   printf("max=%d\n", max);
   lambda=1.0;
   c=2.0/pi;
   dx=pi/(size-1);
   dt=r*dx*dx/lambda;

   for(j=0; j<size; j++) {
        x=dx*j;
        if(x < pi/2.0) 
	  {
            u[j][0]=c*x;
	  }
        else
	  {
	    u[j][0]=c*(pi-x);
	  }
      }

   for(i=0; i<2; i++)                 		/* 境界条件*/
   {
      u[0][i]      = 0.;
      u[size-1][i] = 0.;
   } 

   for(j=0 ; j<size; j++)  
     {
       fprintf(output, "%f\n", u[j][0]);
     }
   fprintf(output, "\n");			/* gnuplotのための空行 */
   
   for(i=1; i<=max; i++)        	/* max時間ステップにわたってループ */
   {
      for(j=1; j<(size-1); j++)         /* 空間に関するループ */
      {
         u[j][1] = u[j][0]+r*(u[j+1][0]+u[j-1][0]-2.0*u[j][0]);
      }

      for(j=0 ; j<size; j++)  
	{
	  fprintf(output, "%f\n", u[j][1]);
	}
      fprintf(output, "\n");			/* gnuplotのための空行 */
 
      for(j=0; j<size; j++) u[j][0]=u[j][1];    /* 新しい時間を古い時間に */
   }

   printf("data stored in heat.dat\n");
   fclose(output);

}