# www.contextfreeart.org/gallery/view.php?id=3731

startshape Clifford_Torus CF::Impure=1 CF::Background=[b 0.95] CF::MinimumSize=0.1 CF::ColorDepth=16 CF::Time=[time 0 1] CF::Frame=0 Frame_Rot=frame()*360 // Frame rotation angle

N = 4*12 // numbers of holes * 12 Betw = 0.1 // Space between hexaeders Deg = 295+Frame_Rot // Rotate this angle Axis = Nrm3((-2,9,-1)) // about this axis Nsun = Nrm3(( 2,5, 3)) // Direction Light (2-sided) CamZ = 8 // (0,0,CamZ) Camera position CF::Size=[[s 1.4 s 4 3]]

// 3 Colorings: i=0, i=1, i=2 Hf(i)=select(i, 115, 85, 55) // Hue Sf(i)=select(i, 0.6, 0.5, 0.6) // Saturnation Bf(i)=select(i, 0.8, 0.9, 0.8) // Brightness Trans=0.33 // Transparency

//—————————————————- MinD = 0.1*CamZ // Minimal distance to Camera Pcam =(0,0,CamZ) // Camera position

//==================================================== // Clifford Torus (Area toward infinity is removed) //==================================================== shape Clifford_Torus { loop i=2*N [] loop j=N []

if (!(mod(i,3)&&mod(j,3))) // Cut outs & No infinit
 if( SawT(i-j,2*N)+SawT(i+j,2*N) <= 1.5*N)
   Hexaeder(i,j)[]

}

shape Hexaeder(i,j) { // Size depending Space between

H =0.5*(1-Betw)    m=-H    p=H  
im=i-H  ip=i+H    jm=j-H  jp=j+H  
Pc=Rotate(Torus(i ,j ,0))  // Centre point
P0=Rotate(Torus(im,jm,m))  P1=Rotate(Torus(ip,jm,m))
P2=Rotate(Torus(im,jp,m))  P3=Rotate(Torus(ip,jp,m))
P4=Rotate(Torus(im,jm,p))  P5=Rotate(Torus(ip,jm,p))
P6=Rotate(Torus(im,jp,p))  P7=Rotate(Torus(ip,jp,p))
if (CamZ-Pc[2]>MinD) { // Prevent too near to Camera
  FaceColor(0,Pc,(P0,P2,P3,P1),i,j)[] // back
  FaceColor(1,Pc,(P4,P0,P1,P5),i,j)[] // bottom
  FaceColor(2,Pc,(P3,P7,P5,P1),i,j)[] // right
  FaceColor(3,Pc,(P3,P2,P6,P7),i,j)[] // top
  FaceColor(4,Pc,(P4,P6,P2,P0),i,j)[] // left
  FaceColor(5,Pc,(P7,P6,P4,P5),i,j)[] // front

} }

shape FaceColor(F, vector3 Pcen, vector12 Quad, i,j){

c = select(F, 0,1,1,1,1,2)        // Face coloring
Ax= mod(2*i,N)==0||mod(2*j,N)==0  // boolean Axis
B = if(Ax, Bf(c)/2.2, Bf(c))      // Darken the Axis
A = -Trans*!Ax
FaceLight(Pcen,Quad)[h Hf(c) sat Sf(c) b B a A]

}

shape FaceLight(vector3 Pcen, vector12 Quad) {

P0=Quad[0,3] P1=Quad[3,3] P2=Quad[6,3] P3=Quad[9,3]
Pmid=Scl3(0.25,Add3(Add3(P0,P1),Add3(P2,P3)))
Dmid=Sub3(Pmid,Pcen)  // Direction Centre->Midpoint
Nmid=Nrm3(Dmid)       // Face Normal approximation

Dcam=Sub3(Pcam,Pmid)  // Direction Midpoint->Camera
Ncam=Nrm3(Dcam)       // If transparant or Camera-
if (Trans || Dot3(Nmid,Ncam)>0) { //  directed face 
  Diff=abs(Dot3(Nmid,Nsun))          // Diffuse 
  Nrfl=Sub3(Scl3(2*Diff,Nmid),Nsun)  // Reflection
  Spec=max(0,Dot3(Ncam,Nrfl))        // Specular 
  SAT =-0.25*Diff^1.5 - 0.50*Spec^1.5
  BR  = 1.00*Diff^1.5 + 0.25*Spec^6 - 0.6 

  Dis =Len3(Sub3(Pmid,Pcam))  // Distance to Pcam
  Scl =Len3(Dmid)/(Dis/CamZ)  // Scale stroke width
  Quadrangle(Quad,Scl)[sat SAT b BR z -Dis]

} }

