/*
################################################################

Modified version of http://psoup.math.wisc.edu/papers/snslowof.c

Changes by Martin Krzywinski 
http://mkweb.bcgsc.ca/snowflakes

- added command-line parameters and default values
- removed X11 dependency

################################################################

USAGE

snowflake -PARAM VALUE [-q] [-h]

where PARAM is one of

 # Gravner-Griffeath model parameters
 a FLOAT alpha
 b FLOAT beta
 k FLOAT kappa
 g FLOAT gamma
 m FLOAT mu
 r FLOAT rho
 t FLOAT theta
 s FLOAT sigma (noise, set to 0 for deterministic output)

 x INT   random seed, only useful if s is not zero

 n INT   size of canvas (max 2024)
 y INT   maximum size of flake (default n)
 i INT   max number of iterations

and

 -q quiet, do not report iteration progress
 -h help

If a parameter is not specified, a default value will be used. The default values and
suggested ranges for parameters are

  a = 0.208625 // 0         - 0.3
  b = 1.105423 // 1.05      - 2.2
  k = 0.020542 // 0         - 0.05
  g = 0.000077 // 0.0000515 - 0.0001
  m = 0.017995 // 0         - 0.01
  r = 0.603418 // 0.3       - 0.9
  t = 0.008316 // 0.02      - 0.57
  s = 0        // -0.5      - 0
  n = 500      // 50        - 1000
  i = 25000    // >0

*/

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

#define T 1
#define F 0

// must set ulimit -s 30000
// to avoid segfaults for values of nrmax and ncmax > 1024
// even with segfaults happen for > 2030

#define nrmax 2024
#define ncmax 2024

double adif[nrmax][ncmax]; /* diffusion field              */ 
int    apic[nrmax][ncmax]; /* indicator of snowflake sites */
double afr[nrmax][ncmax];  /* boundary mass                */ 
double alm[nrmax][ncmax];  /* crystal mass                 */
int    ash[nrmax][ncmax];  /* rings pallette               */

double beta, kappa, mu, theta, alpha, gam, sigma, rho, rhorinit; 
int    seed;

int nr,nc;
int max_flake_size, flake_margin;
int change, centeri,centerj;
int parupdate; 
int parash; 
int rold, rnew, rinit; 
int frchange; 
int twelvesided; 
int noac,pq;
int stop;

double myrand() {
  double drand48();
  return drand48();
}

int norminf(i,j) int i,j; {
  if (i<0) i=-i; 
  if (j<0) j=-j; 
  if (i>j) {
    return i; 
  } else {
    return j;
  }
}

int seminorm(i,j) int i,j; {
  int k; 
  k = i + j; 
  if (k >= 0) {
    return k; 
  } else { 
    return -k;
  }
}

int check_state() {
  int i,j;
  double tol = 0.0001;
  int reached_tol = 1;
  for (i=0; i<=9; i++) {
    for(j=0; j<=9; j++) {
      if(adif[i][j] > tol) reached_tol = 0;
      //printf("%.5lf|", adif[i][j]);
    }
    //printf("\n"); 
  }
  //printf("\n");
  return reached_tol;  
}

