
/* Julia set for z^2+c; raster image */
/* Write to stdout a Postscript file */
/* Default is color image */

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

/* CAUTION is 2^{1/n} */
#define LOG2 .69314718055994530941
#define SQRT2 1.414213562373095
#define LARGE 200.0
#define LENGTH 200
#define CAUTION 0.9330329915368074
#define CIRCLEPTS 1000
#define EPS 0.000001
#define RAYEPS 1e-20
#define JFAC 1.0

/* Defaults */
#define C {0, 1}
#define CENTER {0.0,0.0}
#define RADIUS 1.5
#define PCLIM 0
#define SUBDIVIDE 8
#define ITERLIM 50  

/* Colors */
#define FILLED_JULIA 0
#define JULIA 1
#define PC_SET 2 
#define RAY 3
#define LEVEL 4
#define RESERVED 5
#define FREE (256-RESERVED)

/* Rational numbers */
typedef struct {int num, den;} rational;

static int iterlim=ITERLIM, pclim=PCLIM, subdivide=SUBDIVIDE;
static complex c = C;
double pixrad;
static complex center = CENTER;
double radius = RADIUS;
static double jfac = JFAC;

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

{
	int i,n;

	static int nlev=0, nang = 0;
	double alev[100], lev[100];
	rational ang[100];

	int color(), idim, jdim;
	double gray;

	/* Default color map */
	set_colors(); 	
	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 '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 'l':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&lev[nlev]);
		nlev++;
		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",&subdivide);
		break;

		case 'x':
		if(argc <= (i+3)) usage();
		sscanf(argv[++i],"%d",&ang[nang].num);
		sscanf(argv[++i],"%d",&ang[nang].den);
		sscanf(argv[++i],"%lf",&alev[nang]);
		nang++;
		break;

		default:
		usage();
	}
	}
	scan_window(center,radius);
	scan_pix_radius(&pixrad);
	scan(color); 
	if(pclim > 0) draw_pc(); 
	for(i=0; i<nang; i++)
		draw_ray(alev[i], ang[i]);
	for(i=0; i<nlev; i++)
		draw_level(lev[i]);
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  quad [options]\n");
	fprintf(stderr,"  -C [c.x c.y] for z^2 + c]\n");
	fprintf(stderr,"  -c [center of window]\n");
	fprintf(stderr,"  -d [n m] nxm raster\n");
	fprintf(stderr,"  -f [(default 1.0) nhd of J to draw]\n");
	fprintf(stderr,"  -i [iteration limit]\n");
	fprintf(stderr,"  -j [just draw Julia set]\n");
	fprintf(stderr,"  -k gray [fill Julia set]\n");
	fprintf(stderr,"  -l [level] - image of |z|=2^level under\n");
	fprintf(stderr,"     Riemann mapping\n");
	fprintf(stderr,"  -p [n] show n points of post-critical orbit\n");
	fprintf(stderr,"  -s subdivision of basin of infinity\n");
	fprintf(stderr,"  -r [radius of window]\n");
	fprintf(stderr,"  -x [ang-num ang-denom level] - image of external ray\n");
	fprintf(stderr,"     at rational ang in [0,1], from level down\n");
	fprintf(stderr,"Postscript file is written to stdout.\n");
	exit(1);
}

/* F : \cx-K(f) -> {Re z > 0} is the multivalued
   function given by lim 2^{-n} log f^n(z).
   Re F(z) = rate, 
   (Re F(z))/|F'(z)| = dist.
   Returns 1 if z escaped, 0 otherwise.
*/

escape(z,c,iterlim,rate,dist)
	complex z,c;
	int iterlim;
	double *rate, *dist;
{
	int i;
	double fp, r, s;

	s = 1.0;
	fp = 1.0;
	for(i=0; i<iterlim; i++)
	{	r = cx_abs(z);
		if(r > LARGE)
		{	*rate = log(r)/s;
			*dist = r*log (r)/fp;
			return(1);
		}
		fp = fp * 2.0 * r;
		z = add(mult(z,z),c);
		s = s * 2.0;
	}
	return(0);
}

