partial C++ source code included. To use the code you should include outer loops to cover all voxels in x,y,z coordinates (double variables) and adjust the method by which the volume texture is stored to suit your needs
The best rendition of a mandelbrot is obtained by using a pole vector of {1.0, 0.0, 0.0}
Or maybe you're just interested in the maths.
// 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;
I would be in particular interested in a non-trigonometric approach to this, i.e. a solution to the equations
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
The first equation describes the required length of the resulting vector Z^power
The second equation describes the required angle between the pole vector and Z^power
The third equation describes the plane in which the solution needs to fall into.
Note that x,y,z are the components of the Z^power vector which we try to determine. Of the original vector Z we only
need its length |Z| which is known.
I would expect this to return two results (roots of a quadratic equation), one of which is the point I need.
Can someone help me out with Mathematica or Maple and solve this equation system for x,y,z please? I tried to do this on a sheet of paper but it got too messy and so I went for implementation of the trigonometric version first.