void initialize() {
  int t1,t2;
  int i,j,k;
  double x;
  double drand48();
  double x1, y1; 

  pq        = 0;
  stop      = F;
  parupdate = 0; 

  if(! seed) {
    srand48(time(NULL));
    t1   = time(&t2);    
    srand48(time(NULL));
    seed = t1%1000;
  }
  for (i=1; i < seed; i++){ 
    x = drand48();
  }

  centeri = nr / 2;
  centerj = nc / 2; 

  rold=0; 
  rnew=0; 

  for (i=0;i<nr;i++){
    for (j=0; j<nc;j++){
      x = myrand();
      if (twelvesided == 0){
	if ((norminf(i-centeri,j-centerj)<=rinit)   &&
	    (seminorm(i-centeri, j-centerj)<=rinit) &&
	    (x<=rhorinit)) {
	  adif[i][j] = 0.0; 
	  apic[i][j] = 1; 
	  afr[i][j]  = 1; 
	  ash[i][j]  = 0; 
	  alm[i][j]  = 0.0; 
	  k = norminf(i-centeri, j-centerj); 
	  if (k>rnew) {
	    rnew=k; 
	  }
	} else {
	  adif[i][j] = rho; 
	  apic[i][j] = 0; 
	  afr[i][j]  = 0.0;
	  ash[i][j]  = 0; 
	  alm[i][j]  = 0.0; 
	}
      } else {
	x1=(double)(i-centeri)/rinit; 
	y1=(double)(j-centerj)/rinit; 
	if (  ( ((i-centeri) == -(j-centerj))&&(i-centeri<=0)&&(i-centeri>=-rinit) ) ||
	      ( (i-centeri>=0) && (i-centeri<=rinit) && (j-centerj==0) ) ||
	      ( (j-centerj<=0) && (j-centerj>=-rinit) && (i-centeri==0) ) ) {
	  adif[i][j] = 0.0; 
	  apic[i][j] = 1; 
	  afr[i][j]  = 0.0; 
	  ash[i][j]  = 0; 
	  alm[i][j]  = 1.0; 
	  k = norminf(i-centeri, j-centerj); 
	  if (k>rnew) {
	    rnew=k;
	  }
	} else { 
	  adif[i][j] = rho; 
	  apic[i][j] = 0; 
	  afr[i][j]  = 0.0;
	  ash[i][j]  = 0;
	  alm[i][j]  = 0.0; 
	}
      }
    }
  } 
  rold   = rnew;
  parash = 1;   
}
 
void dynamicsdif()
{   

  double b[nrmax][ncmax];
  double x,y;
  int i,j,k;
  int id,iu,jl,jr;
  int jend; 
  int part;
  int count; 
  double b1, b2; 
  double masscorrection; 
  int nrhalf; 

  for (i=0;i<nr;i++)
    for (j=0;(j<nc);j++){
      b[i][j]=0.0; 
    }
    
  for (i=0;i<nr;i++){
    for (j=0;j<nc;j++){
      if (apic[i][j]==0){
	id=(i+1)%nr;iu=(i+nr-1)%nr; 
	jr=(j+1)%nr; jl=(j+nr-1)%nr;  
	count=0; 
	if (apic[id][j]==0) count++; 
	if (apic[iu][j]==0) count++;
	if (apic[i][jl]==0) count++;  
	if (apic[i][jr]==0) count++;
	if (apic[iu][jr]==0) count++;
	if (apic[id][jl]==0) count++;  
    
	if (count==0) b[i][j]=adif[i][j]; 
	else{
	  b[i][j]=(1.0-(double)count/7.0)*adif[i][j]+
	    (adif[id][j]*(1.0-apic[id][j])+adif[iu][j]*(1.0-apic[iu][j])
	     +adif[i][jl]*(1.0-apic[i][jl])+adif[i][jr]*(1.0-apic[i][jr])
	     +adif[iu][jr]*(1.0-apic[iu][jr])+adif[id][jl]*(1.0-apic[id][jl]))/7.0; 
	}                    
      }  
    }
  }
     
  for (i=0;i<nr;i++){
    for (j=0;(j<nc);j++){
      if (apic[i][j]==0) adif[i][j]=b[i][j]; 
    }
  }  
}

void dynamicspop()
{   
  double b[nrmax][ncmax];
  double x,y;
  int i,j,k;
  int id,iu,jl,jr;
  int part;
  int count;
  double offset;  

  for (i=0;i<nr;i++){
    for (j=0;j<nc;j++){
      x=myrand();                           
      if (x<0.5) {
	adif[i][j]=adif[i][j]*(1+sigma);
      } else {
	adif[i][j]=adif[i][j]*(1-sigma);
      }
    }       
         
  }           
}
  
