//package kcdc.applets.genart;

public class Quaternion
{
	public double e = 0.0; // JM: the data
	public double i = 0.0;
	public double j = 0.0;
	public double k = 0.0;  

	public Quaternion() {}

	public String toString()
	{
		StringBuffer buffer = new StringBuffer();
		buffer.append("{");
		buffer.append(e);
		buffer.append(", ");
		buffer.append(i);
		buffer.append(", ");
		buffer.append(j);
		buffer.append(", ");
		buffer.append(k);
		buffer.append("}");
		return buffer.toString();
	}

	public double normsq()
	{ return this.e * this.e + this.i * this.i + this.j * this.j + this.k * this.k; }

	// JM: all of the following methods set "this" to the result of
	// JM: the given calculation.  

	// JM: basic quaternion ops that can have this equal to inputs

	// JM: sets this fields to a
	public void set(Quaternion a)
	{
		this.e = a.e;
		this.i = a.i;
		this.j = a.j;
		this.k = a.k;
	}

	public void qplus(Quaternion a, Quaternion b)
	{
		this.e = a.e + b.e;
		this.i = a.i + b.i;
		this.j = a.j + b.j;
		this.k = a.k + b.k;
	}

	public void qneg(Quaternion a)
	{
		this.e = -a.e;
		this.i = -a.i;
		this.j = -a.j;
		this.k = -a.k;
	}

	public void qsub(Quaternion a, Quaternion b)
	{
		this.e = a.e - b.e;
		this.i = a.i - b.i;
		this.j = a.j - b.j;
		this.k = a.k - b.k;
	}

	public void qinv(Quaternion a)
	{
		double d;

		d = a.normsq();
		if (d<=1.0e-9)
		{
			// JM: arbitrary behavior for dividing by 0 or near 0
			this.e = 0.0;
			this.i = 0.0;
			this.j = 0.0;
			this.k = 0.0;
		}
		else
		{
			d = 1.0/d;
			this.e = a.e*d;
			this.i = -a.i*d;
			this.j = -a.j*d;
			this.k = -a.k*d;
		}
	}

	public void qconj(Quaternion a)
	{
		this.e = a.e;
		this.i = -a.i;
		this.j = -a.j;
		this.k = -a.k;
	}


	// JM: Quaternion ops that can not have arguments
	// JM: equal to "this" (use scratch slots)

	public void qmult(Quaternion a, Quaternion b)
	{
		this.e = a.e*b.e - a.i*b.i - a.j*b.j - a.k*b.k;
		this.i = a.e*b.i + a.i*b.e + a.j*b.k - a.k*b.j;
		this.j = a.e*b.j - a.i*b.k + a.j*b.e + a.k*b.i;
		this.k = a.e*b.k + a.i*b.j - a.j*b.i + a.k*b.e;
	}

	public void qdiv(Quaternion a, Quaternion b, Quaternion t1)
	{
		t1.qinv(b);
		this.qmult(a,t1);
	}

	// JM: the unit Quaternions under multiplication form an important group
	// JM: ( isomorphic to SU(2) (2x2 complex matracies with determinant 1) )
	// JM: normalizing Quaternions injects us into the group
	public void qnorm(Quaternion a)
	{
		double d;

		d = a.normsq();
		if (d <= 1.0e-9)
		{
			this.e = 0.0;
			this.i = 0.0;
			this.j = 0.0;
			this.k = 0.0;
		}
		else
		{
			d = 1.0/Math.sqrt(d);
			this.e = a.e * d;
			this.i = a.i * d;
			this.j = a.j * d;
			this.k = a.k * d;
		}
	}

