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

typedef struct {
  size_t nx;
  size_t ny;
  size_t gl_width;
  double flux;
  double gl_k;
  double gl_c;
  double epsilon;
  double delt;
  const char *savefile;
  const char *loadfile;
} cfg_t;

typedef struct {
  double **p;
  double *a;
  size_t nx, ny;
} array2_t;

cfg_t cfg = {
  .nx = 100,
  .ny = 100,
  .gl_width = 4,
  .flux = 0.1,
  .gl_k = 0.001,
  .gl_c = 8.0,
  .epsilon = 23.0e-4,
  .delt = 0.012432423423,
  .savefile = NULL,
  .loadfile = NULL,
};

int
setconfig(int argc, const char **argv)
{
  int opt;
  while ((opt = getopt(argc, (char **)argv, "hx:y:w:f:k:c:e:t:v:l:")) != -1) {
    switch (opt) {
    case 'x': cfg.nx = atoi(optarg); break;
    case 'y': cfg.ny = atoi(optarg); break;
    case 'w': cfg.gl_width = atoi(optarg); break;
    case 'f': cfg.flux = atof(optarg); break;
    case 'k': cfg.gl_k = atof(optarg); break;
    case 'c': cfg.gl_c = atof(optarg); break;
    case 'e': cfg.epsilon = atof(optarg); break;
    case 't': cfg.delt = atof(optarg); break;
    case 'v': cfg.savefile = optarg; break;
    case 'l': cfg.loadfile = optarg; break;
    case 'h':
    default:
      printf("Usage: %s [opts]\n", argv[0]);
      printf("\t-x%u\t\twidth of sim field\n", (unsigned)cfg.nx);
      printf("\t-y%u\t\theight of sim field\n", (unsigned)cfg.ny);
      printf("\t-w%u\t\twidth of glass bar\n", (unsigned)cfg.gl_width);
      printf("\t-f%-10g\tper-pixel flux\n", cfg.flux);
      printf("\t-k%-10g\trelative conductivity of glass\n", cfg.gl_k);
      printf("\t-c%-10g\trelative heat capacity of glass\n", cfg.gl_c);
      printf("\t-e%-10g\ttermination criteria\n", cfg.epsilon);
      printf("\t-t%-10g\ttime step\n", cfg.delt);
      printf("\t-v%s\tsave field at exit\n",
	cfg.savefile ? cfg.savefile : "(NULL)");
      printf("\t-l%s\tload field at beginning\n",
	cfg.loadfile ? cfg.loadfile : "(NULL)");
      printf("\t-h\t\tthis help\n");
      exit(EXIT_FAILURE);
    }
  }
}

array2_t *
array2_new(size_t nx, size_t ny)
{
  size_t i;
  array2_t *buf;
  buf = malloc(sizeof(array2_t));
  if (!buf) goto err1;
  buf->nx = nx;
  buf->ny = ny;
  buf->p = malloc(sizeof(double *) * nx);
  if (!buf->p) goto err2;
  buf->a = malloc(sizeof(double) * nx * ny);
  if (!buf->a) goto err3;
  for (i = 0; i < nx; i++) {
    buf->p[i] = (buf->a + i * ny);
  }
  return buf;
err3:
  free(buf->p);
err2:
  free(buf);
err1:
  return NULL;
}

void
array2_zero(array2_t *ap)
{
  size_t i;
  for (i = 0; i < (ap->nx * ap->ny); i++) {
    ap->a[i] = 0.0;
  }
}

int
array2_saveraw(array2_t *x, const char *filename)
{
  FILE *fp;
  size_t nelem = x->nx * x->ny;
  if ((fp = fopen(filename, "w")) == NULL) {
    perror(filename);
    goto err0;
  }
  if (fwrite(&(x->nx), sizeof(x->nx), 1, fp) < 1) goto err1;
  if (fwrite(&(x->ny), sizeof(x->ny), 1, fp) < 1) goto err1;
  if (fwrite(x->p, sizeof(double), nelem, fp) < nelem) goto err1;
  fclose(fp);
  return 0;
err1:
  fputs("write error\n", stderr);
  fclose(fp);
err0:
  return -1;
}

int
array2_loadraw(array2_t *x, const char *filename)
{
  FILE *fp;
  size_t fsize[2], nelems;
  if ((fp = fopen(filename, "r")) == NULL) {
    perror(filename);
    goto err0;
  }
  if (fread(&fsize, sizeof fsize, 1, fp) < 1) {
    fputs("read error or end-of-file\n", stderr);
    goto err1;
  }
  if (fsize[0] != x->nx) {
    fprintf(stderr, "nx %u != %u\n", (unsigned)fsize[0], (unsigned)x->nx);
    goto err1;
  }
  if (fsize[1] != x->ny) {
    fprintf(stderr, "ny %u != %u\n", (unsigned)fsize[1], (unsigned)x->ny);
    goto err1;
  }
  nelems = fsize[0] * fsize[1];
  if (fread(x->a, sizeof(double), nelems, fp) < nelems) {
    fputs("read error or end-of-file (body)\n", stderr);
    goto err1;
  }
  fclose(fp);
  return 0;
err1:
  fclose(fp);
err0:
  return -1;
}

