/***	ii.c - quaternion Julia set inverse iteration visualization program
 ***	compile as:
 ***			cc ii.c -o ii -lGLU -lGL -laux -lX11 -lm
 ***
 ***	plus any other flags you may need.
 ***	run as:
 ***			ii c-real c-imaginary radius
 ***
 ***    where c-real and c-imaginary is the complex constant c and radius is
 ***	a guess at the radius of the quaternion Julia set in the i-j plane
 *** 	(1.0 is usually a reasonable guess).
 ***	try: ii 0.2809 0.53 1.0
 ***
 ***	John C. Hart, April 10, 1996.
 ***	(adapted from code written for Kleiser-Walczak, 1992-1993.)
 ***/

#include <stdio.h>
#include <math.h>
#include <GL/gl.h>
#include <aux.h>

#define LENGTH 5      /* length of particle tails */

typedef struct {
        float r,i,j,k;
} Quaternion;

Quaternion qsqrt(),qsub(),qsum(),qsqr(),qmult(),snorm();
int dwell();

Quaternion c,		/* complex constant, as in z^2 + c */
	   eio;		/* e^(i theta) = cos theta + i sin theta */
float *mylast;
int addr;		/* unique deterministic address of each iteration */
int time,period,risetime,falltime,endtime,xres,yres,particles;
FILE *out = 0;		/* output filename */

main(argc,argv)
int argc;
char *argv[];
{
        Quaternion q,	/* iteration variable q |-> f(q) */
		   oc;	/* input complex constant c */
        float phi,	/* angle of point on initial circle */
	      radius,	/* radius of initial circle */
	      theta,	/* angle julia set is rotated in C */
	      rate;

	char fn[128];	/* filename */

	/***
	 *** set defaults
	 ***/

        oc.r = oc.i = oc.j = oc.k = 0.0;
        radius = 1.0;
        theta = phi = 0.0;
	time = 0;

	/***
	 *** parse command line
	 ***/

        if (argc != 4) {
                        fprintf(stderr,"%s real imag radius\n",argv[0]);
			exit(1);
        }
	oc.r = atof(argv[1]);
	oc.i = atof(argv[2]);
	radius = 1.0;
	particles = 256;

	/***
	 *** setup last array
	 ***/

	mylast = (float *)malloc(particles*LENGTH*3*sizeof(float));

	if (!mylast) {
		printf("cannot malloc last\n");
		exit(-1);
	}

	/***
	 *** set up the graphics
	 ***/

	auxInitDisplayMode(AUX_DOUBLE | AUX_RGBA);
	auxInitPosition(0,0,512,512);
	auxInitWindow(argv[0]);

	gluPerspective(90.0,1.0,0.01,4.0);
	gluLookAt(0.5,2.0,1.0,0.0,0.0,0.0,0.0,0.0,1.0);	/* j is up */

	glEnable(GL_DEPTH_TEST);

	/***
	 *** animate the point trajectories
	 ***/

	while (1) {
                eio.r = fcos(theta);
                eio.i = fsin(theta);
                eio.j = 0.0;
                eio.k = 0.0;

		/***
		 *** multiply eio*c only once at the beginning of iteration
		 *** (since q |-> sqrt(eio*(q-eio*c)))
		 ***/

                c = qmult(eio,oc);

		/***
		 *** begin with point on initial circle
		 ***/

                q.r = 0.0;
                q.i = radius*fsin(phi);
                q.j = radius*fcos(phi);
                q.k = 0.0;

		/***
		 *** iterate q and plot points
		 ***/

		glClear(GL_COLOR_BUFFER_BIT | GL_DEPTH_BUFFER_BIT);

                addr = 0;
                itertree(q,7,mylast);

		auxSwapBuffers();

		/***
		 *** change the theta increment to speed up/slow down
		 ***/

                theta += 0.01;
		time++;
        }

	while (1) sleep(60);
}

