
/* Bifurcation set in a family */
/* We study the behavior of cp[0] */

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

/* Test for convergence */
#define EPS 0.01

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

/* Picture center and radius */
#define CENTER {-0.75,0.0}
#define RADIUS 2.5

/* Colors */
#define JULIA 0
#define PC_SET 0  /* Just for compatibliity */
#define RESERVED 1
#define FREE (256-RESERVED)
#define COLOR 0
#define GRAY 1
#define BW 2
#define CMAP COLOR
#define SPREAD 0

/* Color scheme and some of its variables */
int na, nit;
static int mit=1;
static int cmap=CMAP, rotate=0;
static double spread=SPREAD;
#include "scheme.c"

static int iterlim=ITERLIM;
static char name[20]="POLY_2";
static complex param1 = {1.0,0.0}, param2 = {-1.0,0.0};
static double eps=EPS;
complex iter();

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

{
	int i;

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

/* Default function */
	for(i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1])
	{
		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(argc <= (i+1)) usage();
		sscanf(argv[++i],"%s",name);
		break;

		case 'g':
		cmap = GRAY;
		break;

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

		case 'p':
		if(argc <= (i+2)) usage();
		sscanf(argv[++i],"%lf",&param2.x);
		sscanf(argv[++i],"%lf",&param2.y);
		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();
	}
	}
	scan_window(center,radius);
	set_att_iter(iterlim);
	scan(color);
	set_map();
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  family [options]\n");
	fnusage();
	fprintf(stderr,"  -c [center of window]\n");
	fprintf(stderr,"  -d [n m] nxm raster\n");
	fprintf(stderr,"  -e [eps] test for convergence\n");
	fprintf(stderr,"  -g black and white picture\n");
	fprintf(stderr,"  -f [name] name of function\n");
	fprintf(stderr,"  -i [iteration limit]\n");
	fprintf(stderr,"  -p [set extra parameter for function]\n");
	fprintf(stderr,"  -r [radius of window]\n");
	fprintf(stderr,"  -s [basin spread; 0-1.0]\n");
	fprintf(stderr,"Raster file is written to stdout.\n");
	exit(1);
}

/* Color point according to the least i>0 such that c0 and ci */
/* have the same attractor.  Otherwise return black (JULIA) */
color(z)
	complex z;
{
	int i,j,k,nc;
	complex at[100], *cps, w;

	pick_map(name,z,param2);
	attractor(at,&na);
	critical_points(&cps,&nc);
	na = na - 1;

	nit = FREE/na;
	w = cps[0];

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