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

#include <time.h> 

#include <X11/Xlib.h>
#include <X11/Xutil.h>

#define T 1
#define F 0

#define nrmax 1002
#define ncmax 1002

#define kappamax 64  

int sp;   

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*/

/* parameters*/

double beta, kappa, mu, theta, alpha, gam, sigma; 

int nr,nc;
  
int change, centeri,centerj;

int parupdate; 
 
int parash; 
 
int rold, rnew; 

int rinit; 
double rhorinit;
 
int frchange; 
int twelvesided; 


double rho;

int stop;

char po[30]; 

Display *td;
Window tw;
GC tgc;
XEvent te;
KeySym tk;
XSizeHints th;
int ts;
unsigned long tf,tb;
char tbuf[8];
int kc;
int fin;

char timestring[]="time:";
char activestring[]="active area:";

int noac,pq;

Colormap cmap;
XColor clr[kappamax];

XColor clron[128];
XColor clroff[128];
XColor othp[20];

int flags={DoRed|DoGreen|DoBlue};

char in[]="sn";
char wn[]="digital snowflake";

int red[125], green[125],blue[125];

FILE *picf;

FILE *prof; 

char infile[30]; 

char outfile[30];

char graphicsfile[30];

char comments[100]; 

void bluecolors33()
{
  int i;
  red[0]= 71; green[0]=204; blue[0]=231;
  red[1]= 70; green[1]=200; blue[1]=230;
  red[2]= 69; green[2]=196; blue[2]=229;
  red[3]= 68; green[3]=192; blue[3]=228;
  red[4]= 67; green[4]=188; blue[4]=227;
  red[5]= 66; green[5]=184; blue[5]=226;
  red[6]= 65; green[6]=180; blue[6]=225;
  red[7]= 64; green[7]=176; blue[7]=224;
  red[8]= 63; green[8]=172; blue[8]=223;
  red[9]= 62; green[9]=168; blue[9]=222;
  red[10]= 61; green[10]=164; blue[10]=221;
  red[11]= 60; green[11]=160; blue[11]=220;
  red[12]= 59; green[12]=156; blue[12]=219;
  red[13]= 58; green[13]=152; blue[13]=218;
  red[14]= 57; green[14]=148; blue[14]=217;
  red[15]= 56; green[15]=144; blue[15]=216;
  red[16]= 55; green[16]=140; blue[16]=215;
  red[17]= 54; green[17]=136; blue[17]=214;
  red[18]= 53; green[18]=132; blue[18]=213;
  red[19]= 52; green[19]=128; blue[19]=212;
  red[20]= 51; green[20]=124; blue[20]=211;
  red[21]= 50; green[21]=120; blue[21]=210;
  red[22]= 49; green[22]=116; blue[22]=209;
  red[23]= 48; green[23]=112; blue[23]=208;
  red[24]= 47; green[24]=108; blue[24]=207;
  red[25]= 46; green[25]=104; blue[25]=206;
  red[26]= 45; green[26]=100; blue[26]=205;
  red[27]= 44; green[27]=96; blue[27]=204;
  red[28]= 43; green[28]=92; blue[28]=203;
  red[29]= 42; green[29]=88; blue[29]=202;
  red[30]= 41; green[30]=84; blue[30]=201;
  red[31]= 40; green[31]=80; blue[31]=200;
  red[32]= 10; green[32]=20; blue[32]=100;
}

