/* Mandelbrot set z^2+c; raster image */
/* Write postscript image */

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

#define LOG2 .69314718055994530941
#define SQRT2 1.414213562373095

#define INTM 0
#define BDRYM 1
#define RESERVED 2
#define FREE (256-RESERVED)
#define SHOWBDRY 1
#define NHD 1.0

/* Defaults */
/*Mapping z^d + c */
#define DEGREE 2

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

/* Subdivision of basin of infinity */
#define SUBDIVIDE 2.0

/* Center, radius and dimensions of image */
#define CENTER {-0.8,0.0}
#define RADIUS 1.4
#define PIXELS 200

static int degree   = DEGREE,
	   iterlim  = ITERLIM,
	   showbdry = SHOWBDRY;
static double subdivide=SUBDIVIDE, nhd = NHD;

double pixrad;

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

{
	int i,n;

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

	set_colors();
	for(i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1])
	{
		case 'b':
		showbdry = 0;
		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],"%d",&degree);
		break;

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

		case 'j':
		set_gray(1.0);
		break;

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

		case 'n':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&nhd);
		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",&subdivide);
		break;

		default:
		usage();
	}
	}
	scan_window(center,radius);
	scan_pix_radius(&pixrad);
	pixrad *= SQRT2;
	pixrad *= nhd;
	scan(color);
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  mandel [options]\n");
	fprintf(stderr,"  -b color boundary\n");
	fprintf(stderr,"  -c win center [0.0 0.0]\n");
	fprintf(stderr,"  -d raster dims [%d %d]\n",PIXELS,PIXELS);
	fprintf(stderr,"  -e deg [2]\n");
	fprintf(stderr,"  -i iter lim [%d]\n",ITERLIM);
	fprintf(stderr,"  -j :just draw bdry]\n");
	fprintf(stderr,"  -k gray :fill interior\n");
	fprintf(stderr,"  -n nhd of M [%.1lf]\n",NHD);
	fprintf(stderr,"  -r win radius [%.1lf]\n",RADIUS);
	fprintf(stderr,"  -s basin subdiv [%.0lf]\n",SUBDIVIDE);
	fprintf(stderr,"PS file is written to stdout.\n");
	exit(1);
}

set_gray(g)
        double g;
{
        int i;

        scan_gray(INTM,g);
        scan_gray(BDRYM,0.0);
        for(i=0; i<FREE; i++)
                scan_gray(RESERVED+i,1.0);
}


set_colors()
{
        int i;

        scan_rgb(INTM,255,255,255);
        scan_rgb(BDRYM,0,0,0);
        for(i=0; i<FREE; i++)
                scan_hls(RESERVED+i,
        (i*360.0/FREE),50.0,100.0);
}

color(c)
	complex c;
{
	complex z, zp, zdm1, zd;
	double d, l, logd, r2;
	int i, j;

/*	Quick tests:  main cardioid, period 2 */
	if(test1(c)) return(0);
	if(test2(c)) return(0);

	logd=log((double) degree);

	z.x = z.y = 0;    /* Critical point */
	zp.x = zp.y = 0;  /* dz_n/dc        */
	for(i=1; i<iterlim; i++)
	{	
		zdm1.x=1.0; zdm1.y=0.0;
		for(j=1; j<degree; j++) zdm1=mult(z,zdm1);
		zd=mult(z,zdm1);
		z = add(zd,c); 
/* dz_n/dc = 1 + (dz^{d-1})(dz_{n-1}/dc) */
		zdm1.x *= degree;
		zdm1.y *= degree;
		zp = mult(zp,zdm1);
		zp.x += 1;
		r2 = z.x*z.x + z.y*z.y;
		if(r2 > 64)
		{	
			l = log(r2);
			d = cx_abs(divide(z,zp)) * l;
			if(d < pixrad && showbdry) return(1);
			i = subdivide*(i - log(l)/logd);
			return(2+(254+i)%254);
		}
	}
	return(0);
}

/* Test if c is in the main cardioid */
/* Is |sqrt(1-4c)+/-1| < 1? */
test1(c)
	complex c;
{
	complex w, z;

	w.x = 1.0-4*c.x;
	w.y = 0.0-4*c.y;
	w   = cx_sqrt(w);
	z   = w;
	w.x = w.x - 1.0;
	z.x = z.x + 1.0;
	return(cx_abs(w) < 1.0 || cx_abs(z) < 1.0);
}

/* Test if c is in the period 2 component */
/* Is |c+1| < 1/4? */
test2(c)
	complex c;
{
	complex w;

	w.x = c.x+1.0;
	w.y = c.y;
	return(cx_abs(w) < 0.25);
}
