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

lcca.c					Eric Weeks, 12-7-98

12-07-98: started from pamp2ppm.c, to get color scheme, PPM code, etc

Eric R. Weeks
email:                weeks@physics.emory.edu
home page:            http://www.physics.emory.edu/~weeks/
link to this program: http://www.physics.emory.edu/~weeks/pics/lcca1.html

This program is public domain, but (1) please leave my name and
contact information attached (in this header) and (2) I made no
attempts to make this user-friendly or document it.  It is a fairly
simple program, though, so you might find it useful.

To compile:    
     gcc -o lcca lcca.c -lm

This program formatted for a tab stop of 5

 *******************************************************************/
#include <stdio.h>
#include <math.h>
#include <stdlib.h>
#include <time.h>

#define PI 3.1415926535897932384626433833
#define MAX 1005

float array[MAX][MAX];
int nx,ny,blackwhite=0,truncate=0,coloring=0;
int stretchv;
float stretchf;

void printcolor(float stream,int flag);
void rand_init();





float pot(float theta1,float theta2)  {
	float dt;
	dt = ((theta1>theta2) ? (theta1-theta2) : (theta2-theta1));
	dt = ((dt > 0.5) ? (dt-1.0) : (dt));
	/* fprintf(stderr,"%f \n",dt); */
	return(dt*dt);
}

float funct(int x,int y)
{
	float newtheta,theta,v1,v2,v3,vmin,dt,t1,t2,t3;
	int x1,y1,x2,y2,t;

	x1 = ((x==0) ? (nx-1) : (x-1));
	x2 = ((x==(nx-1)) ? (0) : (x+1));
	y1 = ((y==0) ? (ny-1) : (y-1));
	y2 = ((y==(ny-1)) ? (0) : (y+1));
	vmin = 99999.9;
	for (newtheta=0.0;newtheta<1.0;newtheta+=0.1)  {
		v1=pot(array[x1][y],newtheta)+pot(array[x2][y],newtheta);
		v1+=pot(array[x][y1],newtheta)+pot(array[x][y2],newtheta);
		v1+=pot(array[x1][y1],newtheta)+pot(array[x1][y2],newtheta);
		v1+=pot(array[x2][y1],newtheta)+pot(array[x2][y2],newtheta);
		if (v1 < vmin)  {
			vmin=v1;
			t1=newtheta;
		}
	}
	for (newtheta=t1-0.1;newtheta<t1+0.1;newtheta+=0.005)  {
		v1=pot(array[x1][y],newtheta)+pot(array[x2][y],newtheta);
		v1+=pot(array[x][y1],newtheta)+pot(array[x][y2],newtheta);
		v1+=pot(array[x1][y1],newtheta)+pot(array[x1][y2],newtheta);
		v1+=pot(array[x2][y1],newtheta)+pot(array[x2][y2],newtheta);
		if (v1 < vmin)  {
			vmin=v1;
			t1=newtheta;
		}
	}
	if (truncate==1) { t1=t1+0.08; }

	return(t1);
}


main(argc,argv)
int argc;
char **argv;
{
	int t,r,n,x,y,maxt;
	int red,green,blue;
	int c;					/* for command-line argument */
	int atoi();
	extern int optind;
	extern char *optarg;
	int n1,n2,n3,n4,n5,expand;
 
	nx=100;
	ny=100;
	maxt=1;
	expand=1;
	stretchf=1.0;
	stretchv=255;
	while ((c = getopt(argc, argv, "hx:y:n:Te:Bc:s:S:")) != EOF)
        switch (c)
        {
	case 'h':
fprintf(stderr,"cellular automata, simulates liquid crystals\n");
fprintf(stderr,"Usage:  lcca [opts] > out.ppm\n");
fprintf(stderr,"\nOptions:\n");
fprintf(stderr,"  -h    : this help message\n");
fprintf(stderr,"  -x #  : change x-size [%d]\n",nx);
fprintf(stderr,"  -y #  : change y-size [%d]\n",ny);
fprintf(stderr,"  -n #  : change number of iterations [%d]\n",maxt);
fprintf(stderr,"  -e #  : expand factor [%d]\n",expand);
fprintf(stderr,"  -T    : truncate color\n");
fprintf(stderr,"  -B    : black/white rather than color\n");
fprintf(stderr,"  -c #  : coloring option\n");
fprintf(stderr,"  -s #  : stretch colors towards this value [%d]\n",stretchv);
fprintf(stderr,"  -S #  : stretch factor [%.1f] (<= 1)\n",stretchf);
			exit(1);
			break;
	case 'x': nx = atoi(optarg);
			if (nx>MAX)  nx = MAX;
			break;
	case 'y': ny = atoi(optarg);
			if (ny>MAX)  ny = MAX;
			break;
	case 'n': maxt=atoi(optarg);
			break;
	case 'T': truncate=1;
			break;
	case 'e': expand=atoi(optarg);
			break;
	case 'B': blackwhite=1;
			break;
	case 'c': coloring=atoi(optarg);
			break;
	case 's': stretchv=atoi(optarg);
			break;
	case 'S': stretchf=atof(optarg);
			break;
        }

	/* argc -= (optind-1) ; argv += (optind-1) ; */
	/* fp = (argc > 1) ? fopen(*++argv, "r") : stdin; */

	rand_init();
	for (x=0;x<nx;x++)  {
		for (y=0;y<nx;y++)  {
			array[x][y] = drand48();
		}
	}
	fprintf(stderr,"done initializing\n");

	x = 0; y = 0;
	for (n1=0;n1<2;n1++)  {
		fprintf(stderr,"%d / 2\n",n1);
		for (t=0;t<maxt;t++)  {
			for (r=0;r<60000;r++) {
				array[x][y]=funct(x,y);
				x = x + 7;
				if (x >= nx)  {
					x -= nx;
					y = y + 1;
					if (y >= ny)  {
						y-=ny;
					}
				}
			}
		}
		if ((n1==0)&&(expand>1))  {
			for (n2=nx-1;n2>=0;n2--)  {
				for (n3=ny-1;n3>=0;n3--)  {
					for (n4=0;n4<expand;n4++)  {
						for (n5=0;n5<expand;n5++)  {
							array[n2*expand+n4][n3*expand+n5]=array[n2][n3];
						}
					}
				}
			}
			nx*=expand; ny*=expand;
		}
	}


 
	if (blackwhite==0)  printf("P6\n");
	else  printf("P5\n");
	printf("# created with program lcca (ERW)\n");
	printf("%d %d\n255\n",nx,ny);
	for (n2=0;n2<ny;n2++)  {
		for (n1=0;n1<nx;n1++)  {
			printcolor(array[n1][n2],1);
		}
	}

}
/* ------------------------------------------------------------ */
/* ------------------------------------------------------------ */
  

