Let us imagine a parametric function X=f(r,phi,theta)

(A surface at a threshold r and two angles phi, theta.)

Let there be G = [r,theta,phi] = g(X) = f^-1(X), either explicit or implicit. In case f^-1 is not a function, g could be altered to return the solution with minimal r.

With possible adjustment of the constant C=c(X_0), iterate:

G <- f^-1(X)

X <- f(r^p,p*phi,p*theta) + C

Mandelbulb for f(G) as a sphere, MandelDonut for a torus. What else?

Trifold MandelKnot

The true Trifold Knot can be expressed probably only as a parametric function.

To get X from r,phi,theta, an initial guess is acquired via an empirical function.

Following Newton-Rhapson:

G_n+1 = G_n + J^-1 * (X-f(G_n)); where J is the Jacobian: partial df_i/dG_j

Previous step should work (sometimes) as initial guess. However something isn't working, I'll keep trying.

To get the idea:

Initial guess at parameters r,phi,theta

100 phi = atan2(x,y)

C

t2 = cos(1.5*phi)/3.8

t1 = pi/12.*cos(1.5*phi-t2) - phi/2. + pi/2.

t2 = -pi/12.*cos(1.5*phi+t2) - phi/2. + 1.5*pi

C

X1=3.*(/sin(t1)+2.*sin(2.*t1),cos(t1)-2.*cos(2.*t1),-sin(3.*t1)/)

X2=3.*(/sin(t2)+2.*sin(2.*t2),cos(t2)-2.*cos(2.*t2),-sin(3.*t2)/)

qnX1=(x-X1(1))**2+(y-X1(2))**2+(z-X1(3))**2

qnX2=(x-X2(1))**2+(y-X2(2))**2+(z-X2(3))**2

IF (qnX1.LT.qnX2) THEN

r = sqrt(qnX1)

phi = t1

theta=-atan2(z-X1(3),sqrt(x**2+y**2)-sqrt(X1(1)**2+X1(2)**2))

ELSE

r = sqrt(qnX2)

phi = t2

theta=-atan2(z-X2(3),sqrt(x**2+y**2)-sqrt(X2(1)**2+X2(2)**2))

END IF

j=1

Iterate to get parameters r,phi,theta

110 sphi(1)=sin(phi) !sin/cos lookup table

sphi(2)=sin(2.*phi)

sphi(3)=sin(3.*phi)

cphi(1)=cos(phi)

cphi(2)=cos(2.*phi)

cphi(3)=cos(3.*phi)

C

v0=1./sqrt(32.*cphi(3)+68.) !inverted norms

w0=1./sqrt(8.)/

/ sqrt(795.+36.*cos(9.*phi)+217.*cos(6.*phi)+652.*cphi(3))

C

v1=(-2.*sphi(1)+8.*sphi(2))*v0

v2=-(2.*cphi(1)+8.*cphi(2))*v0

C

w1=-12.*cphi(3)*(cphi(1)+4.*cphi(2))*w0

w2= 12.*cphi(3)*(sphi(1)-4.*sphi(2))*w0

w3= (-68.-32.*cphi(3))*w0

C

cth=cos(theta)

sth=sin(theta)

C New iteration

Xn = (/a*sphi(1)+2*a*sphi(2)+r*(v1*cth+w1*sth),

, a*cphi(1)-2*a*cphi(2)+r*(v2*cth+w2*sth),

, -a*sphi(3)+r*w3*sth/)

dX = (/Xn(1)-x,Xn(2)-y,Xn(3)-z/)

qndX = dX(1)**2 + dX(2)**2 + dX(3)**2

j=j+1

IF ((qndX.LT.1.D-24).OR.(j.LT.30)) GOTO 113

Calculating Jacobian + some tweaks to get it stable.

C dX/dr

Jac(1,1) = v1*cth+w1*sth

Jac(2,1) = v2*cth+w2*sth

Jac(3,1) = w3*sth

C dX/dtheta

Jac(1,3) = -r*(v1*sth+w1*cth)

Jac(2,3) = -r*(v2*sth+w2*cth)

Jac(3,3) = r*w3*cth

C

dv0dphi = -v0*48*sphi(3)

dv1dphi = (-2*cphi(1)+16*cphi(2))*v0 - v1*v0*dv0dphi

dv2dphi = ( 2*sphi(1)+16*sphi(2))*v0 - v2*v0*dv0dphi

C

dw0dphi = 4*w0*(-324*sin(9*phi)-1302*sin(6*phi)-1956*sin(3*phi))

dw1dphi = (36*sphi(3)*(cphi(1)+4*cphi(2))-12*cphi(3)*

* (-sphi(1)-8*sphi(2)))*w0-w1*dw0dphi*w0

