Welcome to Fractal Forums

const h=0.5, k = 0.866025
a = -x; b = +h*x +k*y; c = +h*x - k*y;
// NOTE that a = -b-c !
ai = int(a); bi = int(b); ci = int(c);
a = abs(a-ai); b = abs(b-bi); c = abs(c-ci);
return (min(a,b,c)^2 + z*z - tuberad);

const h=0.5, k = 0.866025
a = -x; b = +h*x +k*y; c = +h*x - k*y;
ai = int(a); bi = int(b); ci = int(c);
a = abs(a-ai); b = abs(b-bi); c = abs(c-ci);
// also tried this (uncomment to see the effect) - varies pattern but it's pointless?
// ai = int(a+h); bi = int(b+h); ci = int(c+h);
if (ai mod 3 == 0) a = 1; //delete x tubes - else try a=0 to make a filled pattern
if (bi mod 3 == 0) b = 1; //delete skew tubes
if (ci mod 3 == 0) c = 1; //delete other tubes
// also tried this and did not work... (replace the previous 3 lines with this)
//if (ai mod 3 == 0) b = 1; 
//if (bi mod 3 == 0) c = 1;
//if (ci mod 3 == 0) a = 1;
return (min(a,b,c)^2 + z*z - tuberad);

const h=0.5, k = 0.866025
a = -x; b = +h*x +k*y; c = +h*x - k*y;
ai = int(a); bi = int(b); ci = int(c);
a = abs(a-ai); b = abs(b-bi); c = abs(c-ci);
if ((ai + bi mod 2) mod 3 == 0) a = 1;
if ((bi + ci mod 2) mod 3 == 0) b = 1;
if ((ci + ai mod 2) mod 3 == 0) c = 1;
return (min(a,b,c)^2 + z*z - tuberad);

static c[3]={-sqrt(3)/2,-0.5,sqrt(3)/4};
(x,y){

   reduce(&x,&y);
   fold(x,y,c);
   x=abs(x);
   y=abs(y);
   fold(x,y,c);
   d=fold(x,y,c);
      
   return d;
   //the 3D distance is: d3D=sqrt(z^2+d^2) . it is applicable whenever you have a 2D DE field that you want to see in 3D. Another trick is to do an intersection with a thin plate ;-)
}
reduce(&x,&y){//infinite tiling of the plane
   tile(x,1.5);
   tile(y,sqrt(3)/2);
}
fold(&x,&y,p[3]){
   t=x*p[0]+y*p[1]+p[2];
   d=t;
   t=2*min(t,0);
   x-=t*p[0];y-=t*p[1];
   return d;
}
tile(&x,a){//infinite fold
   x-=a;
   x=abs(x-(2*floor(x*0.5/a)+1)*a);
}

(x,y) // Infinite cubes illusion. Plot with Evaldraw!
m=0.5;f=sqrt(12);w=f/3;q=f/2; 
// find the tubes and store in a,b,c then find the int part of inclined axis
x = abs(x);
a = x;
b = ABS(x*.5+y*0.8660254037);
c = ABS(x*.5-y*0.8660254037);
ai = int(a+.5)-0;
bi = int(b+.5)-0;
ci = int(c+.5)-0;
a= abs(a-ai); b = abs(b-bi); c = abs(c-ci);
// a,b,c -> tubes now let's cut them to show only a part of them. This is a good way ...
u = abs(-y+q*ai);u=abs(u-int(u/f+m)*f)-w;a = max(a,u); // u is ortho to y-tubes then we cut using max() operator
u = abs(y*.5-x*0.8660254037+q*bi);u=abs(u-int(u/f+m)*f)-w;b = max(b,u);
u = abs(y*.5+x*0.8660254037+q*ci);u=abs(u-int(u/f+m)*f)-w;c = max(c,u);
// choose the min() so we stick together all parts
(sqrt((min(min(a,b),c))^2) - .01)^.2;

//DE for duals of uniform planar 2-3-6 triangle group tilings
//use with evaldraw : http://advsys.net/ken/download.htm
static c[3]={-sqrt(3)/2,-0.5,sqrt(3)/4};
(x,y){
   //fold transformation into 2-3-6 triangle group tiling. It works just like a Kaleidoscope.
      reduce(x,y);
      fold(x,y,c);
      x=abs(x);
      y=abs(y);
      fold(x,y,c);
   //end fold. It can be used as a 
   
   dx=y;//distance to x axis
   dy=x;//distance to y axis
   dc=x*c[0]+y*c[1]+c[2];//distance to third folding line
   
   d=dx;//hexagones
   //d=dy;//triangles
   //d=dc;//star like shapes
   
   //Any combinations of dx,dy and dc gives the 7 duals of uniform 2-3-6 tilings
   //For example:
   //d=min(dx,dy);
   //or any DE wich will be duplicated all around
   return d;
   
   //the 3D distance is:
   // d3D=sqrt(z^2+d^2) . circle cross section.
   // d3D=max(abs(z),d) . square.
   // d3D=abs(d)+abs(z) . diamond.
   //it is applicable whenever you have a 2D DE field that you want to see in 3D. Another trick is to do an intersection with a thin plate ;-)
}
reduce(&x,&y){//infinite tiling of the plane
   tile(x,1.5);
   tile(y,sqrt(3)/2);
}
fold(&x,&y,p[3]){
   t=x*p[0]+y*p[1]+p[2];
   d=t;
   t=2*min(t,0);
   x-=t*p[0];y-=t*p[1];
   return d;
}
tile(&x,a){//infinite fold
   x-=a;
   x=abs(x-(2*floor(x*0.5/a)+1)*a);
}

static a[2]={1,0};
static b[2]={0.5,0.5*sqrt(3)};
(x,y){
   inthex(&x,&y);
   abs(sin(3*x+y*10))//gen a color
}
inthex(&x,&y){//rouds to nearest hexagonal grid "integer" coordinates
   //convert to "sheared" coodinates (x,y) -> (l,m) by solving the system:
   //x=l*a[0]+m*b[0];
   //y=l*a[1]+m*b[1];

   idet=1/(a[0]*b[1]-b[0]*a[1]);
   m=(-x*a[1]+y*a[0])*idet;
   l=(x*b[1]-y*b[0])*idet;
   
   m0=floor(m);l0=floor(l);//coordinates of the lower left nearest "integer" corner
   
   //find nearest "integer" corner
   d2=(l0*a[0]+m0*b[0]-x)^2+(l0*a[1]+m0*b[1]-y)^2;
   nl=l0; nm=m0;nd=d2;
   d2=((l0+1)*a[0]+m0*b[0]-x)^2+((l0+1)*a[1]+m0*b[1]-y)^2;
   if(d2<nd) {nl=l0+1; nm=m0;nd=d2;}
   d2=(l0*a[0]+(m0+1)*b[0]-x)^2+(l0*a[1]+(m0+1)*b[1]-y)^2;
   if(d2<nd) {nl=l0; nm=m0+1;nd=d2;}
   d2=((l0+1)*a[0]+(m0+1)*b[0]-x)^2+((l0+1)*a[1]+(m0+1)*b[1]-y)^2;
   if(d2<nd) {nl=l0+1; nm=m0+1;nd=d2;}
   
   x=nl; y=nm;
}