/* Julia set on the sphere */

/* Color scheme:  
  0- outside the circle;
  1- bounding circle;
  2- points which do not converge;
  3+nit*i+n -converged to ith attractor in n iterations;
  nit = FREE/(number of attractors).

*/

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

/* Colors */
#define SPACE 0
#define CIRCLE 1
#define JULIA 2
#define PC_SET 2 /* Just for compatibility */
#define RESERVED 3
#define FREE (256-RESERVED)
#define COLOR 0
#define GRAY 1
#define BW 2
#define CMAP COLOR
#define SPREAD 0  /* Tones in basins */
#define ROTATE 0  /* Attractor numbering */

/* Color info */
static int cmap = CMAP, rotate = ROTATE;
static double spread = SPREAD;
static int mit=1;


/* Test for convergence */
#define EPS 0.01

/* Number of points on the bounding circle */
#define NCIRCLE 4000

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

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

/* Function */
complex (*f)(), iter();

/* Attractor */
complex at[100]; 
int na, nit;

/* Parameter */
complex param;

/* Twist picture by this Mobius transformation */
static complex a = {1,0}, b = {0,0},
               c = {0,0}, d = {1,0};

static int iterlim=ITERLIM;
static double eps=EPS;

#include "scheme.c"

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

{
	int fun, i;

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

/* Default function */
	pick_map("S3_4");

	for(i=1; i<argc; i++)
	{	if(argv[i][0] != '-') usage();
		switch(argv[i][1])
	{

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

                case 'f':
                if(arg_map(argc,argv,&i)) usage();
                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 'r':
		if(argc <= (i+1)) usage();
		sscanf(argv[++i],"%lf",&radius);
		break;
                 
                case 's':
                if(argc <= (i+1)) usage();
                sscanf(argv[++i],"%lf",&spread);
                break;

		case 't':
		if(argc <= (i+8)) usage();
		sscanf(argv[++i],"%lf",&a.x);
		sscanf(argv[++i],"%lf",&a.y);
		sscanf(argv[++i],"%lf",&b.x);
		sscanf(argv[++i],"%lf",&b.y);
		sscanf(argv[++i],"%lf",&c.x);
		sscanf(argv[++i],"%lf",&c.y);
		sscanf(argv[++i],"%lf",&d.x);
		sscanf(argv[++i],"%lf",&d.y);
		break;

		default:
		usage();
	}
	}
	attractor(at,&na);
	nit = FREE/na;
	scan_window(center,radius);
	scan(color);
	scan_rgb(SPACE,255,255,255);
	set_map();
	draw_circle();
	scan_end(argc,argv);
}

usage()
{
	fprintf(stderr,"Usage:  sphere [options]\n");
	fnusage();
	fprintf(stderr,"  -a [i] rotate color map]\n");
	fprintf(stderr,"  -c [center of window]\n");
	fprintf(stderr,"  -d [n m] nxm raster\n");
	fprintf(stderr,"  -e [eps] test for convergence\n");
	fprintf(stderr,"  -i [iteration limit]\n");
	fprintf(stderr,"  -j just draw Julia set\n");
	fprintf(stderr,"  -r [radius of window]\n");
	fprintf(stderr,"  -s [0-1.0] tones in basins]\n");
        fprintf(stderr,"  -t [a.x a.y ... d.y] twist by {{a,b},{c,d}}\n");
	fprintf(stderr,"Raster file is written to stdout.\n");
	exit(1);
}

color(z)
	complex z;
{
	int i,j,k;

	if(cx_abs(z) >= 1.0) return(SPACE);
	z = mobius(a,b,c,d,disk_to_sphere(z));

	for(i=0; i<iterlim; i++)
	{	z = iter(z);
		for(j=0; j<na; j++)
		if(cx_abs(sub(z,at[j])) < eps)
		{	i = i % nit;
			k=RESERVED+nit*j+i;
			if(i>=mit) mit=i+1;
			return(k);
		}
	}
	return(JULIA);
}

/* Color the circle */
draw_circle()
{
	complex z;
	int i;

	for(i=0; i<NCIRCLE; i++)
	{	z=polar(1.0,i*2*PIVAL/NCIRCLE);
		scan_set(z,1);
	}
}
