Welcome to Fractal Forums

TriplexNewton(QUATERNION){
parameter real Error;

quaternion z0,C,z3,cz,logc,logc1,glub,bluh,f,df;
real x1,y1,z1,r,cosz,a,b,c;

	quaternion expt(quaternion q)
	{
		x1 = part_r(q);
		y1 = part_i(q);
		z1 = part_j(q);
		
		r = exp(x1);
		cosz = cos(z1);
		a = cos(y1)*cosz;
		b = sin(y1)*cosz;
		c = sin(z1);
		
		q = quaternion(a,b,c,0)*r;
		return q;
	}
	
	quaternion lnt(quaternion q)
	{
		x1 = part_r(q);
		y1 = part_i(q);
		z1 = part_j(q);
		
		r = (x1*x1 + y1*y1 + z1*z1);
		a = ln(r)/2;
		b = atan2(x1+flip(y1));
		c = asin(z1/r);
		q = quaternion(a,b,c,0);
		return(q);
	}
	
	void init(void)
	{
		C = pixel;
		/*
			The usual algorithm starts with z=0.
			To avoid ln(0), I start with z = c/(c-1)
		*/
		logc = lnt(pixel);
		logc1 = lnt(pixel-(1,0,0,0));
		
		z = expt(logc-logc1);
	}
	void loop(void)
	{
		z0=z;
		bluh = lnt(z);
				
		z3 = expt(3*bluh);
		cz = expt(bluh+logc);
		f = z3 +cz - z - C; 
		
		df = expt(bluh*2)*3 + C -1;
		bluh = lnt(f);
		glub = lnt(df);
		z = z0 - expt(bluh-glub);
	}
	bool bailout(void)
	{
		return(|z-z0|>Error);
	}
	void description(void)
	{
		this.title="TriplexNewton";
		
		Error.default=0.001;
	}
}