Oh, my deeply sleeping thread!
I revive Thee from the dead...
Ok, seriously, after staying away from FractalForums for a while because of RealLife(TM), I came back and (obviously) found this thread. I had a little exchange of mails with twinbee about the formula involved and he has asked me to post the results here. So, here it comes:
twinbee defined (a "few" posts back)
double r = sqrt(x*x + y*y + z*z );
double yang = atan2(sqrt(x*x + y*y) , z ) // that would be theta in std polar coordinates
double zang = atan2(y , x); // that would be phi in std polar coordinates
so I would suppose he has implicitly
x = r*sin(yang)*cos(zang)
y = r*sin(yang)*sin(zang)
z = r*cos(yang)
and I would have expected (doubling the angles!)
newx = (r*r) * sin(yang*2)*cos(zang*2)
newy = (r*r) * sin(yang*2)*sin(zang*2)
newz = (r*r) * cos(yang*2)
but he defines
newx = (r*r) * sin( yang*2 + 0.5*pi ) * cos(zang*2 +pi);
newy = (r*r) * sin( yang*2 + 0.5*pi ) * sin(zang*2 +pi);
newz = (r*r) * cos( yang*2 + 0.5*pi );
which can be simplified by taking into account the symmetries of sin() and cos() to
newx = - (r*r) * cos(yang*2) * cos(zang*2)
newy = - (r*r) * cos(yang*2) * sin(zang*2)
newz = - (r*r) * sin(yang*2)
which is not exactly equal to doubling the angles.
... but it leads to interesting pictures.
One DOES NOT need atan2(), sin() and cos() to implement these formulae,
because they can be simplified a lot by using the following identities for any
angle a :
cos(a)*cos(a)+sin(a)*sin(a)=1 (well known, I suppose)
cos(2*a) = cos(a)*cos(a)-sin(a)*sin(a) (less well known, it seems
)
sin(2*a) = 2*cos(a)*sin(a)
I'll spare you the details, but you end on:
newx = ( x*x + y*y - z*z )*( x*x - y*y) / ( x*x + y*y )
newy = 2 * ( x*x + y*y - z*z )*x*y / ( x*x + y*y )
newz = - 2 * z * sqrt( x*x + y*y )
no trigonometric functions at all, just additions, multiplications, divisions and a squareroot.
There is NO pole on the z-axis, BUT there might be numerical problems because of the
denominator, solvable e.g. by taking
if( abs(y) < really_small_value )
newx = x*x-z*z
newy = 0
newz = -2*z*sqrt(x*x)
else (view above)
(end of simplification)
Interesting side effect: when I first read about "doubling both angles" I wanted
to try that for myself. When doing geometry instead of math, I prefer measuring
the angle theta not against the z-axis, but against the x-y-plane, so in that case
phi = atan2( y, x ) // just like before
theta = atan2( z, sqrt(x*x + y*y) ) // exchange of arguments, angle is positive above, negative below x-y-Plane
then
x = r*cos(theta)*cos(phi)
y = r*cos(theta)*sin(phi)
z = r*sin(theta)
and simply doubling the angles
newx = r*r*cos(2*theta)*cos(2*phi)
newy = r*r*cos(2*theta)*sin(2*phi)
newz = r*r*sin(2*theta)
simplifying this analogous to above gives
newx = ( x*x + y*y - z*z )*( x*x - y*y) / ( x*x + y*y )
newy = 2 * ( x*x + y*y - z*z )*x*y / ( x*x + y*y )
newz = 2 * z * sqrt( x*x + y*y )
IDENTICAL to twinbee's stuff except for the sign in z !
Since I'm sitting at an OLD Mac (350MHz) and try to use POV-Ray to produce pictures I can't show any yet; that takes TIME! But the first little "thumbnails" I produced show that this change of sign dramatically changes the resulting set!
Forgive my ranting - I hope somebody might find these "simplifications" useful - twinbee thought so, at least!
Happy iterating...
Karl