itertree(q,depth,last)
Quaternion q;
int depth;
float last[][LENGTH][3];
{
        int i,cut,uppercut;
        float len,dud;
	float x,y,z,t;

	if (time) {
		/***
		 *** move tail up
		 ***/

		for (i = LENGTH-1; i >= 1; i--) {
			last[addr][i][0] = last[addr][i-1][0];
			last[addr][i][1] = last[addr][i-1][1];
			last[addr][i][2] = last[addr][i-1][2];
		}
	}

	last[addr][0][0] = q.r;
	last[addr][0][1] = q.i;
	last[addr][0][2] = fsqrt(q.j*q.j + q.k*q.k);

	if (!time) {
		/***
		 *** initialize tail
		 ***/

		for (i = 1; i < LENGTH; i++) {
			last[addr][i][0] = last[addr][0][0];
			last[addr][i][1] = last[addr][0][1];
			last[addr][i][2] = last[addr][0][2];
		}
	}

	/***
	 *** plot
	 ***/

        glColor3f(1.0,1.0,1.0);
	glBegin(GL_LINE_STRIP);
	for (i = 0; i < LENGTH; i++) {
		glVertex3fv(last[addr][i]);
	}
	glEnd();

	/***
	 *** reflect about complex plane (fsqrt above is two valued)
	 ***/

	glBegin(GL_LINE_STRIP);
	for (i = 0; i < LENGTH; i++) {
		last[addr][i][2] = -last[addr][i][2];
		glVertex3fv(last[addr][i]);
		last[addr][i][2] = -last[addr][i][2];
	}
	glEnd();

	/***
	 *** keep addr unique
	 ***/

        addr++;

	/***
	 *** compute q = sqrt(eio*(q-eio*c))
	 *** and further iterate on q and -q
	 ***/

        if (depth) {
                q = qsqrt(qmult(eio,qsub(q,c)));
                itertree(q,depth-1,last);

                q.r = -q.r;
                q.i = -q.i;
                q.j = -q.j;
                q.k = -q.k;
                itertree(q,depth-1,last);
        }
}

Quaternion qsqrt(q)
Quaternion q;
{
        float rho,irho,srho,m,im,a,b,qjk;

        rho = fsqrt(q.r*q.r + q.i*q.i + q.j*q.j + q.k*q.k);
        irho = (rho == 0.0) ? 1000000.0 : 1.0/rho;
        q.r *= irho;
        q.i *= irho;
        q.j *= irho;
        q.k *= irho;
        m = fsqrt(q.i*q.i + q.j*q.j + q.k*q.k);
        im = (m == 0.0) ? 1000000.0 : 1.0/m;
        a = fsqrt(0.5 + 0.5*q.r);
        b = fsqrt(0.5 - 0.5*q.r);
        srho = fsqrt(rho);
        q.r = srho*a;
        q.i *= srho*b*im;
        q.j *= srho*b*im;
        q.k *= srho*b*im;

        return(q);
}

Quaternion qsub(a,b)
Quaternion a,b;
{
        a.r -= b.r;
        a.i -= b.i;
        a.j -= b.j;
        a.k -= b.k;

        return(a);
}

Quaternion qsum(a,b)
Quaternion a,b;
{
        a.r += b.r;
        a.i += b.i;
        a.j += b.j;
        a.k += b.k;

        return(a);
}

Quaternion qmult(a,b)
Quaternion a,b;
{
        Quaternion c;

        c.r = a.r*b.r - a.i*b.i - a.j*b.j - a.k*b.k;
        c.i = a.r*b.i + a.i*b.r + a.j*b.k - a.k*b.j;
        c.j = a.r*b.j + a.j*b.r + a.k*b.i - a.i*b.k;
        c.k = a.r*b.k + a.k*b.r + a.i*b.j - a.j*b.i;

        return(c);
}

Quaternion qsqr(a)
Quaternion a;
{
        Quaternion b;

        b.r = a.r*a.r - a.i*a.i - a.j*a.j - a.k*a.k;
        b.i = 2.0*a.r*a.i;
        b.j = 2.0*a.r*a.j;
        b.k = 2.0*a.r*a.k;

        return(b);
}