/***********************************************************************
 *
 * File: fire.c (21-Jun-2001)
 *              (21-Jun-2005)
 *
 * Waldbrand-Simulation (Gittermodell fuer selbstorganisierte
 * Kritizitaet), nach B. D. Malamud & D. L. Turcotte, Computing in
 * Science and Engineering 2/3, pp. 42-15 (2000)
 *
 * nx=ny............Gitterdimension
 * f_m..............Zuendfrequenz [nach je (1/f_m - 1) Baumpflanzungen
 *                  wird ein Zuendholz weggeworfen]
 * nsteps...........Anzahl der Zuendversuche pro Lauf
 * "fire_in.dat"....Startkonfiguration fuer Lauf
 * "fire_out.dat"...Endkonfiguration des Laufs
 *
 **********************************************************************/

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

#define min(a,b) ((a)<=(b) ? (a):(b))

/* Globale Variablen */

FILE *fp;
int *ihist,*nlistx,*nlistx_old,*nlisty,*nlisty_old;
int **lattice;
int istep,ns,nsteps,nx,ny;

/**********************************************************************/

void init() {     /* Lauf initialisieren */

/**********************************************************************/
  
  unsigned short int iran[3];
  int ix,mode;

  /* Eingabeaufforderung */

  mode=0;
  while(mode<1 || mode>3) {
    printf(" Select mode=1,2,3 to\n");
    printf("   1...start new simulation\n");
    printf("   2...re-initialize simulation\n");
    printf("   3...continue simulation\n\n");
    printf("   mode=");
    scanf("%d",&mode);
  }

  if(mode==1) {     /* Neue Simulation starten */
    
    printf("  nx=ny=");
    scanf("%d",&nx);
    ny=nx;

    ns=0;
    while(ns<2) {
      printf("  1/f_m=");
      scanf("%d",&ns);
    }

    printf(" nsteps=");
    scanf("%d",&nsteps);

    ihist=(int*)malloc((nx*ny+1)*sizeof(int));
    lattice=(int**)malloc((nx+2)*sizeof(int*));
    for(ix=0;ix<=nx+1;ix=ix+1) {
      lattice[ix]=(int*)malloc((ny+2)*sizeof(int));
    }

    istep=0;
    
  } else if(mode==2) {     /* Simulation re-initialisieren */

    printf(" nsteps=");
    scanf("%d",&nsteps);
    
    fp=fopen("fire_in.dat","r");

    fread(&nx,sizeof(int),1,fp);
    ny=nx;

    ihist=(int*)malloc((nx*ny+1)*sizeof(int));
    lattice=(int**)malloc((nx+2)*sizeof(int*));
    for(ix=0;ix<=nx+1;ix=ix+1) {
      lattice[ix]=(int*)malloc((ny+2)*sizeof(int));
    }
    
    fread(&istep,sizeof(int),1,fp);
    fread(&ns,sizeof(int),1,fp);
    fread(iran,sizeof(unsigned short int),3,fp);
    
    for(ix=1;ix<=nx;ix=ix+1) {
      fread(&lattice[ix][1],sizeof(int),ny,fp);
    }
    
    fclose(fp);

    seed48(iran);     /* RNG mit seed von File initialisieren */
    
    istep=0;
    
  } else if(mode==3) {     /* Simulation fortsetzen */
    
    printf(" nsteps=");
    scanf("%d",&nsteps);
    
    fp=fopen("fire_in.dat","r");
    
    fread(&nx,sizeof(int),1,fp);
    ny=nx;

    ihist=(int*)malloc((nx*ny+1)*sizeof(int));
    lattice=(int**)malloc((nx+2)*sizeof(int*));
    for(ix=0;ix<=nx+1;ix=ix+1) {
      lattice[ix]=(int*)malloc((ny+2)*sizeof(int));
    }
    
    fread(&istep,sizeof(int),1,fp);
    fread(&ns,sizeof(int),1,fp);
    fread(iran,sizeof(unsigned short int),3,fp);

    for(ix=1;ix<=nx;ix=ix+1) {
      fread(&lattice[ix][1],sizeof(int),ny,fp);
    }

    fread(&ihist[1],sizeof(int),nx*ny,fp);
    
    fclose(fp);

    seed48(iran);     /* RNG mit seed von File initialisieren */

  }
  
  nlistx=(int*)malloc((nx*ny+1)*sizeof(int));
  nlistx_old=(int*)malloc((nx*ny+1)*sizeof(int));
  nlisty=(int*)malloc((nx*ny+1)*sizeof(int));
  nlisty_old=(int*)malloc((nx*ny+1)*sizeof(int));

}

