Clarification of the bugman generalised formulas

[ZsquaredplusC]

I have a quick and hopefully easy question.

I have got the trig formulas working, but am having problems with the non-trig formulas here
http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/msg8680/#msg8680

z[0,1,2]=z[x,y,z]
c[0,1,2]=complex "+c"

Currently using the trig formulas...
    z[0]=power(R,pow)*cos(pow*ph)*cos(pow*th)+c[0]
    z[1]=power(R,pow)*cos(pow*ph)*sin(pow*th)+c[1]
    z[2]=power(R,pow)*sin(pow*ph)+c[2]

But when I use the non trig (the following is for positive power 2)...
     rxy=sqrt(z[0]*z[0]+z[1]*z[1])
     a=1-z[2]*z[2]/rxy*rxy

     ztmp[0]=(z[0]*z[0]-z[1]*z[1])*a
     ztmp[1]=2*z[0]*z[1]*a
     ztmp[2]=2*z[2]*rxy

     z[0]=ztmp[0]+c[0];
     z[1]=ztmp[1]+c[1];
     z[2]=ztmp[2]+c[2];

I get no results.  And the iteration loop seems to hang like it is stuck waiting for the bounding volume to be exceeded or the epsilon closeness value to be triggered?
Is there something obviously wrong here?

Z starts at 1,0,0

Any ideas?  It would be great to get these speedups and save the trig calculations

[ZsquaredplusC]

OK, maybe I wasn't clear or my code wasn't clear.  I will try again.  Please help if you have code working using Paul's (bugman) non trig versions.

The usual trigonometric formulas...

zx=power(R,pow)*cos(pow*ph)*cos(pow*th)+cx
zy=power(R,pow)*cos(pow*ph)*sin(pow*th)+cy
zz=power(R,pow)*sin(pow*ph)+cz

Now, when trying to replace the above trig versions with Paul's non-trig versions from
http://www.fractalforums.com/3d-fractal-generation/true-3d-mandlebrot-type-fractal/msg8680/#msg8680
I get the following for power +2

In the mathematica output (sorry I cannot do latex) he shows...
rxy=sqrt(sqr(zx)+sqr(zy))
and then for the power 2 version the rest is
{x,y,z}^2={{x^2-y^2}a,2xya,2zrxy},a=1-z^2/rxy^2

So I am assuming (please correct me if I am wrong) that you first calculate rxy and a and then feed them through the formula for power 2.

This is how I attempted in the first post of this thread with the code;

Firstly work out the values for rxy and a...
rxy=sqrt(zx*zx+zy*zy)
a=1-sqr(zy)/sqr(rxy)

Then use a temp varialable to track the new xyz...
zxtmp=(sqr(zx)-sqr(zy))*a
zytmp=2*zx*zy*a
zztmp=2*zz*rxy

Then assign the new xyz values plus the C constant...
zx=zxtmp+cx
zy=zytmp+cy
zz=zztmp+cz

And this is also for the Z=Z^pow+C calculations right?  The distance estimation calculations still need to use the log etc?

Am I interpreting the mathematica output correctly?  The only change I have made from working trig functions to non-working non-trig formulas is as above.

Thanks for any insight. 

[David Makin]

I get the following for power +2

In the mathematica output (sorry I cannot do latex) he shows...
rxy=sqrt(sqr(zx)+sqr(zy))
and then for the power 2 version the rest is
{x,y,z}^2={{x^2-y^2}a,2xya,2zrxy},a=1-z^2/rxy^2

So I am assuming (please correct me if I am wrong) that you first calculate rxy and a and then feed them through the formula for power 2.

This is how I attempted in the first post of this thread with the code;

Firstly work out the values for rxy and a...
rxy=sqrt(zx*zx+zy*zy)
a=1-sqr(zy)/sqr(rxy)

Then use a temp varialable to track the new xyz...
zxtmp=(sqr(zx)-sqr(zy))*a
zytmp=2*zx*zy*a
zztmp=2*zz*rxy

Then assign the new xyz values plus the C constant...
zx=zxtmp+cx
zy=zytmp+cy
zz=zztmp+cz
Thanks for any insight. 

You've used sqr(zy) instead of sqr(zz) in "a=1-sqr(zy)/sqr(rxy)"

[ZsquaredplusC]

Thanks David.  That error was just in my typing.  It is the correct sqr(zz) as in the first post.
Must be some other error in the my surrounding code then.

[David Makin]

Thanks David.  That error was just in my typing.  It is the correct sqr(zz) as in the first post.
Must be some other error in the my surrounding code then.

In case it helps at all here's my UF iteration (for the -sine version):

            zri = sqr(zri)*(1.0 - sqr(zj)/(r=|zri|)) + cri
            zj = -2.0*zj*sqrt(r) + cj

zri and cri are complex and |zri| is x^2+y^2 where zri = (x+i*y), the other terms are reals.

[David Makin]

I just realised that you didn't mention testing for zero.

Assuming you calculated the square of the magnitude for bailout purposes then try the following:

Before the iteration loop precalculate the magnitude^2 (i.e. x^2+y^2+z^2)

At the beginning of the iteration loop if the magnitude^2 is zero then set the value to the constant i.e. z[0]=c[0], z[1]=c[1],z[2]=c[2] otherwise do the calculation as normal.
Obviously at the end of the loop calculate the new value of the magnitude^2 which is then used for bailout and at the start of the next loop.