void braquecolors64()
{
  int i;
  red[0]= 130; green[0]=166; blue[0]=167;
  red[1]= 140; green[1]=176; blue[1]=186;
  red[2]= 156; green[2]=193; blue[2]=200;
  red[3]= 163; green[3]=204; blue[3]=212;
  red[4]= 167; green[4]=213; blue[4]=212;
  red[5]= 169; green[5]=207; blue[5]=215;
  red[6]= 168; green[6]=207; blue[6]=211;
  red[7]= 157; green[7]=199; blue[7]=205;
  red[8]= 145; green[8]=162; blue[8]=155;
  red[9]= 122; green[9]=137; blue[9]=151;
  red[10]= 114; green[10]=128; blue[10]=136;
  red[11]= 101; green[11]=130; blue[11]=142;
  red[12]= 102; green[12]=130; blue[12]=157;
  red[13]= 96; green[13]=129; blue[13]=162;
  red[14]= 96; green[14]=130; blue[14]=165;
  red[15]= 98; green[15]=131; blue[15]=166;
  red[16]= 130; green[16]=166; blue[16]=167;
  red[17]= 140; green[17]=176; blue[17]=186;
  red[18]= 156; green[18]=193; blue[18]=200;
  red[19]= 163; green[19]=204; blue[19]=212;
  red[20]= 167; green[20]=213; blue[20]=212;
  red[21]= 169; green[21]=207; blue[21]=215;
  red[22]= 168; green[22]=207; blue[22]=211;
  red[23]= 157; green[23]=199; blue[23]=205;
  red[24]= 137; green[24]=175; blue[24]=189;
  red[25]= 130; green[25]=166; blue[25]=174;
  red[26]= 118; green[26]=152; blue[26]=164;
  red[27]= 118; green[27]=153; blue[27]=157;
  red[28]= 118; green[28]=158; blue[28]=160;
  red[29]= 123; green[29]=164; blue[29]=166;
  red[30]= 136; green[30]=177; blue[30]=177;
  red[31]= 146; green[31]=191; blue[31]=197;
  red[32]= 106; green[32]=140; blue[32]=143;
  red[33]= 116; green[33]=162; blue[33]=160;
  red[34]= 142; green[34]=183; blue[34]=185;
  red[35]= 184; green[35]=201; blue[35]=205;
  red[36]= 189; green[36]=224; blue[36]=229;
  red[37]= 214; green[37]=248; blue[37]=247;
  red[38]= 224; green[38]=232; blue[38]=240;
  red[39]= 201; green[39]=228; blue[39]=234;
  red[40]= 187; green[40]=218; blue[40]=222;
  red[41]= 170; green[41]=194; blue[41]=203;
  red[42]= 138; green[42]=175; blue[42]=176;
  red[43]= 115; green[43]=153; blue[43]=162;
  red[44]= 101; green[44]=137; blue[44]=151;
  red[45]= 83; green[45]=126; blue[45]=132;
  red[46]= 67; green[46]=106; blue[46]=111;
  red[47]= 70; green[47]=87; blue[47]=93;
  red[48]= 162; green[48]=197; blue[48]=198;
  red[49]= 166; green[49]=201; blue[49]=203;
  red[50]= 161; green[50]=200; blue[50]=203;
  red[51]= 153; green[51]=189; blue[51]=200;
  red[52]= 137; green[52]=175; blue[52]=185;
  red[53]= 121; green[53]=163; blue[53]=167;
  red[54]= 108; green[54]=147; blue[54]=151;
  red[55]= 106; green[55]=141; blue[55]=153;
  red[56]= 98; green[56]=131; blue[56]=161;
  red[57]= 96; green[57]=131; blue[57]=161;
  red[58]= 100; green[58]=136; blue[58]=155;
  red[59]= 111; green[59]=132; blue[59]=145;
  red[60]= 112; green[60]=134; blue[60]=144;
  red[61]= 105; green[61]=134; blue[61]=138;
  red[62]= 101; green[62]=134; blue[62]=138;
  red[63]= 105; green[63]=130; blue[63]=141;
}

void offcolors()
{
  int i;
  for (i=0; i<64; i++) {
    red[i]  = 129+2*i; 
    green[i]= 129+2*i; 
    blue[i] = 129+2*i;
  }
}

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 shape12(x,y)
     double x, y; 
{
  if (x<0) x=-x; 
  if (y<0) y=-y; 
  if ( (sqrt(2.0)*y<1.0) && (y+sqrt(3.0)*x<=sqrt(2.0)) && 
       (sqrt(2.0)*x<1.0) && (x+sqrt(3.0)*y<=sqrt(2.0)) ) 
    return 1; else return 0; 
}

int shapecircle(x,y)
     double x, y; 
{
  if ( x*x+y*y<=1)  
    return 1; else return 0; 
  
}

int chi(i)
     int i;
{if (i==0) return 0; else return 1;}
   
double sqr (x)
     double x; 
{
  return x*x; 
}

