Welcome to Fractal Forums

triplex tadd(triplex a, triplex b) {
	triplex n;
	n.x = a.x + b.x;
	n.y = a.y + b.y;
	n.z = a.z + b.z;
	return n;
}
triplex tsub(triplex a, triplex b) {
	triplex n;
	n.x = a.x - b.x;
	n.y = a.y - b.y;
	n.z = a.z - b.z;
	return n;
}
triplex tdiv(triplex a, double b) {
	triplex n;
	n.x = a.x / b;
	n.y = a.y / b;
	n.z = a.z / b;
	return n;
}
triplex tmul(triplex a, double b) {
	triplex n;
	n.x = a.x * b;
	n.y = a.y * b;
	n.z = a.z * b;
	return n;
}


triplex tpow(triplex a, double p) {
	triplex n;

	double r     = sqrt(a.x*a.x + a.y*a.y + a.z*a.z) ;
	double phi   = atan2( pow(a.x*a.x + a.y*a.y, 0.5) , a.z) ;
	double theta = atan2(a.y , a.x) ;

	r=pow(r,p);	
	phi = phi*p;
	theta = theta * p;

	n.x = r * sin(phi) * cos(theta);
	n.y = r * sin(phi) * sin(theta);
	n.z = r * cos(phi);
	return n;
}

triplex tpow(triplex a, double p) {

	... same ...

	double phi   = atan2( a.z , pow(a.x*a.x + a.y*a.y, 0.5) ) ;

	... same ...

	n.x = r * cos(phi) * cos(theta);
	n.y = r * cos(phi) * sin(theta);
	n.z = r * sin(phi);

	... same ...

}

triplex tpowPOSZ(triplex &a, double p) {
	triplex n;

    	double r    = sqrt(a.x*a.x + a.y*a.y + a.z*a.z) ;
	double theta = atan2(a.y , a.x) ;
	double phi = atan2( a.z , sqrt(a.x*a.x + a.y*a.y) ) ;

	r=pow(r,p);	phi *=p;	theta *=p;

	n.x = r * cos(theta) * cos(phi) ;
	n.y = r * cos(phi) * sin(theta);
	n.z = r * sin(phi);
	return n;
}

Ztanzbulb(QUATERNION) {
parameter real bailout;
parameter int steps;
real pixelr,pixeli,pixelj;
real x1,x2,y1,y2,z1,z2,rho1,rho2,A,B,C;
quaternion res, res1,c;
int m;

quaternion inv(quaternion q)
{
	q = q-quaternion(0,2*part_i(q),0,0);
	q = q/|q|;
	return(q);
}

quaternion mult(quaternion q, quaternion p)
{
	x1 = part_r(q);
	y1 = part_i(q);
	z1 = part_j(q);
	rho1 = sqrt(x1*x1+y1*y1);
	
	x2 = part_r(p);
	y2 = part_i(p);
	z2 = part_j(p);
	rho2 = sqrt(x2*x2+y2*y2);
	
	if (rho1 == 0){
		p =(z1*p);
	} else if (rho2 == 0){
		p = z2*p;
	} else{
		A = (x1*x2-y1*y2)*(z1/rho1+z2/rho2);
		B = (x1*y2+x2*y1)*(z1/rho1+z2/rho2);
		C = (rho1*rho2-z1*z2);
		p = (quaternion(A,B,C,0));
	}
	return(p);
}
quaternion div(quaternion q, quaternion p)
{
	x1 = part_r(q);
	y1 = part_i(q);
	z1 = part_j(q);
	rho1 = sqrt(x1*x1+y1*y1);
	
	x2 = part_r(p);
	y2 = part_i(p);
	z2 = part_j(p);
	rho2 = sqrt(x2*x2+y2*y2);
	
	if (rho1 == 0){
		p =inv(p)*z1;
	} else if (rho2 == 0){
		p = p/z2;
	} else{
		A = (x1*x2+y1*y2)*(z1/rho1-z2/rho2)/|p|;
		B = (x1*y2-x2*y1)*(z1/rho1-z2/rho2)/|p|;
		C = (rho1*rho2+z1*z2)/|p|;
		p = (quaternion(A,B,C,0));
	}
	return(p);
}

quaternion tangent(quaternion q)
{
	res = 0.0;
	for (m = steps+1; m>1; m=m-1){
		res1 = mult(z,z)*(-1.0/(2*m-3)/(2*m-1));
		res = div(res1,res+quaternion(0,0,1,0));
	}
	return(div(z,res+quaternion(0,0,1,0)));
}

	void init(void)
	{
		pixelr=part_r(pixel);
		pixeli=part_i(pixel);
		pixelj=part_j(pixel);
		c=quaternion(pixelr,pixeli,pixelj,0);
		z=c;
	}
	void loop(void)
	{
		z=mult(z,tangent(z))+c;
	}
	bool bailout(void)
	{
		return(|z|<bailout);
	}

	void description(void)
	{
		this.title = "ztanzbulb";

		bailout.caption = "Bailout Value";
		bailout.default = 4.0;
		bailout.min = 1.0;
		bailout.hint = "Defines the bailout radius: As soon as a pixel falls outside a circle with this radius, the iteration stops.";
		
		steps.caption = "Steps";
		steps.default = 2;
		steps.min = 1;
		steps.hint = "Controls how deep into the continued fraction the tangent function goes.";
		
	}
}