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

/* =================================================================================
    The distribution of heat over time is described by the so called heat equation:
    df/dt = alpha*(d2f/dx2 + d2f/dy2)   for a function f(x,y,t).
    This formula may be discretized in a regular grid G(:,:) by computing the new
       value G1(x,y) in a point (x,y) at each time step as:
       G1(x,y) = G(x,y) + Ax * ( G(x+1,y) + G(x-1,y) - 2.0 * G(x,y))
  		+ Ay * ( G(x,y+1) + G(x,y-1) - 2.0 * G(x,y))
    For each point in the grid the next value depends on the values of the four up
       and down, left and right adjacent points.
 
    This solution is based on a source code written by 
       Blaise Barney, D. Turner, George L. Gusciora
  ****************************************************************************
  Explanation of constants and variables
    Nrows                          =  x dimension of problem grid 
    Ncols                          =  y dimension of problem grid
    Tsteps                         =  number of time Tsteps
    Ax, Ay                         =  used in heat equation
    Grid                           =  array for grids
  ---------------------------------------------------------------------------
*/
#define Nrows 2000
#define Ncols 2000
#define Tsteps 1000
float *Grid;

#define Ax 0.1
#define Ay 0.1

// to address C arrays in Fortran-like column order
#define index2D(a,b,LD1)  a + (b*LD1)
#define index3D(a,b,c,LD1,LD2)  a + (b*LD1) + (c*LD1*LD2)
      
void compute(int rows, int cols, float *g1, float *g2)
{
/*
    Compute next step using matrices g1, g2 of dimension (nr,nc)
*/

      int i, j, nr=Nrows;

      for ( j = 1; j < cols-1; j++ ) {
        for ( i = 1; i < rows-1; i++ ) {
            g2[index2D(i,j,nr)] = g1[index2D(i,j,nr)] 
             + Ax * ( g1[index2D((i+1),j,nr)] + g1[index2D((i-1),j,nr)] - 2.0 * g1[index2D(i,j,nr)]) 
             + Ay * ( g1[index2D(i,(j+1),nr)] + g1[index2D(i,(j-1),nr)] - 2.0 * g1[index2D(i,j,nr)]);
        }
      }
}


void inidat()
{
/*
      Initialize Grid data
*/
      int i,j;

      for ( i = 0; i < Nrows*Ncols*2; i++ ) Grid[i] = 0.0;

      for ( i = 0; i < Nrows; i++ ) {
        for ( j = Ncols/2-1; j <= Ncols/2+1; j++ ) {
          Grid[index3D(i,j,0,Nrows,Ncols)] = (float)i*(float)(Nrows-i-1) + (float)j*(float)(Ncols-j-1);
        }
      }
      
      for ( i = Ncols/2-5; i < Ncols/2+5; i++ ) {
        for ( j = Ncols/2-5; j <= Ncols/2+5; j++ ) {
          Grid[index3D(i,j,0,Nrows,Ncols)] = (float)i*(float)(Nrows-i-1) + (float)j*(float)(Ncols-j-1);
        }
      }
     
     return; 
}


int MINVAL(int s, float* a)
{
   int e;
   float v;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v > a[e] ) v = a[e];
   
   }
   
   return(v);
}

int MAXVAL(int s, float* a)
{
   int e;
   float v;

   v = a[0];

   for ( e = 0; e < s; e++ )
   {
       if ( v < a[e] ) v = a[e];
   
   }
   
   return(v);
}
   
   

void RealData2pgm(int s1, int s2, float *rdata, char* name) 
{
   /*
      Simple subroutine to dump integer data in a PGM format
   */
   int i, j;
   FILE* ouni;
   int  vp, vs;
   float  rmin, rmax;
   char  fname[80];
   /*
      Write on unit 700 with  PGM format
   */
   strcpy(fname,name);
   strcat(fname,".pgm\0");
   ouni = fopen(fname,"w");
   if ( ! ouni ) {
      fprintf(stderr,"!!!! Error write access to file %s\n",fname);
   }
   /*  Magic code */
   fprintf(ouni,"P2\n");
   /*  Dimensions */
   fprintf(ouni,"%d %d\n",s1, s2);
   /*  Maximum value */
   fprintf(ouni,"255\n");
   /*  Values from 0 to 255 */
   rmin = MINVAL(s1*s2,rdata); rmax = MAXVAL(s1*s2,rdata);
   vs = 0;
   for ( i = 0; i < s1; i++ ) {
      for ( j = 0; j < s2; j++ ) {
         vp = (int) ( (rdata[i+(j*s1)] - rmin) * 255.0 / (rmax - rmin) );
         vs++;
         fprintf(ouni,"%4d ", vp);
         if ( vs >= 10 ) {
            fprintf(ouni," \n");
            vs = 0;
         }
      }
      fprintf(ouni," ");
      vs = 0;
   }
   fclose(ouni);
   return;
}

int main( int argc, char *argv[] )   /* Program Heat2D */
{

      char fname[80];
      int it, sp;

//    Allocate and initialize grid
      Grid = (float*)malloc(sizeof((float)1.0)*Nrows*Ncols*2);
      if ( Grid == NULL ) { 
              fprintf(stderr," Error allocating Grid(%d,%d,2)",Nrows,Ncols);
              return(-1);
       }
       inidat();

//          Dump initial configuration
       it = 0;
       sprintf(fname,"tmap%4.4d",it);
       RealData2pgm(Nrows,Ncols,Grid,fname);
 
       printf("Start computing ... \n");       
          sp=0;
          for ( it = 1; it <= Tsteps; it++ ) {
       
//           Calculate  grid points
           compute(Nrows,Ncols,
                    &Grid[index3D(0,0,sp,Nrows,Ncols)],
                    &Grid[index3D(0,0,(1-sp),Nrows,Ncols)]);
            sp=1-sp;

         }

//       Dump final configuration
            it = Tsteps;
            sprintf(fname,"tmap%4.4d",it);
            RealData2pgm(Nrows,Ncols,&Grid[index3D(0,0,sp,Nrows,Ncols)],fname);

         printf("  ... computed!\n");
         return(0);
 
}