/* Analytically continue Riemann mapping */
/* F : unit disk -> \cx-K(f) */
/* Assuming F(r,ang) is approx. z, find exact value*/

complex riemann(r,ang,z)
	double r;
	rational ang;
	complex z;
{
	complex orb[LENGTH], neworb[LENGTH];
	int i, n;

	orb[0] = z;
	n = 0;
	for (i=0; i<LENGTH-1; i++)
	{	if(cx_abs(orb[i]) > LARGE) break;
		orb[i+1] =
		add(mult(orb[i],orb[i]),c);
		n++;
		r = 2*r;
		ang.num = (2*ang.num)%ang.den;
	}
	neworb[n] = polar(pow(2.0,r),2*PIVAL*ang.num/ang.den);

	for(i=n; i>0; i--)
	{	neworb[i-1] =
		contsqrt(sub(neworb[i],c),orb[i-1]);
	}
	return(neworb[0]);
}

set_gray(g)
	double g;
{
	int i;

	scan_gray(FILLED_JULIA,g);
	scan_gray(JULIA,0.0);
	scan_gray(PC_SET,0.0);
	scan_gray(RAY,0.0);
	scan_gray(LEVEL,0.0);
	for(i=0; i<FREE; i++)
		scan_gray(RESERVED+i,1.0);
}

set_colors()
{
	int i;

	scan_hls(FILLED_JULIA,0.0,70.0,0.0);
	scan_rgb(JULIA,0,0,0);
	scan_rgb(PC_SET,0,255,0);
	scan_rgb(RAY,0,0,0);
	scan_rgb(LEVEL,0,0,0);
	for(i=0; i<FREE; i++)
		scan_hls(RESERVED+i,
	((FREE-i)*360.0/FREE),50.0,100.0);
}

color(z)
	complex z;
{
	complex grad;
	double rate, dist;
	int i;

	if(escape(z,c,iterlim,&rate,&dist))
	{	if(dist < jfac*SQRT2*pixrad) 
			return(JULIA);
		i = subdivide * log(rate)/LOG2;
		while (i<0) i += FREE;
		return(RESERVED + i%FREE);
	}
	return(FILLED_JULIA);
}

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

	z.x = z.y = 0;			/* Critical point */

	for(i=1; i<pclim; i++)
	{	scan_set(z,PC_SET);
		z = add(mult(z,z),c); 
	}
}

/* Draw an external ray */
draw_ray(level,ang)
	double level;
	rational ang;
{
	complex lastz, z;
	double r;

	r = level;
	while(pow(2.0,r) < LARGE) r *= 2;
	z = polar(pow(2.0,r),2*PIVAL*ang.num/ang.den);
	do
	{	r = r*CAUTION;
		lastz = z;
		z = riemann(r,ang,lastz);
		if(r < level*0.99) 
			scan_line(lastz,z,RAY);
	} while(r > RAYEPS);
}


/* Color a level curve of escape */
draw_level(level)
	double level;
{
	double r;
	complex lastz, w, z;
	static rational zero={0,1};
	rational ang;

	r = level;
	while(pow(2.0,r) < LARGE) r *= 2;
	z = polar(pow(2.0,r),0.0);
	do
	{	r = r * CAUTION;
		lastz = z;
		z = riemann(r,zero,lastz);
	} while(r > level);

	lastz = z;
	ang.den = CIRCLEPTS;
	for(ang.num=0; ang.num<=CIRCLEPTS; ang.num++)
	{	z = riemann(level,ang,lastz);
		if(ang.num>0) scan_line(z,lastz,LEVEL);
		lastz = z;
	}
}

/* Set the color on a line joining z to w */
scan_line(z,w,col)
	complex z, w;
	int col;
{	
	int i, n;
	complex u;

	n = 4*cx_abs(sub(z,w))/pixrad;
	if(n == 0)
	{	scan_set(z,col);
		return;
	}
	for(i=0; i<=n; i++)
	{	u.x = (i*z.x + (n-i)*w.x)/n; 	
		u.y = (i*z.y + (n-i)*w.y)/n; 	
		scan_set(u,col);
	}
}

