/* Julia set in the plane */

/* Color scheme:  
  0- did not converge;
  1+nit*i+n -converged to ith attractor in n iterations;
  nit = (COLORS-1)/(number of attractors).
*/

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

/* Colors */
#define JULIA 0
#define PC_SET 1
#define RESERVED 2
#define FREE (256-RESERVED)

/* Color map data */
#define COLOR 0
#define GRAY 1
#define BW 2
#define CMAP COLOR
#define SPREAD 0  /* Tones in basins */
#define ROTATE 0  /* Attractor numbering */
#define PCLIM 0   /* Iterates of pc set to draw */

/* Test for convergence */
#define EPS 0.01

/*Maximum number of iterations*/
#define ITERLIM 50  

/* Picture center and radius */
#define CENTER {0.0,0.0}
#define RADIUS 5.0

/* Attractor */
complex at[100];
int na;

/* No. of iter. to distinguish per attractor */
/* Max which actually occurred */
int nit;
static int mit = 1;
static int iterlim=ITERLIM, pclim = PCLIM;

/* Distance from attractor */
static double eps = EPS;

/* Map to iterate */
complex iter(), (*f)();

/* Color info */
static int cmap = CMAP, rotate = ROTATE;
static double spread = SPREAD;

#include "scheme.c"

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

{
	int i;

	static complex center = CENTER;
	double radius = RADIUS;
	int color(), idim, jdim;

/* Default function */
pick_map("P3_4",&f);

	for(i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1])
	{
		case 'a':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&rotate);
		break;

		case 'c':
		if(argc <= (i+2)) usage();
		sscanf(argv[++i],"%lf",&center.x);
		sscanf(argv[++i],"%lf",&center.y);
		break;

		case 'd':
		if(argc <= (i+2)) usage();
		sscanf(argv[++i],"%d",&idim);
		sscanf(argv[++i],"%d",&jdim);
		scan_dim(idim,jdim);
		break;

		case 'e':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&eps);
		break;

		case 'f':
		if(arg_map(argc,argv,&i))
			usage();
		break;

		case 'g':
		cmap = GRAY;
		break;

		case 'i':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&iterlim);
		break;

		case 'j':
		cmap = BW;
		break;

		case 'p':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%d",&pclim);
		break;

		case 'r':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&radius);
		break;

		case 's':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&spread);
		break;

		default:
		usage();
	}
	}
	attractor(at,&na);
	nit = FREE/na;
	scan_window(center,radius);
	scan(color);
	if(pclim > 0) draw_pc();
	set_map();
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  plane [options]\n");
	fnusage();
fprintf(stderr,"  -a i [rotate attractor list by i\n");
fprintf(stderr,"  -g [grayscale image]\n");
fprintf(stderr,"  -c x y [center of window]\n");
fprintf(stderr,"  -d n m [nxm image]\n");
fprintf(stderr,"  -e [eps] test for convergence\n");
fprintf(stderr,"  -i i [iteration limit]\n");
fprintf(stderr,"  -j just draw Julia set\n");
fprintf(stderr,"  -p iters [postcritical set]\n");
fprintf(stderr,"  -r r [window radius]\n");
fprintf(stderr,"  -s s [basin spread; 0-1.0]\n");
fprintf(stderr,"Postscript file is written to stdout.\n");
	exit(1);
}

color(z)
	complex z;
{
	int i,j,k;

	for(i=0; i<iterlim; i++)
	{	z = iter(z);
		for(j=0; j<na; j++)
	if(cx_abs(sub(z,at[j])) < eps)
	{	i = i%nit;
		k=RESERVED+nit*j+i;
		if(i>=mit) mit = i+1;	
		return(k);
	}
	}
	return(JULIA);
}

/* Color the post-critical set */
draw_pc()
{
        complex z;
        int i, j;
	complex *cps;
	int ncps;

	critical_points(&cps,&ncps);

        for(i=0; i<ncps; i++)
	{	z = cps[i];
	        for(j=0; j<pclim; j++)
	        {       scan_set(z,PC_SET);
	                z = iter(z);
	        }
	}
}