path Quadrangle(vector12 Quad, Scl) {

MOVETO(Proj(Quad[0,3]))  LINETO(Proj(Quad[3,3]))
LINETO(Proj(Quad[6,3]))  LINETO(Proj(Quad[9,3]))
CLOSEPOLY() FILL[]
STROKE(0.012*Scl, CF::RoundCap)[b -0.7]

} //==================================================== // Parametric model of Clifford Torus in 3D //—————————————————- vector3 Torus(i,j,k)=let(

u=i*180/N; v=j*180/N; w=k*180/N;
Cu=cos(u); Cv=cos(v); Cw=cos(w);
Su=sin(u); Sv=sin(v); Sw=sin(w);
p0=Cu*Cv ; p1=Cu*Sv ; p2=Su*Cv ; p3=Su*Sv ;
n0=-p3   ; n1=p2    ; n2=p1    ; n3=-p0   ;
r1=Cw*p1+Sw*n1; r2=Cw*p2+Sw*n2; r3=Cw*p3+Sw*n3;
r0=Cw*p0+Sw*n0;       // (1,0,0,0) point at infinity
d =1-min(r0,0.99);    // Prevent infinity  
// A pleasant X,Y,Z orientation of the torus
X=(r2+r1)/sqrt(2); Y=r3; Z=(r2-r1)/sqrt(2);
(X/d, Y/d, Z/d)       // Stereographic projection

)

vector2 Proj(vector3 V)=let( // Canvas Projection

dZ= max(0.05*CamZ,CamZ-V[2]); // force in front
(V[0]*CamZ/dZ, V[1]*CamZ/dZ)  //          of Camera

)

vector3 Rotate(vector3 V)=let( // Rotate Deg Degrees

W=cos(Deg/2);  Sin=sin(Deg/2);   // about Axis 
X=Sin*Axis[0]; Y=Sin*Axis[1];  Z=Sin*Axis[2];
V0=V[0]; V1=V[1]; V2=V[2];
( V0 -2*((Y*Y+Z*Z)*V0 -(X*Y-Z*W)*V1 -(Z*X+Y*W)*V2)
, V1 -2*((Z*Z+X*X)*V1 -(Y*Z-X*W)*V2 -(X*Y+Z*W)*V0)
, V2 -2*((X*X+Y*Y)*V2 -(Z*X-Y*W)*V0 -(Y*Z+X*W)*V1)

) )

//==================================================== // Some functions //—————————————————- number SawT(X,N)= // Triangular Saw

abs(X-N*floor(X/N)-N/2)

number Len3(vector3 A)= // Length

sqrt(A[0]^2+A[1]^2+A[2]^2)

vector3 Nrm3(vector3 A)=let( // Normalize

L=Len3(A); ( A[0]/L, A[1]/L, A[2]/L) )

vector3 Scl3(S,vector3 A)= // Scale

( S*A[0], S*A[1], S*A[2] )

vector3 Add3(vector3 A, vector3 B)= // Add vectors

( A[0]+B[0] , A[1]+B[1] , A[2]+B[2] )

vector3 Sub3(vector3 A, vector3 B)= // Subtract vector

( A[0]-B[0] , A[1]-B[1] , A[2]-B[2] )

number Dot3(vector3 A, vector3 B)= // Dot product

( A[0]*B[0] + A[1]*B[1] + A[2]*B[2] )