dw2dphi = (-36*sphi(3)*(sphi(1)-4*sphi(2))+12*cphi(3)*

* (cphi(1)-8*cphi(2)))*w0-w2*dw0dphi*w0

dw3dphi = 96*sphi(3)*w0-w3*w0*dw0dphi

C dX/dphi

Jac(1,2) = a*cphi(1)+4*a*cphi(2)+r*cth*dv1dphi+r*sth*dw1dphi

Jac(2,2) = -a*sphi(1)+4*a*sphi(2)+r*cth*dv2dphi+r*sth*dw2dphi

Jac(3,2) = -3*a*cphi(3)+r*sth*dw3dphi

C

invJac(1,1) = Jac(2,2)*Jac(3,3)-Jac(3,2)*Jac(2,3)

invJac(1,2) = Jac(3,2)*Jac(1,3)-Jac(1,2)*Jac(3,3)

invJac(1,3) = Jac(1,2)*Jac(2,3)-Jac(2,2)*Jac(1,3)

detJac = 1./(Jac(1,1)*invJac(1,1)+Jac(2,1)*invJac(1,2)+

+ Jac(3,1)*invJac(1,3))

C Inverse Jacobian

invJac(1,1) = invJac(1,1)*detJac

invJac(1,2) = invJac(1,2)*detJac

invJac(1,3) = invJac(1,3)*detJac

invJac(2,1) = (Jac(2,3)*Jac(3,1)-Jac(2,1)*Jac(3,3))*detJac

invJac(2,2) = (Jac(1,1)*Jac(3,3)-Jac(1,3)*Jac(3,1))*detJac

invJac(2,3) = (Jac(1,3)*Jac(2,1)-Jac(1,1)*Jac(2,3))*detJac

invJac(3,1) = (Jac(2,1)*Jac(3,2)-Jac(2,2)*Jac(3,1))*detJac

invJac(3,2) = (Jac(1,2)*Jac(3,1)-Jac(1,1)*Jac(3,2))*detJac

invJac(3,3) = (Jac(1,1)*Jac(2,2)-Jac(1,2)*Jac(2,1))*detJac

dr = invJac(1,1)*dX(1)+invJac(1,2)*dX(2)+

+ invJac(1,3)*dX(3)

dphi = invJac(2,1)*dX(1)+invJac(2,2)*dX(2)+

+ invJac(2,3)*dX(3)

dtheta = invJac(3,1)*dX(1)+invJac(3,2)*dX(2)+

+ invJac(3,3)*dX(3)

IF (dtheta.LT.0.1) dr=0.5*dr*cos(dtheta)

IF (abs(dtheta).GT.1.5) dtheta=1.5*sign(1.,dtheta)

theta=mod(theta+dtheta,2*pi)

IF (dr+r.LT.0.) dr=0

r=r+dr

IF (r.GT.25) r=25

IF (abs(dphi).GT.0.1) dphi=0.1*sign(1.,dphi)

phi=phi+dphi

GOTO 110

And the very fractal

113 IF (ni.EQ.1) THEN

c=r

END IF

xc=r*(v1*cth+w1*sth)

yc=r*(v2*cth+w2*sth)

zc=r*w3*sth

r=r**3

IF ((r.GT.10.).OR.(ni.EQ.nimax)) GOTO 115

phi=phi*3.

theta=theta*3.

sphi(1)=sin(phi) !sin/cos lookup table

sphi(2)=sin(2.*phi)

sphi(3)=sin(3.*phi)

cphi(1)=cos(phi)

cphi(2)=cos(2.*phi)

cphi(3)=cos(3.*phi)

C

v0=1./sqrt(32.*cphi(3)+68.) !inverted norms

w0=1./sqrt(8.)/

/ sqrt(795.+36.*cos(9.*phi)+217.*cos(6.*phi)+652.*cphi(3))

C

v1=(-2.*sphi(1)+8.*sphi(2))*v0

v2=-(2.*cphi(1)+8.*cphi(2))*v0

C

w1=-12.*cphi(3)*(cphi(1)+4.*cphi(2))*w0

w2= 12.*cphi(3)*(sphi(1)-4.*sphi(2))*w0

w3= (-68.-32.*cphi(3))*w0

C

cth=cos(theta)

sth=sin(theta)

C New fractal iteration

x = a*sphi(1)+2*a*sphi(2)+r*(v1*cth+w1*sth) + xc

y = a*cphi(1)-2*a*cphi(2)+r*(v2*cth+w2*sth) + yc

z = -a*sphi(3)+r*w3*sth + zc

C

ni=ni+1

GOTO 100