/* Plots limit sets of Kleinian groups */

#include <math.h>
#include <stdio.h>
#include <stdlib.h>
#include "cx.h"
#include "circles.h"
#include "lim.h"

/* Defaults */
static double agree=AGREE, eps=EPS, huge=HUGENESS;
static int depth = DEPTH, mode=PLANE, verbose=1;
static int nfinal=0, ngens=0, nstack=0, maxc=MAXC;
static size_t circle_bytes;

circle final[MAXC], stack[MAXC];
matrix gens[MAXG];

/* Main program */
int main(argc,argv)
	int argc;
	char *argv[];
{
	int i;
	complex c;
	double gray, lw, r;
	double xmin, ymin, xmax, ymax;

	circle_bytes = sizeof(circle);

	for (i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1]) 
	{
	case 'a':
		if(argc <= ++i) usage();
		sscanf(argv[i],"%lf",&agree);
		break;
	case 'b':
		ps_box();
		break;
	case 'c':
		if(argc <= (i+3)) usage();
		sscanf(argv[++i],"%lf",&c.x);
		sscanf(argv[++i],"%lf",&c.y);
		sscanf(argv[++i],"%lf",&r);
		ps_clip_circle(c,r);
		break;
	case 'd':
		if(argc <= ++i) usage();
		sscanf(argv[i],"%d",&depth);
		break;
	case 'e':
		if(argc <= ++i) usage();
		sscanf(argv[i],"%lf",&eps);
		break;
	case 'h':
		if(argc <= ++i) usage();
		sscanf(argv[i],"%lf",&huge);
		break;
	case 'l':
		if(argc <= ++i) usage();
		sscanf(argv[i],"%lf",&lw);
		set_linewidth(lw);
		break;
	case 'm':
		if(argc <= ++i) usage();
		sscanf(argv[i],"%d",&maxc);
		if(maxc>MAXC)
		{	maxc=MAXC;
			fprintf(stderr,
	"Maximum number of circles is %d\n",maxc);
		}
		break;
	case 'q':
		verbose = 0;
		break;
	case 's':
		mode=SPHERE;
		break;
	case 't':
		mode=SEETHRU;
		if(argc <= ++i) usage();
		sscanf(argv[i],"%lf",&gray);
		break;
	case 'w':
		if(argc <= (i+4)) usage();
		sscanf(argv[++i],"%lf",&xmin);
		sscanf(argv[++i],"%lf",&ymin);
		sscanf(argv[++i],"%lf",&xmax);
		sscanf(argv[++i],"%lf",&ymax);
		ps_win(xmin,ymin,xmax,ymax);
		break;
	default:
		usage();
	}
	}

	read_stack();
	generate();
	if(verbose) fprintf(stderr,
		"Final circle count: %d\n",nfinal); 

	ps_open(argc,argv);
	switch(mode)
	{	case PLANE:
		ps_plot_circles(final,nfinal);
		break;

		case SPHERE:
		ps_plot_circles_s2(final,nfinal);
		break;

		case SEETHRU:
		ps_see_thru(final,nfinal,gray);
	}
	ps_close();
}

void usage()
{
	fprintf(stderr,"Limit sets of Kleinian groups\n");
        fprintf(stderr,"Usage:  lim [options] < circles > postscript\n");
        fprintf(stderr,"Op [Argument]              (default)\n");
        fprintf(stderr,"-a [agreement]             (%d)\n",AGREE);
        fprintf(stderr,"-b boxing                  (none)\n");
        fprintf(stderr,"-c [c.x c.y r] clip circle (none)\n");
        fprintf(stderr,"-d [depth]                 (%d)\n",DEPTH);
        fprintf(stderr,"-e [eps  cutoff radius]    (%.3lf)\n",EPS);
        fprintf(stderr,"-h [huge cutoff radius]    (%.0lf)\n",HUGENESS);
        fprintf(stderr,"-l [linewidth]             (%.2lf)\n",LINEWIDTH);
        fprintf(stderr,"-m [max no. of circles]    (%d)\n",MAXC);
        fprintf(stderr,"-q quieter\n");
        fprintf(stderr,"-s plot on sphere          (planar)\n");
        fprintf(stderr,"-t [gray, 0-1.0] transparent sphere]\n");
        fprintf(stderr,"-w [xmin ymin xmax ymax] window]\n");
        exit(1);
}

void read_stack()
{
	int i;
	matrix t, ti;
	circle circ;
	char code[20];

	nstack = ngens = 0;
	while(EOF != scanf("%s",code))
	switch(code[0])
	{	
	case 'c':
		scanf("%lf %lf %lf",
			&circ.c.x,&circ.c.y,
			&circ.r);
		push_circle(circ);
		break;
	case 'm':
		scanf("%lf %lf %lf %lf %lf %lf %lf %lf",
			&gens[ngens].a.x, &gens[ngens].a.y,
			&gens[ngens].b.x, &gens[ngens].b.y,
			&gens[ngens].c.x, &gens[ngens].c.y,
			&gens[ngens].d.x, &gens[ngens].d.y
			);
		gens[ngens]=make_sl2(gens[ngens]);
		ngens++;
		break;
	case 'r':
		scanf("%lf %lf %lf",
			&circ.c.x,&circ.c.y,
			&circ.r);
		gens[ngens] = 
			circle_to_matrix(circ);
		ngens++;
		break;
	case 't':
		scanf("%lf %lf %lf %lf %lf %lf %lf %lf",
			&t.a.x,&t.a.y,&t.b.x,&t.b.y,
			&t.c.x,&t.c.y,&t.d.x,&t.d.y
			);
		t = make_sl2(t);
		conjugate_by(t);
		break;
	case 'u':
		scanf("%lf %lf %lf",
			&circ.c.x,&circ.c.y,
			&circ.r);
		t = circle_to_matrix(circ);
		conjugate_by(t);
		break;
	default:
		fprintf(stderr,"Unrecognized code: %s\n",code);
	}
}