void
array2_saveppm(array2_t *x, const char *filename)
{
  int i, j;
  double min = HUGE_VAL;
  double max = -HUGE_VAL;
  double factor;
  FILE *fp;
  for (i = 0; i < (x->nx * x->ny); i++) {
    if (x->a[i] < min) min = x->a[i];
    if (x->a[i] > max) max = x->a[i];
  }
  factor = 255.0 / (max - min);
  fp = fopen(filename, "w");
  fprintf(fp, "P6 %u %u 255\n", (unsigned)x->nx, (unsigned)x->ny);
  for (j = x->ny; j >= 0; j--) {
    for (i = 0; i < x->nx; i++) {
      int c;
      c = (int)(factor * (x->p[i][j] - min) * 4.0) % 256;
      fputc(255 - c, fp);
      fputc(255 - c, fp);
      fputc(c, fp);
    }
  }
  fclose(fp);
  printf("# save %s %lg %lg\n", filename, min, max);
}

size_t
array2_increment(array2_t *x, array2_t *dx, double epsilon)
{
  size_t c = 0;
  size_t i;
  for (i = 0; i < (x->nx * x->ny); i++) {
    if (dx->a[i] >= epsilon) c++;
    x->a[i] += dx->a[i];
  }
  return c;
}

void
setconst(array2_t **pinvc, array2_t **pkx, array2_t **pky)
{
  size_t i, j;
  size_t jbottom = cfg.ny/2 - cfg.gl_width;
  size_t jtop = cfg.ny/4 * 3;
  size_t ileft = cfg.nx/2 - cfg.gl_width;
  size_t iright = cfg.nx/2 + cfg.gl_width;
  array2_t *tmpk = array2_new(cfg.nx, cfg.ny);
  array2_t *invc = *pinvc = array2_new(cfg.nx, cfg.ny);
  array2_t *kx = *pkx = array2_new(cfg.nx, cfg.ny);
  array2_t *ky = *pky = array2_new(cfg.nx, cfg.ny);
  for (i = 0; i < invc->nx; i++) {
    for (j = 0; j < invc->ny; j++) {
      if ((i > ileft) && (i < iright) && (j > jbottom) && (j < jtop)) {
	invc->p[i][j] = cfg.gl_c * cfg.delt;
	tmpk->p[i][j] = cfg.gl_k;
      } else {
	invc->p[i][j] = 1.0 * cfg.delt;
	tmpk->p[i][j] = 1.0;
      }
    }
  }
  for (i = 0; i < (tmpk->nx - 1); i++) {
    for (j = 0; j < tmpk->ny; j++) {
      kx->p[i][j] = 0.5 * (tmpk->p[i][j] + tmpk->p[i+1][j]);
    }
  }
  for (i = 0; i < tmpk->nx; i++) {
    for (j = 0; j < (tmpk->ny - 1); j++) {
      ky->p[i][j] = 0.5 * (tmpk->p[i][j] + tmpk->p[i][j+1]);
    }
  }
}

void
boundary(array2_t *t)
{
  size_t i, j;
  for (i = 0; i < t->nx; i++) {
    t->p[i][0] = t->p[i][t->ny - 1] = 0.0;
  }
  for (j = 0; j < t->ny; j++) {
    t->p[0][j] = t->p[t->nx - 1][j] = 0.0;
  }
}

void
gettend(array2_t *t, array2_t *dt, array2_t *invc, array2_t *kx, array2_t *ky)
{
  size_t i, j;
  /* east-west flux */
  for (i = 0; i < (t->nx - 1); i++) {
    for (j = 0; j < t->ny; j++) {
      double dif = t->p[i+1][j] - t->p[i][j];
      double flux = dif * kx->p[i][j];
      dt->p[i][j] += flux * invc->p[i][j];
      dt->p[i+1][j] -= flux * invc->p[i+1][j];
    }
  }
  /* north-south flux */
  for (i = 0; i < t->nx; i++) {
    for (j = 0; j < (t->ny - 1); j++) {
      double dif = t->p[i][j+1] - t->p[i][j];
      double flux = dif * ky->p[i][j];
      dt->p[i][j] += flux * invc->p[i][j];
      dt->p[i][j+1] -= flux * invc->p[i][j];
    }
  }
  /* forcing */
  {
    size_t jbottom = t->ny/2 - cfg.gl_width - 1;
    size_t jtop = t->ny/2 + cfg.gl_width + 1;
    size_t ileft = t->nx/2 - cfg.gl_width - 1;
    size_t iright = t->nx/2 + cfg.gl_width + 1;
    for (i = ileft; i <= iright; i++) {
      dt->p[i][jbottom] += cfg.flux;
    }
    for (j = jbottom + 1; j <= jtop; j++) {
      dt->p[ileft][j] += cfg.flux;
      dt->p[iright][j] += cfg.flux;
    }
  }
}

int
run(void)
{
  array2_t *t = array2_new(cfg.nx, cfg.ny);
  array2_t *dt = array2_new(cfg.nx, cfg.ny);
  array2_t *invc;
  array2_t *kx;
  array2_t *ky;
  unsigned i;
  if (cfg.loadfile) {
    array2_loadraw(t, cfg.loadfile);
  } else {
    array2_zero(t);
  }
  setconst(&invc, &kx, &ky);
  array2_saveppm(invc, "invc.ppm");
  for (i = 0; 1; i++) {
    size_t n;
    char fnam[100];
    array2_zero(dt);
    gettend(t, dt, invc, kx, ky);
    n = array2_increment(t, dt, cfg.epsilon);
    boundary(t);
    if ((i % 100) == 0) {
      /*
      sprintf(fnam, "t%04u.ppm", (unsigned)(i / 100));
      array2_saveppm(t, fnam);
      */
      printf("# i=%u n=%u\n", (unsigned)i, (unsigned)n);
    }
    if (n == 0) {
      array2_saveppm(t, "t.ppm");
      if (cfg.savefile) {
	array2_saveraw(t, cfg.savefile);
      }
      return 0;
    }
  }
}

int
main(int argc, const char **argv)
{
  setconfig(argc, argv);
  return run();
}