int dynamicsunfre()
{   
  double x,y,afrij;
  int i,j,k;
  int id,iu,jl,jr;
  int part;
  int count; 
  double offset; 
  double difmass; 
  
  int ilo,iup,jlo,jup;
  
  ilo=centeri-rnew-1; 
  iup=centeri+rnew+1; 
  jlo=centerj-rnew-1;
  jup=centerj+rnew+1;

  frchange = F; 

  if( (ilo < 0) || (jlo < 0) )  { return T; }      //ilo = 0;
  if( 2*rnew > max_flake_size-2 ) { return T; }

  for (i=ilo;i<=iup;i++){   
    for (j=jlo;j<=jup;j++){
      if (apic[i][j]==0) { 
	afrij      = afr[i][j]; 
	y          = afrij*mu;        	  
	// printf("%d %d\n",i,j);
	afr[i][j]  = afr[i][j]-y; 
	adif[i][j] = adif[i][j]+y;
	afrij      = alm[i][j]; 
	if (afrij>0.0){
	  y          = afrij*gam;        	  
	  alm[i][j]  = alm[i][j]-y; 
	  adif[i][j] = adif[i][j]+y;
	}	    	      
      }	   
    } 
  }   
  return F;
}

int dynamicsfre()
{   

  int bpic[nrmax][ncmax];
  
  double x,y,afrij;
  int i,j,k;
  int id,iu,jl,jr;
  int part;
  int count; 
  double offset; 
  double nfrsum; 
  double frmass, difmass;
  
  int ilo,iup,jlo,jup;
 
  ilo=centeri-rnew-1; 
  iup=centeri+rnew+1; 
  jlo=centerj-rnew-1; 
  jup=centerj+rnew+1;  
  frchange=F; 
  
  if( (ilo<0) || (jlo<0) ) { return T; }
  //if(ilo<0) ilo = 0;
  //if(jlo<0) jlo = 0;
  //if(iup>nrmax-1) iup = nrmax-1;
  //if(jup>nrmax-1) jup = ncmax-1;

  for (i=ilo;i<=iup;i++)    
    for (j=jlo;j<=jup;j++){
      bpic[i][j]=apic[i][j];
    }
  for (i=ilo;i<=iup;i++){   
    for (j=jlo; j<=jup;j++){
      if (apic[i][j]==0){
	
	id=(i+1) % nr;iu=(i+nr-1) % nr; 
	jr=(j+1) % nc; jl=(j+nc-1)%nc; 
	count=0; 
	if (apic[id][j]==1) count++; 
	if (apic[iu][j]==1) count++;
	if (apic[i][jl]==1) count++;  
	if (apic[i][jr]==1) count++;
	if (apic[iu][jr]==1) count++;
	if (apic[id][jl]==1) count++;
	
	if (count>=1){
	  
	  difmass=adif[i][j]+adif[id][j]*(1-apic[id][j])+adif[iu][j]*(1-apic[iu][j])
	    +adif[i][jl]*(1-apic[i][jl])+adif[i][jr]*(1-apic[i][jr])
	    +adif[iu][jr]*(1-apic[iu][jr])+adif[id][jl]*(1-apic[id][jl]);	       
	  
	  if (count<=2) { 
	    if (afr[i][j]>=beta){ bpic[i][j]=1;} 
	  }
 
	  if (count>=3){ 
	    if ( (afr[i][j]>=1.0) || ((difmass<=theta)&&(afr[i][j]>=alpha)) ) {
	      bpic[i][j]=1;                    
	    } 
	  }
	  if (count>=4) bpic[i][j]=1;   
	}	    	     
      }	            
    }       
  }
  for (i=ilo;i<=iup;i++){    
    for (j=jlo;j<=jup;j++){
          
      if (apic[i][j]!=bpic[i][j]){
	apic[i][j]=bpic[i][j]; 

	alm[i][j]+=afr[i][j]; afr[i][j]=0.0; 	     
	k=norminf(i-centeri, j-centerj);
	if (k>rnew) rnew=k; if (rnew>2*nr/3) stop=T; 
	ash[i][j]=parash; 
	frchange=T; 
      }
    }
  }
  parupdate=1-parupdate; 
  if (rnew-rold==1){ parash=parash+1; rold=rnew;}
  return F;
}

