/*** q.c quaternion math library routines for computing quaternion Julia sets
 ***
 *** Author: John C. Hart
 *** Date: June 26, 1996 (although I've had these around for at least 6 years)
 ***
 *** Look carefully and you'll find some interesting fractal shape definitions.
 ***/

#include <math.h>

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

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

double can(double x, double y, double z)
{
	Quaternion q,q_minus_1;
	Quaternion q1,q2;
	double theta,mxy;
	int i;

	if (z < 0.5 || z > 2.5)
		return(1.0);

	q.r = 0.5;
	q.i = 0.0;
	q.j = 0.0;
	q.k = 0.0;

	mxy = x*x + y*y;

	q1.r = x*cos(3.14*z);	/* var */
	q1.i = 0;
	q1.j = y*sin(3.14*z);	/* var */
	q1.k = 0.0;

	q2.r = x;	/* var */
	q2.i = y;	/* var */
	q2.j = 0.0;
	q2.k = 0.0;

	for (i = 30; i > 0; i--) {

		q_minus_1.r = q.r - 1;
		q_minus_1.i = q.i;
		q_minus_1.j = q.j;
		q_minus_1.k = q.k;

		q = qmult(qmult(qmult(q1,q),q_minus_1),q2);

		if (q.r*q.r + q.i*q.i + q.j*q.j + q.k*q.k > 8.0)
			break;
        }

	/* printf("qjs(%g,%g,%g) = %d\n",x,y,z,i); */

	return((double)--i);
}

double turk(double x, double y, double z)
{
	double cr,ci,zr,zi,tmp;
	int i;

	cr = x;
	ci = y;
	zr = z;
	zi = 0.0;

	if (zr < 0.0) return(1.0);

	for (i = 30; i > -1; i--) {
		tmp = zr*zr - zi*zi + cr;
		zi = 2.0*zr*zi + ci;
		zr = tmp;

		if (zr*zr + zi*zi > 2.0)
			break;
	}

	return((double) i);
}


Quaternion c,		/* complex constant, as in z^2 + c */
	   eio;		/* e^(i theta) = cos theta + i sin theta */

extern double global_real, global_imag, global_theta;

double qjs(double x, double y, double z)
{
	Quaternion oc,c,eio,emio,q;
	double theta;
	int i;

	q.r = x;
	q.i = y;
	q.j = z;
	q.k = 0.0;

	oc.r = global_real;
	oc.i = global_imag;
	oc.j = oc.k = 0.0;

	theta = global_theta;

	eio.r = cos(theta);
	eio.i = sin(theta);
	eio.j = 0.0;
	eio.k = 0.0;

	emio.r = cos(-theta);
	emio.i = sin(-theta);
	emio.j = 0.0;
	emio.k = 0.0;

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

	c = qmult(eio,oc);

	for (i = 30; i > 0; i--) {
		q = qsum(qmult(emio,qsqr(q)),c);

		if (q.r*q.r + q.i*q.i + q.j*q.j + q.k*q.k > 8.0)
			break;
        }

	/* printf("qjs(%g,%g,%g) = %d\n",x,y,z,i); */

	return((double)--i);
}

void start(double *x, double *y, double *z, double inc)
{
	*x = 2.0;
	*y = 2.0;
	*z = 2.0;

	while (qjs(*x,*y,*z) > -1.0 && *x >= 0.0) {
		/* printf("%g %g %g\n",*x,*y,*z); */
		*x -= inc;
		*y -= inc;
		*z -= inc;
	}

	/* printf("starting point %g,%g,%g\n",*x,*y,*z); */
}

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

        rho = sqrt(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 = sqrt(q.i*q.i + q.j*q.j + q.k*q.k);
        im = (m == 0.0) ? 1000000.0 : 1.0/m;
        a = sqrt(0.5 + 0.5*q.r);
        b = sqrt(0.5 - 0.5*q.r);
        srho = sqrt(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);
}