	// JM: the projective version of the unit Quaternions (identifying a with -a)
	// JM: is another important group isomophic to PSU(2) qnormp injects into this
	public void qnormp(Quaternion a)
	{
		this.qnorm(a);

		if (this.e!=0.0)
		{
			if (this.e<0.0)
				this.qneg(this);
			return;
		} 

		if (this.i!=0.0)
		{
			if (this.i < 0.0)
				this.qneg(this);
			return;
		} 

		if (this.j!=0.0)
		{
			if (this.j < 0.0)
				this.qneg(this);
			return;
		} 

		if (this.k!=0.0)
		{
			if (this.k<0.0)
				this.qneg(this);
			return;
		}
	}

	// JM: a coordinate independent rounding to the integer 
	// JM: subring of the Quaternions
	public void qfloor(Quaternion a)
	{
		this.e = Math.floor(a.e);
		this.i = Math.floor(a.i);
		this.j = Math.floor(a.j);
		this.k = Math.floor(a.k);
	}

	// JM: constant yielding ops

	public void qc1()
	{
		this.e = 1.0;
		this.i = 0.0;
		this.j = 0.0;
		this.k = 0.0;
	}

	public void qc2()
	{
		this.e = 0.0;
		this.i = 1.0;
		this.j = 0.0;
		this.k = 0.0;
	}

	public void qc3()
	{
		this.e = 0.0;
		this.i = 0.0;
		this.j = 1.0;
		this.k = 0.0;
	}

	public void qc4()
	{
		this.e = 0.0;
		this.i = 0.0;
		this.j = 0.0;
		this.k = 1.0;
	}

	public void qc5()
	{
		this.e = 1.61803398875;  // JM: golden ratio
		this.i = 0.0;
		this.j = 0.0;
		this.k = 0.0;
	}

	// JM: varible entry points

	public void qcx(double x, double y)
	{
		this.e = x;
		this.i = 0.0;
		this.j = 0.0;
		this.k = 0.0;
	}

	public void qcy(double x, double y)
	{
		this.e = y;
		this.i = 0.0;
		this.j = 0.0;
		this.k = 0.0;
	}

	// JM: extras to help close tree and make eqns reactive
	public void qcx1(double x, double y)
	{
		this.e = x;
		this.i = 0.0;
		this.j = 0.0;
		this.k = 1.0;
	}

	public void qcy1(double x, double y)
	{
		this.e = y;
		this.i = 0.0;
		this.j = 0.0;
		this.k = 1.0;
	}

	public void qcx2(double x, double y)
	{
		this.e = x;
		this.i = 1.0;
		this.j = 0.0;
		this.k = 0.0;
	}

	public void qcy2(double x, double y)
	{
		this.e = y;
		this.i = 0.0;
		this.j = 1.0;
		this.k = 0.0;
	}

	public void qcxy(double x, double y)
	{
		this.e = x;
		this.i = y;
		this.j = 0.0;
		this.k = 0.0;
	}

	public void qcxy2(double x, double y)
	{
		this.e = x;
		this.i = y;
		this.j = x;
		this.k = y;
	}

	// JM: coordinate independent functions
	// JM: not interesting from the Quaternion point of view- but since we
	// JM: derive colors from the coordinates independently these functions should
	// JM: have some visual value.

	public void qisin(Quaternion a)
	{
		this.e = Math.sin(a.e);
		this.i = Math.sin(a.i);
		this.j = Math.sin(a.j);
		this.k = Math.sin(a.k);
	}

	public void qilog(Quaternion a)
	{
		this.e = Math.log(Math.max(a.e,0.000001));
		this.i = Math.log(Math.max(a.i,0.000001));
		this.j = Math.log(Math.max(a.j,0.000001));
		this.k = Math.log(Math.max(a.k,0.000001));
	}

	public void qiexp(Quaternion a)
	{
		this.e = Math.exp(Math.min(Math.max(a.e,-30.0),30.0));
		this.i = Math.exp(Math.min(Math.max(a.i,-30.0),30.0));
		this.j = Math.exp(Math.min(Math.max(a.j,-30.0),30.0));
		this.k = Math.exp(Math.min(Math.max(a.k,-30.0),30.0));
	}