int dynamicsfre1()
{   
  
  double x,y;
  int i,j,k;
  int id,iu,jl,jr;
  int part;
  int count; 
  double offset; 
  double frmass; 
  double blockmass;
  
  int ilo,iup,jlo,jup;
  double epsilon; 
   
  ilo=centeri-rnew-1; 
  iup=centeri+rnew+1; 
  jlo=centerj-rnew-1; 
  jup=centerj+rnew+1;
  frchange=F; 

  if( (ilo<0) || (jlo<0) ) { return T; } //ilo = 0;
  //if(ilo<0) ilo = 0;
  //if(jlo<0) jlo = 0;
  //if(iup>nrmax-1) iup = nrmax-1;
  //if(jup>nrmax-1) jup = ncmax-1;

  for (i=ilo;i<=iup;i++){   
    for (j=jlo;j<=jup;j++){
      if (apic[i][j]==0){  
	id=(i+1) % nr;iu=(i+nr-1) % nr; 
	jr=(j+1) % nc; jl=(j+nc-1)%nc; 
	count=0; 
	if (apic[id][j]==1) count++; 
	if (apic[iu][j]==1) count++;
	if (apic[i][jl]==1) count++;  
	if (apic[i][jr]==1) count++;
	if (apic[iu][jr]==1) count++;
	if (apic[id][jl]==1) count++;
	if (count>=1){
	  offset=(1.0-kappa)*adif[i][j];   
	  afr[i][j]=afr[i][j]+offset;  
	  offset=adif[i][j]-offset; adif[i][j]=0;  alm[i][j]+=offset; 
	}	     
      }        
    }
  }   
  return F;
}

// Grow the flake in four steps
// - diffusion
// - freezing 
// - melting
int dynamics()
{
  dynamicsdif();    
  // these wi*ll return 1 if we need to stop computing
  int f1 = dynamicsfre1();
  int f2 = dynamicsfre();         
  int f3 = dynamicsunfre();   
  if (sigma>0.0) dynamicspop(); 
  //printf("%d %d %d\n",f1,f2,f3);
  return f1 || f2 || f3;
}

// STDOUT report of flake state
// y Y z1 z2 z3 ...
//
// where zi is the value of the pixel at y=Y and x=i
// z < 0 diffuse mass
// z > 0 flake mass
void report_flake_state()
{
  int i,j,kf;
  float  k;
  double y; 
  for (i=0; i <nr; i++){
    printf("y %d ",i);
    for (j=0; j<nc; j++){
      if (apic[i][j] == 0){ 
	// greys
	k = -255.0*(adif[i][j]/rho);                  
	if (k < -255.0) k = -255.0; 
      } else {
	// blues
 	y = alm[i][j] + adif[i][j]; 
	k = (255.0*y-alpha)/(beta-alpha); 
	if (k > 255.0) k = 255.0; 
      }       
      printf("%.0f ",round(k));
    }
    printf("\n");
  }
} 