int plotstate()
{
  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 check() 
{   
 
  int i,j,iup; 
  iup=centeri+rnew+1; 
  for (i=1;i<=iup;i++){   
    for (j=1;((j<=i)&&(i+j<=nr-1));j++){
   
      if ( (apic[i][j]==1) && (adif[i][j]>0.0)) printf("*%d %lf", apic[i][j], adif[i][j]);
	   
    } 
  }
   
}

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

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

  srand48(time(NULL));
  t1=time(&t2);    
  srand48(time(NULL));
  t1=t1%1000;
  printf("seed:%d\n",t1);

  for (i=1;i<t1;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);
      }
    }       
         
  }           
}

void dynamicspop1()
{
  int i,j;
  double offset; 
  if (sigma<0){ 
    for (i=0;i<nr;i++)
      for (j=0;(j<nc);j++){
	if(apic[i][j]==0){ 
	  offset=sigma*adif[i][j]; 
	  adif[i][j]+=offset;
	}
      }
  }
}
 
 
void 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) ilo = 0;
  if(jlo<0) jlo = 0;

  // printf("ilo %d iup %d jlo %d jup %d\n",ilo,iup,jlo,jup);
  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;
	}	    	      
      }	   
    } 
  }   
}

void 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) ilo = 0;
  if(jlo<0) jlo = 0;

  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;}
}

void 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) ilo = 0;
  if(jlo<0) jlo = 0;

  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; 
	}	     
      }        
    }
  }   
}


void dynamics()
{
  int i;
  dynamicsdif();     
  dynamicsfre1();    
  dynamicsfre();         
  dynamicsunfre();   
  if (sigma>0.0) dynamicspop(); 
  //plotstate(); 
}

void picturebig()
{
  int i,j,k,pqn,kf;
   
  double y; 
  char pqc[10];
  for (i=0;i<nr;i++){
    for (j=0; j< nc; j++){
      if (apic[i][j]==0){ 
	// greys
	k = floor(63.0*(adif[i][j]/rho));                  
	XSetForeground(td,tgc,clroff[k].pixel);                 
	XFillRectangle(te.xexpose.display,
		       te.xexpose.window,
		       tgc,
		       j*sp+30,
		       i*sp+60,
		       sp,sp);
      } else {
	// blues
 	y=alm[i][j]+adif[i][j]; 
	k=floor((33.0*y-alpha)/(beta-alpha)); 
	if (k>32) k=32; 
	XSetForeground(td,tgc,clron[k].pixel);
	XFillRectangle(te.xexpose.display,
		       te.xexpose.window,
		       tgc,
		       j*sp+30,
		       i*sp+60,
		       sp,sp);
      }       
    }
  }
  if (pq==0) {
    pqc[0]='0'; 
    pqc[1]='\0';
  } else {
    pqn=pq; 
    for (k=9;k>=0;k--){
      pqc[k]=pqn%10+'0'; 
      pqn=pqn/10; 
    }
    for (kf=0;pqc[kf]=='0';kf++);
    for (k=0;k<=9-kf; k++) pqc[k]=pqc[k+kf];
    pqc[10-kf]='\0';
  }
  XSetForeground(td,tgc,tf);
  XSetBackground(td,tgc,tb);
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,10,45,timestring,strlen(timestring));
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,40,45,pqc,strlen(pqc));
} 
 
