/* Julia set for renormalizations of z^2+c;
   postscript image */

/* Specify iterate f^n to renormalize, and a
   circle around the origin to define a
   quadratic-like map.
*/


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

/* Renormalization f^p */
#define P 2		

/* Renormalization radius */
#define RENORMR 1.2

/*Maximum number of iterations*/
#define ITERLIM 20  
#define LOG2 .69314718055994530941
#define JFAC 0.01
#define SUBDIVIDE 8

/* Default map is Feigenbaum */
#define C {-1.401155189,0}

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

/* Number of points in a circle */
#define PTS 4000

/* Maximum renormalization (f^MAX) */
#define MAX 1000

/* Number of colors to use */
#define CIRCLE 0
#define JULIA 1
#define FILLED_JULIA 2
#define OUTSIDE 3
#define RESERVED 4
#define FREE (256-RESERVED)

/* Color map */
#define CMAP COLOR
#define COLOR 0
#define GRAY 1
#define BW 2

static complex c = C;
static int iterlim=ITERLIM, mit=1;
static int cmap = CMAP, draw=1, p=P;
static double renormr = RENORMR;
static double jfac = JFAC, subdivide = SUBDIVIDE;
static double g = 1.0;

double pixrad;

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

{
	int i,n;

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

	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",&c.x);
		sscanf(argv[++i],"%lf",&c.y);
		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 'f':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&jfac);
		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 'k':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&g);
		cmap = BW;
		break;

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

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

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

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

		case 'u':
		draw=0;
		break;

		default:
		usage();
	}
	}
	scan_window(center,radius);
	scan_pix_radius(&pixrad);
	scan_begin();
	scan(color);
	if(draw) draw_circles();
	set_map();
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  renorm [options]\n");
	fprintf(stderr,"  -C c.x c.y [for z^2 + c]\n");
	fprintf(stderr,"  -c x y [center of window]\n");
	fprintf(stderr,"  -d n m [nxm raste]r\n");
	fprintf(stderr,"  -f f [Julia fuzz factor]\n");
	fprintf(stderr,"  -g [gray tones]\n");
	fprintf(stderr,"  -i n [iteration limit]\n");
	fprintf(stderr,"  -j [Julia set]\n");
	fprintf(stderr,"  -k f [filled Julia set]\n");
	fprintf(stderr,"  -p n [renormalize f^n]\n");
	fprintf(stderr,"  -r f [radius of window]\n");
	fprintf(stderr,"  -R f [renormalization radius]\n");
	fprintf(stderr,"  -s f [subdivision factor]\n");
	fprintf(stderr,"  -u [don't draw circles]\n");
	fprintf(stderr,"Postscript file is written to stdout.\n");
	exit(1);
}

color(z)
	complex z;
{
	int k, n;
	complex fz, fpz;

	n = escape(z,&fz,&fpz);
	if(n == 0) return(OUTSIDE);
	if(n < 0) return(FILLED_JULIA); 
	if(cx_abs(fz)/cx_abs(fpz) < jfac*pixrad) return(JULIA);
	k = RESERVED + (n-1)%FREE;
	if(n>mit) mit = n;
	return(k);
}

escape(z,fz,fpz)
	complex z;
	complex *fz, *fpz;
{
	int i, j;

	*fz = z;
	fpz->x = 1.0; fpz->y = 0.0;
	for(i=0; i<iterlim; i++)
	{	if(cx_abs(*fz) > renormr) return(i);
		for(j=0; j<p; j++) 
		{	*fz = add(mult(*fz,*fz),c); 
			*fpz = mult(*fpz,*fz);
			fpz->x *= 2.0; fpz->y *= 2.0;
		}
	}
	return(-1);
}

draw_circles()
{
	complex corb[MAX], zorb[MAX], z;
	double t;
	int i, j;

	corb[0].x=corb[0].y=0;
	for(i=0; i<p; i++) corb[i+1] =
		add(mult(corb[i],corb[i]),c);

	zorb[p].x = renormr;
	zorb[p].y = 0;
	scan_set(zorb[p],CIRCLE);
	for(i=p-1; i>=0; i--)
		zorb[i] = contsqrt(sub(zorb[i+1],c),
			corb[i]);
	scan_set(zorb[0],CIRCLE);

	for(i=1; i<2*PTS; i++)
	{	t = 2*PIVAL*i/PTS;	
		zorb[p] = polar(renormr,t);
		scan_set(zorb[p],CIRCLE);
		for(j=p-1; j>=0; j--)
			zorb[j] =
		contsqrt(sub(zorb[j+1], c), zorb[j]);
		scan_set(zorb[0],CIRCLE);
	}
}

set_map()
{
	if(cmap == COLOR) set_colors();
	if(cmap == GRAY) set_grays();
	if(cmap == BW) set_bw();
}

set_bw()
{
        int i, j, k;

        scan_gray(JULIA,0.0);
        scan_gray(FILLED_JULIA,g);
        scan_gray(OUTSIDE,1.0);
        scan_gray(CIRCLE,0.0);
	for(i=0; i<FREE; i++)
        	scan_gray(RESERVED+i,1.0);
}

set_grays()
{
        int i, j, k;
	double g;

        scan_gray(JULIA,0.0);
        scan_gray(FILLED_JULIA,0.01);
        scan_gray(CIRCLE,0.0);
        scan_gray(OUTSIDE,1.0);
	for(i=0; i<mit; i++)
	{	k = RESERVED+i;
		g = i;
		g = 0.5+0.5*g/mit;
		scan_gray(k,g);
	}
}

set_colors()
{
        int i, j, k;
	double h, l, s;

        scan_rgb(JULIA,0,0,0);
        scan_hls(FILLED_JULIA,0.0,70.0,0.0);
        scan_rgb(CIRCLE,0,0,0);
        scan_rgb(OUTSIDE,255,255,255);
	for(i=0; i<mit; i++)
	{	k = RESERVED+i;
		h = i*360.0/mit;
		l = 50.0;
		s = 100.0;
		scan_hls(k,h,l,s);
	}
}
