#include <stdio.h>
#include <math.h>
/* time.h is a standard include file for timing with clock() */
#include <time.h>
/* discretization error dominates, so let's make epsilon bigger... */
/* #define EPSILON 0.00001 ...Jacobi iterations can be less precise */
#define EPSILON 0.0001
#define N 100
#define time_steps 100
/* the following should be static values, valid for all times */
#define maxval 2.0
#define minval (-(maxval))

int main (int argc, char *argv[])
 {
  int i,j;
  int step;
  double time;
  double eps, enew;
  double time_max = 6.28;
  double alpha = 0.06;
  double dx = 1.0/N;
  double dy = 1.0/time_steps;
  double dt = time_max/time_steps;
  double dxinv = 1.0/dx;
  double dyinv = 1.0/dy;
  double dtinv = 1.0/dt;
  double divinv = 1.0/(dtinv + 2 * alpha * (dxinv * dxinv + dyinv * dyinv));
/* declare t for interior points only */
  double t[N-2][N-2];
  double told[N][N];
/* double minval, maxval; */
/* for use with POSIX clock() function */
  float stime, cputime;
  char fname[40];
  FILE *out;

  /* start the timer */
  stime = (float)clock();
  // initialize interior values
  for (i=1; i<(N-1); i++)
   for (j=1; j<(N-1); j++)
    told[i][j] = 0.0;  
   // set initial boundary conditions
   for (i=0; i<N; i++)
    {
     told[i][0]   = 0.0; // left
     told[i][N-1] = 0.0; // right
    }
   for (j=0; j<N; j++)
    told[N-1][j] = 0.0; // bottom
  // for all time steps
  for (step = 1; step <= time_steps; step++)
   {
    time = step * (time_max/time_steps);
    // reset top boundary condition each timestep
    for (j=0; j<N; j++)
     told[0][j]   =  2.0 * sin(time); // top 
    do
     {
      eps = 0.0;
      for (i=1; i<(N-1); i++)
       for (j=1; j<(N-1); j++)
        t[i-1][j-1]=((told[i][j+1]+told[i][j-1])*alpha*dyinv*dyinv+(told[i+1][j]+told[i-1][j])*alpha*dxinv*dxinv+(told[i][j]*dtinv))*divinv;
      for (i=1; i<(N-1); i++)
       {
        for (j=1; j<(N-1); j++)
         {
          enew = fabs(t[i-1][j-1] - told[i][j]);
          if (enew > eps) { eps = enew; }
         }
       }
      /* bug fix: loop should only be over the interior points */
      for (i=1; i<N-1; i++)
       for (j=1; j<N-1; j++)
        told[i][j] = t[i-1][j-1];
     }
    while(eps > EPSILON);

    // Dump raster data to a file
/* know global max and min already for all times...
   normalizing at each timestep can be deceptive...
    minval = 0.0;
    maxval = 0.0;
    for (i=0; i<N; i++)
     {
      for (j=0; j<N; j++)
       {
        if (t[i][j] < minval) { minval = t[i][j]; }
        if (t[i][j] > maxval) { maxval = t[i][j]; }
       }
     }
...don't need the loop above, comment it out! */
    /* use a well-defined image standard, Netpbm, and the Unix directory convention */
    sprintf(fname,"Output/heat%03d.pgm",step); 
    out = fopen(fname,"w+b");
    /* need a header to indicate 8-bit grayscale */
    fprintf(out, "P2 %d %d 255\n",N,N);
    for (i=0; i<N; i++)
     for (j=0; j<N; j++)
       /* bug fix: if boundaries are to be included, use told, not t */
       fprintf(out,"%d ",(int)(((told[i][j]-minval)*255.0)/(maxval - minval)));
    fclose(out);
    printf("Time step: %d\n",step);
   } // for all time steps  
  /* compute cpu time */
  cputime = ((float)clock()-stime)/CLOCKS_PER_SEC;
  printf("%d time steps in %.2f seconds\n",step-1,cputime);
  return 0;
 }