void picturerings()
{
  int i,j,k,pqn,kf;

  char pqc[10];
  
  for (i=0;i<nr;i++){
    for (j=0; j< nc; j++){
	 
      if (apic[i][j]==0){ 
                                              
	k=floor(63.0*(adif[i][j]/(rho)));                  
	XSetForeground(td,tgc,clroff[k].pixel);                 
	XFillRectangle(te.xexpose.display,te.xexpose.window,
		       tgc,j*sp+30,i*sp+60,sp,sp);
      }
      else {              
	k=ash[i][j];
	k=k%kappamax;
	XSetForeground(td,tgc,clr[k].pixel);
	XFillRectangle(te.xexpose.display,te.xexpose.window,
		       tgc,j*sp+30,i*sp+60,sp,sp); 
	if (alm[i][j]>1+0.5*(beta-1.0)){
	  if (alm[i][j]>=1+0.2*(beta-1.0)) k=12;
	  if (alm[i][j]>=1+0.5*(beta-1.0)) k=13;
	  if (alm[i][j]>=1+0.7*(beta-1.0)) k=14;
	  if (alm[i][j]>=beta) k=15; 
                   
	  XSetForeground(td,tgc,othp[k].pixel);
	  XFillRectangle(te.xexpose.display,te.xexpose.window,
			 tgc,j*sp+30,i*sp+60,sp,sp); 
	}
      }        
    }        
  }       
                    
  if (pq==0) {pqc[0]='0'; pqc[1]='\0';}
  else{
    pqn=pq; 
    for (k=9;k>=0;k--){
  
      pqc[k]=pqn%10+'0'; 
      pqn=pqn/10; 
    }
    for (kf=0;pqc[kf]=='0';kf++);
                   
    for (k=0;k<=9-kf; k++) pqc[k]=pqc[k+kf];
    pqc[10-kf]='\0';
  }
              
                    
  XSetForeground(td,tgc,tf);
  XSetBackground(td,tgc,tb);
               
  XDrawImageString(te.xexpose.display,te.xexpose.window,
		   tgc,10, 45,timestring,strlen(timestring));
  XDrawImageString(te.xexpose.display,te.xexpose.window,
		   tgc,40, 45,pqc,strlen(pqc));
   
 
} 

void drawbuttons()
{
  char quitstring[]="QUIT";
  char pausestring[]="pause";
  char playstring[]="play";
  char savestring[]="save";  
  char readstring[]="read"; 
  XSetForeground(td,tgc,tf);
  XSetBackground(td,tgc,tb);
  XDrawRectangle(te.xexpose.display,te.xexpose.window,tgc,
		 20, 50, nc*sp + 20, nr*sp+20);
  XDrawRectangle(te.xexpose.display,te.xexpose.window,tgc,
		 10, 10, 50, 20);
  XDrawRectangle(te.xexpose.display,te.xexpose.window,tgc,
		 65, 10, 50, 20);
  XDrawRectangle(te.xexpose.display,te.xexpose.window,tgc,
		 120, 10, 50, 20);
  XDrawRectangle(te.xexpose.display,te.xexpose.window,tgc,
		 175, 10, 50, 20);
  XDrawRectangle(te.xexpose.display,te.xexpose.window,tgc,
		 230, 10, 50, 20);
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,
		   20,25,quitstring,strlen(quitstring));
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,
		   75,25,pausestring,strlen(pausestring));
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,
		   130,25,playstring,strlen(playstring));
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,
		   185,25,savestring,strlen(savestring));
  XDrawImageString(te.xexpose.display,te.xexpose.window,tgc,
		   240,25,readstring,strlen(readstring));
}
 
void skip()
{
  char dum; 
  dum=getchar(); 
  while (dum!=':') dum=getchar();  
}

void readpicture()
{
  int i,j,k;
  double x;
   
  picf=fopen(infile,"r");   
  for (i=0;i<nr;i++){
    for (j=0;j<nc;j++){
      fscanf(picf,"%lf", &x);            
      adif[i][j]=x;
      fscanf(picf,"%d", &k);            
      apic[i][j]=k;
      fscanf(picf,"%lf", &x);            
      afr[i][j]=x;
      fscanf(picf,"%d", &k);
      ash[i][j]=k;
      fscanf(picf,"%lf", &x);            
      alm[i][j]=x;
    }
  } 
  fscanf(picf,"%d", &k);rold=k; 
  fscanf(picf,"%d", &k);rnew=k;
  fscanf(picf,"%d", &k);pq=k; 
  fclose(picf);
}
 
void savepicture()
{
  int i,j;
  picf=fopen(outfile, "w");
  for (i=0;i<nr;i++){
    for (j=0;j<nc;j++){
      fprintf(picf, "%.10lf %d %.10lf %d %.10lf ", 
	      adif[i][j], 
	      apic[i][j], 
	      afr[i][j], 
	      ash[i][j], 
	      alm[i][j]);
    }
  }
  fprintf(picf, "%d %d ", rold, rnew);
  fprintf(picf, "%d ", pq ); 
  fclose(picf);
}
 