main(argc,argv) int argc; char *argv[]; {

  int i,j,k,nop;
  int reached_end;
  int posx,posy;
  int rootx,rooty;

  char ch;

  int aflag = 0;
  int bflag = 0;
  int index;
  int c;
  int silent = F;
  int maxiter;
  int t0;

  rinit    = 0;
  rhorinit = 1;

  alpha   = 0.208625; // 0 - 0.3
  beta    = 1.105423; // 1.05 - 3
  kappa   = 0.020542; // 0 - 0.05
  gam     = 0.000077; // 0.0000515 - 0.0001
  mu      = 0.017995; // 0 - 0.01
  rho     = 0.603418; // 0.3 - 0.9
  sigma   = 0;        // -0.5 - 0
  theta   = 0.008316; // 0.02 - 0.57
  nr      = 500;
  maxiter = 25000;
  seed    = 0;

  max_flake_size = nr;
  twelvesided    = 0;

  if (rinit<0) {
    rinit        =- rinit;
    twelvesided  = 1;
  }
 
  while ((c = getopt(argc,argv,"a:b:k:g:m:r:s:t:n:y:z:i:x:hq")) != -1)
    switch (c) {
    case 'a':
      alpha = atof(optarg);
      break;
    case 'b':
      beta = atof(optarg);
      break;
    case 'k':
      kappa = atof(optarg);
      break;
    case 'g':
      gam = atof(optarg);
      break;
    case 'm':
      mu = atof(optarg);
      break;
    case 'r':
      rho = atof(optarg);
      break;
    case 's':
      sigma = atof(optarg);
      break;
    case 't':
      theta = atof(optarg);
      break;
    case 'n':
      nr = atoi(optarg);
      break;
    case 'y':
      max_flake_size = atoi(optarg);
      break;
    case 'z':
      flake_margin = atoi(optarg);
      break;
    case 'i':
      maxiter = atoi(optarg);
      break;
    case 'x':
      seed = atoi(optarg);
      break;
    case 'q':
      silent = T;
      break;
    case 'h':
      printf("\n");
      printf("Generate Gravner-Griffeath Snowflakes\n");
      printf("\n");
      printf("*************************************************************************\n");
      printf("\n");
      printf("Adapted by Martin Krzywinski (http://mkweb.bcgsc.ca/snowflakes) from Modeling Snow Crystal Growth II\n");
      printf("by Gravner and Griffeath http://psoup.math.wisc.edu/papers/snslowof.c\n");
      printf("\n");
      printf("*************************************************************************\n");
      printf("\n");
      printf("snowflake -a X  -b X -k X -g X -m X -r X -t X -s X -x Y -n Y -y Y -z Y -i Y\n");
      printf("\n");
      printf("where X is a float and Y is an integer. The parameters and default values:\n");
      printf("\n");
      printf("a 0.208625 alpha\n");
      printf("b 1.105423 beta\n");
      printf("k 0.020542 kappa\n");
      printf("g 0.000077 gamma\n");
      printf("m 0.017995 mu\n");
      printf("r 0.603418 rho\n");
      printf("t 0.008316 theta\n");
      printf("s 0.0      sigma (noise, set to 0 for deterministic output)\n");
      printf("n 500      size of canvas (max 2024)\n");
      printf("y 500      maximum size of flake (default n)\n");
      printf("z 0        flake margin (n-y, set either y or z)\n");
      printf("i 25000    max number of iterations\n");
      printf("x          random seed, only useful if s is not zero\n");
      printf("\n");
      printf("Default values will be used for parameters that are not defined.\n");
      printf("\n");
      printf("*************************************************************************\n");
      printf("\n");
      printf("For examples see http://mkweb.bcgsc.ca/snowflakes\n");
      printf("\n");
      exit(0);
    default:
      abort();
    }

  if(max_flake_size > nr) max_flake_size = nr;
  if(flake_margin) max_flake_size        = nr - 2*flake_margin;
  
  printf("param a %.10f\n",alpha);
  printf("param b %.10f\n",beta);
  printf("param g %.10f\n",gam);
  printf("param k %.10f\n",kappa);
  printf("param m %.10f\n",mu);
  printf("param r %.10f\n",rho);
  printf("param s %.10f\n",sigma);
  printf("param t %.10f\n",theta);
  printf("param n %d\n",nr);
  printf("param y %d\n",max_flake_size);
  printf("param z %d\n",flake_margin);
  printf("param i %d\n",maxiter);

  nc = nr;
  pq = 0; 
  initialize();

  printf("param x %d\n",seed);

  t0 = time(NULL);
  j  = 0;
  while(pq!=-1) {
    noac=0;
    pq++;
    reached_end = dynamics();
    reached_end = reached_end || check_state();
    if(pq > maxiter) reached_end = 1;
    if(reached_end) break;
    if (! silent && pq%100 == 0){
      printf("iter %d\n",pq);
    }
    j++;
  }
  if(j>maxiter) j = maxiter;
  printf("run j %d\n",j);
  printf("run time %d\n",time(NULL)-t0);
  printf("run time_per_iter %.6f\n",(double)(time(NULL)-t0)/j);

  report_flake_state();
}