	public void qimin(Quaternion a, Quaternion b)
	{
		this.e = Math.min(a.e,b.e);
		this.i = Math.min(a.i,b.i);
		this.j = Math.min(a.j,b.j);
		this.k = Math.min(a.k,b.k);
	}

	public void qimax(Quaternion a, Quaternion b)
	{
		this.e = Math.max(a.e,b.e);
		this.i = Math.max(a.i,b.i);
		this.j = Math.max(a.j,b.j);
		this.k = Math.max(a.k,b.k);
	}


	// JM: shift coordinates around (not a Quaternion automorphism)

	public void qrl(Quaternion a)
	{
		double t;
		t = a.e;
		this.e = a.i;
		this.i = a.j;
		this.j = a.k;
		this.k = t;
	}

	public void qrr(Quaternion a)
	{
		double t;
		t = a.k;
		this.k = a.j;
		this.j = a.i;
		this.i = a.e;
		this.e = t;
	}


	// JM: more exotic ops (can not have this equal to an argument)

	// JM: all R-algrebra automorphisms of H (the Quaternions) are of the form
	// JM: qaut1 by corollary of a theorem of Cayley's
	public void qaut1(Quaternion a, Quaternion b, Quaternion t1, Quaternion t2)
	{
		t2.qdiv(b,a,t1);
		this.qmult(a,t2);
	}

	// JM: not an automorphism- but close
	public void qaut2(Quaternion a, Quaternion b, Quaternion t1, Quaternion t2, Quaternion t3)
	{
		t3.qconj(b);
		this.qaut1(a,t3,t1,t2);
	}

	// JM: exponential map -
	// JM:    using first few terms of power series as aproximation
	public void qexp(Quaternion ai, Quaternion p, Quaternion a, Quaternion t3) // JM: temps
	{
		a.set(ai);
		// JM: don't touch big arguments
		while(a.normsq()>=900.0)
		{
			a.e /= 10.0;
			a.i /= 10.0;
			a.j /= 10.0;
			a.k /= 10.0;
		}

		this.set(a);
		p.set(a);
		this.e += 1.0; // JM: zero power term

		for (int ti=2;ti<10;++ti)
		{
			double d;
			d = 1.0/(double)ti;
			t3.qmult(p,a);
			p.e = t3.e*d;
			p.i = t3.i*d;
			p.j = t3.j*d;
			p.k = t3.k*d;
			t3.set(this);
			this.qplus(t3,p);
		}
	}


	// JM: there is a polar repersentation of Quaternions: for every Quaternion a
	// JM: there is an imaginvar Quaternion u s.t. a = |a| exp(u).
	// JM: finding u looks hard.
	// JM: we can use the fact the u.u = -|u|^2 to break exp(u) into two real
	// JM: series:  sum_{k=0}^{\infty} (-|u|^2)^k / (2 k)! +
	// JM: u (sum_{k=0}^{\infty} (-|u|^2)^k / ( 2 k +1)!
	// JM: but this looks like a mess


	// JM: by Hamilton's theorem for every ortogonal mappinger f: Im H->Im H there
	// JM: is a unit Quaternion a s.t. f(b) = a b bar(a) or f(b) = - a b bar(a)
	// JM: qorth1 and qorth2 implement these maps
	public void qorth1(Quaternion a, Quaternion b, Quaternion t1, Quaternion t2, Quaternion t3)
	{
		t1.qnorm(a);
		t2.qconj(t1);
		t3.qmult(b,t2);
		this.qmult(t1,t3);
	}

	public void qorth2(Quaternion a, Quaternion b,
	Quaternion t1, Quaternion t2, Quaternion t3)
	{
		this.qorth1(a,b,t1,t2,t3);
		this.qneg(this);
	}

	// JM: mod (or remainder) over the Quaternions in analogy to traditional mod
	public void qmod(Quaternion a,Quaternion b, Quaternion t1)
	{
		this.qdiv(a,b,t1);
		this.qfloor(this);
		t1.qmult(this,b);
		this.qsub(a,t1);
	}
}