
/* Compute uniformizing map to exterior of the Mandelbrot set */

/* Complex arithmetic package */
#include "cx.h"
#include <stdlib.h>
#include <stdio.h>

/* 	Max orbit size */
/*	Escape radius  */
/* 	Accuracy goal */
/*	Caution (in [0,1]) for updating c */
#define 	R  	1000000
#define 	EPS 	0.0000000001
#define		N 	1000
#define		CAUTION 0.99

complex orbit[N], c;
int norb;

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

	startup();

	if(argc<2) usage();
	for(i=1; i<argc; i++)
	{	if(argc <= (i+1)) usage();
		sscanf(argv[i],"%lf",&t.x);
		sscanf(argv[i+1],"%lf",&t.y);
		i = i+1;

		locate(t);
		printf("%12.8lf %12.8lf\n",c.x,c.y);
	}
}

usage()
{
	printf("Usage:  m t1.x t1.y [t2.x t2.y ...]\n");
	printf("Maps points t, |t|>1, to the outside of the Mandelbrot set\n");
	exit(1);
}

complex sq(z)
{
	return(mult(z,z));
}

complex f(z)
{
	return(add(sq(z),c));
}
   
/* Set orbit[N-i]=f_c^{-i} sq^{N}(z) */
complex preorbit(z)
	complex z;
{
	int i;

	init(z);
	for(i=1; i<norb; i++) z = sq(z);
	orbit[norb-1]=z;
	for(i=norb-1; i>0; i--)
		orbit[i-1]=
		contsqrt(sub(orbit[i],c),orbit[i-1]);
	return(orbit[0]);
}

/* Adjust orbit (larger or smaller) */
init(z)
	complex z;
{
	int i;

	if(norb==0) {orbit[0]=z; norb=1;}

	for(i=1; i<N; i++)
	{	z=orbit[i-1];
		if(cx_abs(z) > R) {norb=i; return;}
		if(i>=norb) {orbit[i]=f(z); norb=i+1;}
	}
	printf("Orbit overflow; limit is %d\n",N);
	exit(1);
}

/* Find value of c */
locate(t)
	complex t;
{
	complex cnew, cold;

	do
	{	cold=c;
		cnew=preorbit(t);
		c.x = (1-CAUTION)*cnew.x + CAUTION*cold.x;
		c.y = (1-CAUTION)*cnew.y + CAUTION*cold.y;
	}
	while (cx_abs(sub(c,cold)) > EPS);
}

startup()
{
	norb=0;
	c.x = c.y = 0;
}
