
/* Bifurcation plot of z^2+c, c real */

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

/* Picture frame [0,1]x[0,1] */
#define CENTER {0.5,0.5}
#define RADIUS 0.5
#define IDIM 1152
#define JDIM 900

/* Test for divergence */
#define LARGE 4.0

/* Number of iterations*/
#define ITERLIM 1000  

/* Range of parameters to examine */
#define CLOW -2.0
#define CHIGH -1.0
#define XLOW -2.1
#define XHIGH 2.1

static int iterlim=ITERLIM;
static double clow  = CLOW;
static double chigh = CHIGH;
static double xlow  = XLOW;
static double xhigh = XHIGH;
static int idim=IDIM;
static int jdim=JDIM;

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

	static complex center = CENTER;
	static int exponents = 0;
	double radius = RADIUS;
	int color(); 

	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",&clow);
		sscanf(argv[++i],"%lf",&chigh);
		break;

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

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

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

		default:
		usage();
	}
	}
	scan_gray(0,1.0);
	scan_gray(1,0.0);
	scan_window(center,radius);
        scan_dim(idim,jdim);
	scan_begin();
	plot_attr();
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  bif [options]\n");
	fprintf(stderr,"  -c [clow chigh]\n");
        fprintf(stderr,"  -d [n m] nxm raster\n");
	fprintf(stderr,"  -i [iteration limit]\n");
	fprintf(stderr,"  -x [xlow xhigh]\n");
	fprintf(stderr,"Postscript file is written to stdout.\n");
	exit(1);
}

plot_attr()
{	
	int i, j;
	double c, x;
	complex z;

	c = clow;
	do
	{	 z.y = (c-clow)/(chigh-clow);
		c += (chigh-clow)/(2*jdim);
		x = 0.0;
		for(i=0; i<iterlim; i++) x = x*x + c;
		for(i=0; i<iterlim; i++)
		{	x = x*x + c;
			z.x = (x-xlow)/(xhigh-xlow); 
			scan_set(z,1);
		}
	} while (z.y < 1.0);
}
