/* diffusion_seq.c: sequential version of 2d diffusion.
 * Author: Stephen F. Siegel <siegel@cis.udel.edu>, August, 2007.
 */
#include <stdlib.h>
#include <stdio.h>
#include <assert.h>
#define SQUARE(x) ((x)*(x))

/* Constants: the following should be defined at compilation:
 *
 *       D = diffusion constant
 *      dt = time step
 *      dx = distance between two lattice points
 *       r = radius of initial circle in center
 *  nsteps = number of time steps
 *      nx = number of points in x direction
 *      ny = number of points in y direction
 */

/* Global variables */
double u[nx][ny];  /* values of the discretized function */

/* init: initializes global variables. */
void init() {
  int i, j;
  
  for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
      u[i][j] =
        (SQUARE(dx*i-0.5*dx*(nx-1)) + SQUARE(dx*j-0.5*dx*(ny-1))
         < r*r ? 100 : 0);
}

/* write_frame: writes current values of u to file */
void write_frame(int time) {
  FILE *file;
  char filename[50];
  char command[50];
  int i,j;
  
  sprintf(filename, "./seqout/out_%d", time);
  file = fopen(filename, "w");
  assert(file);
  for (j = 0; j < ny; j++) {
    for (i = 0; i < nx; i++) fprintf(file, "%8.2f", u[i][j]);
    fprintf(file, "\n");
  }
  fclose(file);
  sprintf(command, "cat %s", filename);
  system(command);
  sleep(3);
}

/* updates u for next time step. */
void update() {
  int i, j;
  double u_new[nx][ny];
  double k = D*dt/(dx*dx);

  for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
      u_new[i][j] =  u[i][j]
        + k*(u[(i+1)%nx][j] + u[(i+nx-1)%nx][j]
             +u[i][(j+1)%ny] + u[i][(j+ny-1)%ny] - 4*u[i][j]);
  for (i = 0; i < nx; i++)
    for (j = 0; j < ny; j++)
      u[i][j] = u_new[i][j];
}

/* main: executes simulation, creates one output file for each time
 * step */
int main(int argc,char *argv[]) {
  int iter;

  init();
  write_frame(0);
  for (iter = 1; iter <= nsteps; iter++) {
    update();
    write_frame(iter);
  }
  return 0;
}
