END OF AN ERA, FRACTALFORUMS.COM IS CONTINUED ON FRACTALFORUMS.ORG

rho2  = sqrt(x**2+y**2)
rho8 = ((rho2-2)**2+z**2)**4
phi    = 8*ATAN2(y,x)
theta = 8*ATAN2(z,rho2-2.)
x = (rho8*COS(theta)+2) * COS(phi) + x_0
y = (rho8*COS(theta)+2) * SIN(phi) + y_0
z = rho8*SIN(theta)                 + z_0

R=2
ni=0
RHO = sqrt(x^2+y^2)
rho = sqrt((RHO-R)^2+z^2)
phi    = ATAN2(y,x)
theta = ATAN2(z,RHO-R)
Cx   = rho*COS(theta) * COS(phi)  //This Cx,Cy,Cz adjustment is new. Now the fractal set has a volume, a well defined orbit trap.
Cy   = rho*COS(theta) * SIN(phi)
Cz   = rho*SIN(theta)

Iterate, while rho<2R+2
------------------
RHO = sqrt(x^2+y^2)
rho = ((RHO-R)^2+z^2)^4
phi   = 8.*ATAN2(y,x)
theta = 8.*ATAN2(z,RHO-R)
x = (rho*COS(theta)+R) * COS(phi) + Cx
y = (rho*COS(theta)+R) * SIN(phi) + Cy
z =  rho*SIN(theta)               + Cz
ni++
------------------
I didn't manage to derivate this, so I used:
step=min(max(1E-8; 0.05/ni/ni*0.5**ni); 0.05)
There hasn't been done much tweaking of this formula, so it's probably pretty inefficient.