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