void savesnowflake()
{

  int i,j,i1,j1,k,pqn,kf;

  char pqc[10]; 
  double totalmass; 
  double y; 
  
  FILE *dum;

  picf=fopen(graphicsfile, "w"); 
  fprintf(picf, "P3\n"); 
   
  fprintf(picf, "#rho:%lf\n", rho);
  fprintf(picf, "#h:%d\n",rinit);
  fprintf(picf, "#p:%lf\n", rhorinit);  
  fprintf(picf, "#beta:%lf\n", beta);
  fprintf(picf, "#alpha:%lf\n", alpha );
  fprintf(picf, "#theta:%lf\n",theta); 
  fprintf(picf, "#kappa:%lf\n", kappa);
  fprintf(picf, "#mu:%lf\n", mu);
  fprintf(picf, "#gam:%lf\n",gam);  
  fprintf(picf, "#sigma:%lf\n",sigma); 
   
  fprintf(picf, "#L:%d\n", nr); 
  fprintf(picf, "#Z:%d\n", sp); 
   
  fprintf(picf, "#: no : no : no \n"); 
   
  fprintf(picf, "#: %s\n", po);  
  fprintf(picf, "#: %s\n", comments); 
   
  fprintf(picf, "%d %d\n", nr, nr); 
  fprintf(picf, "255\n"); 
  printf("\n");
   
  for (i=0;i<nr;i++){
    for (j=0; j<nc; j++){
      i1=i;j1=j;   
      if (pq%2==1){  
	if (apic[i1][j1]==0){ 
	  k=floor(63.0*(adif[i1][j1]/(rho)));                  
	  fprintf(picf,"%d %d %d ", 
		  clroff[k].red*255/65535, clroff[k].green*255/65535,clroff[k].blue*255/65535);            
	} else {
	  y=alm[i1][j1]+adif[i1][j1]; 
	  k=floor((33.0*y-alpha)/(beta-alpha)); if (k>32) k=32; 
	  fprintf(picf,"%d %d %d ", clron[k].red*255/65535,
		  clron[k].green*255/65535,clron[k].blue*255/65535); 
	}   
      } else { 
	if (apic[i1][j1]==0){ 
	  k=floor(63.0*(adif[i1][j1]/(rho)));                  
	  fprintf(picf,"%d %d %d ", 
		  clroff[k].red*255/65535, clroff[k].green*255/65535,clroff[k].blue*255/65535);                        
	} else { 	    
	  if (alm[i1][j1]>1+0.5*(beta-1.0)){
	    if (alm[i1][j1]>=1+0.2*(beta-1.0)) k=12;
	    if (alm[i1][j1]>=1+0.5*(beta-1.0)) k=13;
	    if (alm[i1][j1]>=1+0.7*(beta-1.0)) k=14;
	    if (alm[i1][j1]>=beta) k=15; 
	    fprintf(picf,"%d %d %d ", 
		    othp[k].red*255/65535,othp[k].green*255/65535,othp[k].blue*255/65535);   
	  } else {           
	    k=ash[i1][j1];
	    k=k%kappamax;
	    fprintf(picf,"%d %d %d ", 
		    clr[k].red*255/65535,clr[k].green*255/65535,clr[k].blue*255/65535);
	  }
	}
      }	 
    }
    fprintf(picf, "\n"); 	 	 
  }	
  fclose(picf); 
}				

main(argc,argv)

int argc;
char *argv[];