/**********************************************************************/

void tree() {     /* Baum an zufaelligem Gitterpunkt pflanzen */

/**********************************************************************/

  int ix,iy;
  double ran;
  
  ran=drand48();
  ix=min((int)(ran*nx+1.0),nx);
  ran=drand48();
  iy=min((int)(ran*ny+1.0),ny);

  lattice[ix][iy]=1;
  
}

/**********************************************************************/

void match() {     /* Zuendholz auf zufaelligen Gitterpunkt werfen */

/**********************************************************************/

  int i,ix,iy,nburn,nn,nn_old;
  double ran;
  
  ran=drand48();
  ix=min((int)(ran*nx+1.0),nx);
  ran=drand48();
  iy=min((int)(ran*ny+1.0),ny);

  
  if(lattice[ix][iy]==1) {     /* Baum auf Gitterplatz? */
    
    /* Cluster von Baeumen bestimmen, die Feuer fangen */
    
    nn=1;                  /* Brandherd */
    nlistx[1]=ix;
    nlisty[1]=iy;
    lattice[ix][iy]=0;
    
    nburn=0;
    while(nn>=1) {         /* Nachbarn u. Nachbarn der Nachbarn... */

      nburn=nburn+nn;
      nn_old=nn;
      for(i=0;i<=nx*ny;i=i+1) {
	nlistx_old[i]=nlistx[i];
	nlisty_old[i]=nlisty[i];
      }

      nn=0;
      for(i=1;i<=nn_old;i=i+1) {
	ix=nlistx_old[i];
	iy=nlisty_old[i];
	if(lattice[ix+1][iy]==1) {
	  lattice[ix+1][iy]=0;
	  nn=nn+1;
	  nlistx[nn]=ix+1;
	  nlisty[nn]=iy;
	}
	if(lattice[ix-1][iy]==1) {
	  lattice[ix-1][iy]=0;
	  nn=nn+1;
	  nlistx[nn]=ix-1;
	  nlisty[nn]=iy;
	}
	if(lattice[ix][iy+1]==1) {
	  lattice[ix][iy+1]=0;
	  nn=nn+1;
	  nlistx[nn]=ix;
	  nlisty[nn]=iy+1;
	}
	if(lattice[ix][iy-1]==1) {
	  lattice[ix][iy-1]=0;
	  nn=nn+1;
	  nlistx[nn]=ix;
	  nlisty[nn]=iy-1;
	}
      }
      
    }
    
    ihist[nburn]=ihist[nburn]+1;     /* Histogramm */
      
  }
  
}

/**********************************************************************/

void finish() {     /* Lauf abschliessen */

/**********************************************************************/

  unsigned short int *iran;
  unsigned short int dummy[3];
  int ix;
  
  iran=seed48(dummy);     /* RNG seed auslesen */

  fp=fopen("fire_out.dat","w");
  
  fwrite(&nx,sizeof(int),1,fp);
  fwrite(&istep,sizeof(int),1,fp);
  fwrite(&ns,sizeof(int),1,fp);
  fwrite(iran,sizeof(unsigned short int),3,fp);

  for(ix=1;ix<=nx;ix=ix+1) {
    fwrite(&lattice[ix][1],sizeof(int),ny,fp);
  }
  
  fwrite(&ihist[1],sizeof(int),nx*ny,fp);
  
  fclose(fp);
  
}

/**********************************************************************/

int main() {

/**********************************************************************/

  int i,j;
  
  init();
  
  for(i=1;i<=nsteps;i=i+1) {
    istep=istep+1;
    for(j=1;j<=ns-1;j=j+1) {
      tree();
    }
    match();     /* Nach je (ns-1) Baumpflanzungen ein Zuendholz */
  }

  finish();
  
  return 0;
  
}

/**********************************************************************/