Welcome to Fractal Forums

					// euclidean vector utility class with the required operator overloads
					struct triplex
					{
						double x,y,z;

						// magnitude (length) of this object
						double len() { return sqrt(x*x+y*y+z*z); }

						// return result of rotation around x axis
						triplex rotx(double arg)
						{
							triplex result = { x,
                                                                               cos(arg)*y - sin(arg)*z,
                                                                               sin(arg)*y + cos(arg)*z };
							return result;
						}

						// return result of rotation around y axis
						triplex roty(double arg)
						{
							triplex result = { cos(arg)*x - sin(arg)*z,
                                                                               y,
                                                                               sin(arg)*x + cos(arg)*z };
							return result;
						}

						// return result of rotation around z axis
						triplex rotz(double arg)
						{
							triplex result = { cos(arg)*x - sin(arg)*y,
                                                                               sin(arg)*x + cos(arg)*y,
                                                                               z                                };
							return result;
						}

						// vector addition between this object and another vector b
						triplex operator + (const triplex &b)
						{
							triplex result = {x + b.x, y + b.y, z + b.z};
							return result;
						}

						// vector division by a scalar b
						triplex operator / (const double &b)
						{
							triplex result = { x / b, y / b, z / b };
							return result;
						}

						// vector multiplication with a scalar b
						triplex operator * (const double &b)
						{
							triplex result = { x * b, y * b, z * b };
							return result;
						}

						// dot product between this object and another vector b
						double operator * (const triplex &b)
						{
							return x * b.x + y * b.y + z * b.z;
						}

						// cross product between this object and another vector b
						triplex operator % (const triplex &b)
						{
							triplex result =  { y*b.z - z*b.y, z*b.x - x*b.z, x*b.y - y*b.x };
							return result;
						}

						// raise this object to the power b relative to the pole vector
						triplex operator ^ (const double &b)
						{
							triplex Z_test = *this;

							// storing the magnitude saves us multiple calls to len()
							double magnitude = len();

							// compute normal vector that defines the plane in which the result shall live
							triplex N = (pole() % *this) / magnitude;

							// compute angle theta between current vector and the pole vector
							// which is of range 0 to PI
							double theta = acos( (*this * pole()) / magnitude );

							// generate a result vector that is situated in the plane
							// spanned by origin, pole and this object and has an angle theta
							// to the pole and is of the required magnitude.

							// rotate Z space such that normal vector points "up" (in z direction)
							double N_phi   = atan2(N.y, N.x);
							double N_theta = atan2(-N.x, N.z);
							Z_test = Z_test.rotz(-N_phi);
							Z_test = Z_test.roty(-N_theta);

							// z coord of Z_test is now expected to be 0, i.e.
							// all following operations take place in the x,y plane

							// rotate Z_test by theta*(b-1.0)
							Z_test = Z_test.rotz(theta * (b-1.0));

							// scale magnitude of Z_test with magnitude raised to the power b-1
							// so the result is the original length raised to the power b
							Z_test = Z_test * pow(magnitude, (b-1.0));

							// reverse the previous rotations
							Z_test = Z_test.roty(N_theta);
							Z_test = Z_test.rotz(N_phi);

							return Z_test;
						}

						// return the pole vector (which must have unit length)
						static triplex pole()
						{	
							triplex result = {1.0, 1.0, 1.0};
							return result / result.len();
						}
					};


// perform a loop over x,y,z here and fix how the result is stored according to your needs.

					triplex Z = {x,y,z};
					triplex C = Z;

					double mag;
					int iter = 0;

					// iteration loop. This couldn't be any easier, right?
					while((mag = Z.len()) < 2.0)
					{
						Z = (Z ^ 2.0)  + C;

						if (++iter == 100) break;
					}

					// store the result in a volume texture. This has been tweaked to result
					// in the following graphically displayed right handed coordinate system
					// x points right
					// y points "away"
					// z points up
					mem[volumeSize.width*volumeSize.height*(volumeSize.depth-1-iy) + volumeSize.width*iz + ix] = iter;

1) x^2 + y^2 + z^2 = (power * |Z|)^2
2) x*px + y*py + z*pz = power * |Z| * cos(power * theta)	
3) x*nx + y*ny + z*nz = 0												  
																						
theta is the angle between Z and the pole vector P=[px,py,pz]
N=[nx,ny,nz] is the normal vector  (normalized cross product of pole vector and Z)

[x,y,z] are then the resulting coordinates of Z^power

							// rotate Z space such that normal vector points "up" (in z direction)
							double N_phi   = atan2(N.y, N.x);
							double N_theta = atan2(-N.x, N.z);
							Z_test = Z_test.rotz(-N_phi);
							Z_test = Z_test.roty(-N_theta);

							// rotate Z space such that normal vector points "up" (in z direction)
							double N_phi   = atan2(N.y, N.x);
							triplex Nxz = N.rotz(-N_phi);                // insert this line for correct rotation
							double N_theta = atan2(-Nxz.x, Nxz.z); // and base N_theta calculation on above result
							Z_test = Z_test.rotz(-N_phi);
							Z_test = Z_test.roty(-N_theta);