/* 
************************************************************************
*  heat.c:   Solution of heat equation using with finite differences   *
************************************************************************
*/

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

#define size 11			/* grid size */
#define max 10 			/* iterations */

main()
{ 
   int i,j;
   double pi, c, r, lambda, x, dx, dt, u[size][2];     
   
   FILE *output;			/* save data in file eqheat.dat */
   output = fopen("heat.dat","w");
   
   pi=3.1415926;

   printf("Input: r\n");
   scanf("%lf", &r);
   printf("r=%g\n", r);
   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++)         		/* endpoints always at 0*/
   {
      u[0][i]      = 0.;
      u[size-1][i] = 0.;
   } 
   
   for(i=1; i<=max; i++)        	/* loop over max timesteps */
   {
      for(j=1; j<(size-1); j++)         /* loop over space */
      {
         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++) u[j][0]=u[j][1];    /*shift new time to old time*/
   }

   for(j=0 ; j<size; j++)  
   {
       fprintf(output, "%d\t%f\n", j, u[j][1]);
    }
   fprintf(output, "\n");			/* empty line for gnuplot */

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

}