void printcolor(float datum,int flag)
{
	int red,green,blue;
	float x,y;


	if (flag==0)  {
		/* WHITE */
		red = 255; green = 255; blue = 255;
	} else if (flag>=2)  {
		/* BLACK */
		red = 0; green = 0; blue = 0;
	} else {
		if (blackwhite==0)  {
			x = datum*2.0;
			while (x>1.0)  x -= 2.0;
			while (x< (-1.0))  x += 2.0;
			/***** now -1 < x < +1 *****/
			if (coloring==3)  {
				y=(y+1)*2; if (y>2.0) y=4.0-y;
				y-=1.0;
			}
			if (x<-.9)  {
				/* DK RED to RED */
				y = (x+1.0)/.1;
				red = 255;
				green = (0.2+(0.5*y-0.2))*255;
				blue =  (0.2+(0.5*y-0.2))*255;
			} else if (x<-.5)  {
				/* RED to ORANGE */
				y = (x+0.9)/.4;
				red = 255;  green = (0.5+(0.67-0.5)*y)*255;
				blue = (0.5+(-0.5)*y)*255;
			} else if (x<-.1) {
				/* ORANGE to YELLOW */
				y = (x+0.5)/.4;
				red = 255;  green = (0.67+(1.0-0.67)*y)*255;
				blue = 0;
			} else if (x<0.2)  {
				/* YELLOW to LT GREEN */
				y = (x+0.1)/.3;
				red = (1.0+(0.67-1.0)*y)*255;  green = 255;
				blue = 0;
			} else if (x<.4)  {
				/* LT GREEN to BLUE */
				y = (x-.2)/.2;
				red = (.67+(.5-.67)*y)*255;  green = 255;
				blue = y*255;
			} else if (x<.7)  {
				/* BLUE to VIOLET */
				y = (x-.4)/.3;
				red = 128;  green = (1.0+(0.5-1.0)*y)*255;
				blue = 255;
			} else if (x<0.9)  {
				/* VIOLET to DK PURPLE */
				y = (x-0.7)/0.2;
				red =   (0.5+(0.8-0.5)*y)*255;
				green = (0.5+(0.2-0.5)*y)*255;
				blue =  255;
			} else {
				/* DK PURPLE to RED */
				y = (x-0.9)/0.1;
				red =   (0.8+(1.0-0.8)*y)*255;
				green = (0.2+(0.2-0.2)*y)*255;
				blue =  (1.0+(0.2-1.0)*y)*255;
			}
		} else {
			y = datum;
			if (coloring==1)  {
				y*=2; if (y>1.0) y=2.0-y;
			}
			if (coloring==2)  {
				if (y<0.0) y+=1.0;
				if (y>1.0) y-=1.0;
				y*=2; if (y>1.0) y=2.0-y;
			}
			if (coloring==3)  {
				if (y<0.0) y+=1.0;
				if (y>1.0) y-=1.0;
				y*=2; if (y>1.0) y=2.0-y;
				y*=2; if (y>1.0) y=2.0-y;
			}
			if (coloring==4)  {
				if (y<0.0) y+=1.0;
				if (y>1.0) y-=1.0;
				y*=2; if (y>1.0) y=2.0-y;
				y*=2; if (y>1.0) y=2.0-y;
				y*=2; if (y>1.0) y=2.0-y;
			}
			red = (y)*255;
			green = (y)*255;
			blue = (y)*255;
		}
	}

	red = (int)(((float)red-stretchv)*stretchf+stretchv);
	green = (int)(((float)green-stretchv)*stretchf+stretchv);
	blue = (int)(((float)blue-stretchv)*stretchf+stretchv);

	if (blackwhite==0)  {
		fputc((char)red,stdout);
		fputc((char)green,stdout);
		fputc((char)blue,stdout);
	} else fputc((char)red,stdout);
}


/* initialize random number generator */
void rand_init()  {
        int seed;
        time_t timenow;
        time(&timenow);
        seed=timenow % 100000;
        srand48(seed);
}