{
  int i,j,k,nop;
  int reached_end;
  char ch;
  int posx,posy;

  Window rw,cw;
  int rootx,rooty;
  unsigned int kgb;

  int aflag = 0;
  int bflag = 0;
  int index;
  int c;

  rinit    = 0;
  rhorinit = 1;

  alpha = 0;     // 0 - 0.3
  beta  = 2.0;   // 1.05 - 3
  kappa = 0;     // 0 - 0.05
  gam   = 0.0000515;
  mu    = 0;     // 0 - 0.01
  rho   = 0.58;  // 0.3 - 0.9
  sigma = 0;     // -0.5 - 0
  theta = 0.002; // 0.02 - 0.57

  twelvesided = 0;
  if (rinit<0) {
    rinit =- rinit;
    twelvesided = 1;
  }

  nr = 350;
  sp = 1;
  
  while ((c = getopt(argc,argv,"a:b:k:g:m:r:s:t:n:")) != -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;
    default:
      abort();
    }
  
  printf("a %f\n",alpha);
  printf("b %f\n",beta);
  printf("k %f\n",kappa);
  printf("g %f\n",gam);
  printf("m %f\n",mu);
  printf("r %f\n",rho);
  printf("s %f\n",sigma);
  printf("t %f\n",theta);
  printf("n %d\n",nr);

  nc = nr;

  td = XOpenDisplay("");
  ts = DefaultScreen(td);
  tb = XWhitePixel(td,ts);
  tf = XBlackPixel(td,ts);
   
  th.x = 0;
  th.y = 0;
   
  th.width  = nc*sp+150;
  th.height = nr*sp+150;
   
  th.flags = PPosition | PSize;
   
  tw = XCreateSimpleWindow(td, DefaultRootWindow(td),th.x,th.y,
			   th.width, th.height,7, tf,tb);
   
  XSetWindowBorderWidth(td,tw,100);
  XSetStandardProperties(td,tw,wn,in,None,argv,argc,&th);
   
  cmap=DefaultColormap(td,ts);
   
  braquecolors64();   
  for (i=0;i<kappamax;i++){
    clr[i].red   = red[i]   *65535/255;
    clr[i].green = green[i] *65535/255;
    clr[i].blue  = blue[i]  *65535/255;
    XAllocColor(td,cmap,&clr[i]);
  }
      
  bluecolors33(); 
  // 0 is lightest blue
  for (i=0;i<=32;i++){
    clron[i].red   = red[i]   * 65535/255;   
    clron[i].green = green[i] * 65535/255;
    clron[i].blue  = blue[i]  * 65535/255;
    XAllocColor(td,cmap,&clron[i]);
  } 

  offcolors();  
  // 63 is lightest
  for (i=0;i<=63;i++){
    clroff[63-i].red   = red[i]   * 65535/255;   
    clroff[63-i].green = green[i] * 65535/255;
    clroff[63-i].blue  = blue[i]  * 65535/255;
    XAllocColor(td,cmap,&clroff[63-i]);
  }

  tgc=XCreateGC(td,tw,0,0);

  XSetBackground(td,tgc,tb);
  XSelectInput(td,tw,(ButtonPressMask|ExposureMask));
  XMapRaised(td,tw);
  fin=F;
  XNextEvent(td,&te);
  drawbuttons();
             
  printf("creating init. st.\n"); 
  initialize();
  picturebig();
  plotstate(); 

  pq=0; 

  while (fin==F){
    XNextEvent(td,&te);
    switch(te.type){
      
    case ButtonPress:
      XQueryPointer(td,tw,&rw,&cw,&rootx,&rooty,&posx,&posy,&kgb); 
      if ((posx>=10)&&(posx<=60) &&(posy>=10) && (posy<=30)){
	printf("QUIT\n");
	fin=T;
      }
      else if ((posx>=65)&&(posx<=115) &&(posy>=10) && (posy<=30)){
	printf("pause\n"); 
	if (pq%2==0) picturerings(); else picturebig(); 
      }
      else if ((posx>=120)&&(posx<=170) &&(posy>=10) && (posy<=30)){
	printf("play\n");
	while((XEventsQueued(td,QueuedAfterReading)==0)&&(pq!=-1)&&(stop==F)){
	  noac=0;
	  pq++;
	  dynamics();
	  reached_end = plotstate();
	  if(reached_end) break;
	  if (pq%25==0){
	    picturebig();
	  }
	}
      }
      else if ((posx>=175)&&(posx<=225) &&(posy>=10) && (posy<=30)){
	printf("save to file %s\n", outfile);
	savepicture();
	savesnowflake(); 
      }
      else if ((posx>=230)&&(posx<=280) &&(posy>=10) && (posy<=30)){
	printf("read from file %s\n", infile);
	readpicture(); 
	dynamicspop1();  
	picturebig();
      }
      else {
	printf("step\n");
	noac=0;
	pq++;
	dynamics();
	picturebig();
      }
          
      break;
      
    case Expose:
      if (te.xexpose.count==0){
	picturebig();
	drawbuttons();
      }
      break;
    }
  } 
  XFreeGC(td,tgc);
  XDestroyWindow(td,tw);
  XCloseDisplay(td);  
}