/* 
************************************************************************
*  laplace.c:  Solution of Laplace equation with finite differences    *
*    			                                               *
*  comment: Output data is saved in 3D grid format used by gnuplot     *
************************************************************************
*/
#include <stdio.h>
#include <math.h>

#define maxx 40                         /* number of grid points */
#define maxy 40                         /* number of grid points */

main()
{
   double p[maxx][maxy], pnew[maxx][maxy], dp, dpmax, omega;
   int i, j, iter;
   
   FILE *output;			/* save data in laplace.dat */
   FILE *file;			/* save data in laplace_dp.dat */

   output = fopen("laplace.dat","w");
   file = fopen("laplace_dp.dat","w"); 

   for(i=0; i<maxx; i++)                 /* clear the array  */
   {   
      for (j=0; j<maxy; j++) p[i][j] = 0;
   }

   for(i=0; i<maxx; i++) p[i][0] = 100.0;        /* p[i][0] = 100 V */		

   for(iter=0; iter<1000; iter++)               /* iterations */
   {
      dpmax=0.0;
      for(i=1; i<(maxx-1); i++)                  /* x-direction */
      {
         for(j=1; j<(maxy-1); j++)               /* y-direction */
         {
            pnew[i][j] = 0.25*(p[i+1][j]+p[i-1][j]+p[i][j+1]+p[i][j-1]);
	    dp=pnew[i][j]-p[i][j];
	    if(fabs(dp) > dpmax) dpmax=fabs(dp);
         }
      }

      for(i=1; i<(maxx-1); i++)                  /* x-direction */
      {
         for(j=1; j<(maxy-1); j++)               /* y-direction */
         {
            p[i][j] = pnew[i][j]; 
         }
      }

      fprintf(file, "%d\t%f\n",iter,dpmax); 
   }
   
   for (i=0; i<maxx ; i++)         /* write data gnuplot 3D format */
   {
      for (j=0; j<maxy; j++) 
      {
         fprintf(output, "%f\n",p[i][j]);
      }
      fprintf(output, "\n");	  /* empty line for gnuplot */
   }
   printf("data stored in laplace.dat\n");
   fclose(output);
   fclose(file); 
}