void conjugate_by(t)
	matrix t;
{
	matrix ti;
	int i;

	ti = inverse(t);
	for(i=0; i<nstack; i++)
		stack[i] = image_circle(t,stack[i]);
	for(i=0; i<ngens; i++)
		gens[i] = product(t,product(gens[i],ti));
}

/* Return code: 0) OK 1) too small or far 2) full */
int push_circle(circ)
	circle circ;
{
	int i;

	if(circ.r < eps && circ.r > 0) return(1);
	if(circ.r > huge || cx_abs(circ.c) > huge) return(1);
	if(nstack==maxc) return(2);
	stack[nstack]=circ;
	nstack++;
	return(0);
}

void evolve(c)
	circle c;
{
	circle d;
	int i;

	for(i=0; i<ngens; i++)
	{	d = image_circle(gens[i],c);
		push_circle(d);
	}
}

void evolve_stack()
{
	int i, nold;

	nold = nstack;
	for(i=0; i<nold; i++) evolve(stack[i]);
}

/* Sorting routines */

/* Cut circle of slopes at an irrational number such as PI */
double slope(x)
	double x;
{
	x = x + PI;
	while (x < 0)    x = x + 1.0; 
	while (x >= 1.0) x = x - 1.0;
	return(x);
}
	
int compare_circles(c,d)
	circle *c, *d;
{
	double dr, dx, dy;
	double small;

	/* Lines precede circles */
	if(c->r <= 0 && d->r > 0) return(-1); 
	if(d->r <= 0 && c->r > 0) return(1); 

	small=eps/agree;
	dx = c->c.x-d->c.x;
	if(dx < -small) return(-1);
	if(dx > small) return(1);

	dy = c->c.y-d->c.y;
	if(dy < -small) return(-1);
	if(dy > small) return(1);

	if(c->r <= 0) 		/* Lines */
	{	dr = slope(c->r)-slope(d->r);
		if(dr < -PARAFUZZ) return(-1);
		if(dr >  PARAFUZZ) return(1);
		return(0);
	} else {		/* Circles */
		dr = c->r-d->r;
		if(dr < -small) return(-1);
		if(dr > small) return(1);
		return(0);
	}
}

void sort_circles(stack,nstack)
	circle stack[];
	int nstack;
{
	qsort(stack,nstack,circle_bytes,compare_circles);
}

/* Sort a list of circles and eliminate dups */

void unique_circles(stack,nstack)
	circle stack[];
	int *nstack;
{
	int i, n;

	if(*nstack<2) return;
	sort_circles(stack,*nstack);
	i=n=1;
	while(i<*nstack)
	{	if(0==compare_circles(&stack[i],&stack[i-1]))	i++; 
	else {stack[n]=stack[i]; n++; i++;}
	}
	*nstack = n;
}

/* Move circles from stack to final, skipping dups */
/* Returns 1 if run out of room */
int merge_circles()
{
	int c, i, ifinal, istack, nold;

	unique_circles(stack,&nstack);
	if(nstack+nfinal > maxc) return(1);

	ifinal = istack = 0;
	nold=nfinal;
	while(ifinal < nold && istack < nstack)
	{	c =compare_circles(&final[ifinal],&stack[istack]);
		if(c<0) ifinal++;
		if(c==0) istack++;
		if(c>0)
		{	final[nfinal]=stack[istack];
			nfinal++;
			istack++;
		}
	}
	while(istack < nstack)
	{	final[nfinal]=stack[istack];
		nfinal++;
		istack++;
	}
	for(i=nold; i<nfinal; i++)
		stack[i-nold]=final[i];
	nstack = nfinal-nold;
	sort_circles(final,nfinal);  
	return(0);
}

void generate()
{
	int i, nevo;

	if(verbose) fprintf(stderr,
	"%6d Circles/ %6d Generators\n", nstack,ngens);

/* Put original circles on final list */
	for(i=0; i<nstack; i++) final[i]=stack[i]; 
	nfinal = nstack;
	unique_circles(final,&nfinal);

	for(i=1; i<=depth; i++)
	{	if(nstack==0)
		{	if(verbose)
			fprintf(stderr,"Generation complete.\n");
			break;
		}

		nevo = nstack;
		evolve_stack();

		if(verbose) fprintf(stderr,
		"Gen %3d: final %6d evolving %6d stack %6d\n",
			i,nfinal,nevo,nstack);
		if(merge_circles())
		{	fprintf(stderr,"Stack full\n");
			break;
		}
	}
	if(verbose && i > depth)
		fprintf(stderr,"Depth attained\n");
}
