Interesting question. Putting two torus is not a big deal as soon as you have some 3D routines and a lot of 3D solutions can do it. The remark of Charles Staats is relevant.
I tried to make a sufficiently smooth transition between two torus.
To do this I take an (x,y,z) parametrization of the left part which is possible for x belonging to [-R-a,-R], compute the interval in which y lies and then obtain the z value. The modification consists to add a x dependant function null in -R and having a null derivative in -R-a.
From a mathematical point of view the result is C^1 (I think)
Second Solution. I remember the interpolation in Tchebychev points versus equidistant points. Because for x fixed the curve is closed in some parts of an ellipse, I change the y parametrization so that the points are more well-balanced (it was the reason I put sqrt(y) in the first attempt). So the solution is to take sin(pi/2*y) and then it is possible to use the smooth-spline functionality of Asymptote. Please find the (always dirty and to improve) code
size(200);
import graph3;
currentprojection=perspective(5,4,4);
real R=3;
real a=1;
triple f(pair t) {
triple z;
z= ((R+a*cos(t.y))*cos(t.x),(R+a*cos(t.y))*sin(t.x),a*sin(t.y));
return z;
}
surface s=surface(f,(0,0),(2pi,2pi),Spline);
surface ss=shift((2*R+2*a,0,0))*surface(f,(0,0),(2pi,2pi),Spline);
draw(s,green,render(compression=Low,merge=true));
draw(ss,green,render(compression=Low,merge=true));
triple f11 (pair z)
{
real y,ty,x,yy;
triple zz;
if ((z.y>=0) & (z.y<=1))
{
x=z.x+1/2*(z.x+R)^2/a;
y=sin(pi*(z.y)/2)*sqrt((R+a)^2-x*x);
zz=(-z.x,y,sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}
if ((z.y<0) & (z.y>=-1))
{
x=z.x+1/2*(z.x+R)^2/a;
y=-sin(pi*abs(z.y)/2)*sqrt((R+a)^2-x*x);
zz=(-z.x,y,sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}
if ((z.y>1) && (z.y<=2))
{
x=z.x+1/2*(z.x+R)^2/a;
yy=2-z.y;
y=sin(pi*(yy)/2)*sqrt((R+a)^2-x*x);
zz=(-z.x,y,-sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}
if ((z.y>1) && (z.y<=2))
{
zz=(-z.x,y,-sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}
if ((z.y>=-2) && (z.y<-1))
{
x=z.x+1/2*(z.x+R)^2/a;
yy=-2-z.y;
y=-sin(pi*(abs(yy))/2)*sqrt((R+a)^2-x*x);
zz=(-z.x,y,-sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}
return zz;
}
triple f21(pair z)
{
triple z=f11(z);
return (-z.x,z.y,z.z);
}
int N=20;
surface st=surface(f11,(-R-a,-2),(-R,2),N,2*N,Spline);
surface sstb=shift((2*R+2*a,0,0))*surface(f21,(-R-a,-2),(-R,2),N,2*N,Spline);
surface pont1=surface(st,sstb);
draw(pont1,green, render(compression=Low,merge=true));
You can observe that the two torus are drawn and a 'bridge' is added.
The code is faster and the result looks like
First solution. However I cannot use the surface(...,Spline) routine to have a true smooth surface. I think it is due in fact to some artefacts of the Bézier patch parametrization (De Boor reference). Perhaps that using some monotonic option to Spline is the solution but I have to admit that I do not have the time.
Im y opinion, with the standard surface(....,100,100), the result is satisfactory.
I am sure that it is possible to improve it.
Please find the (dirty and to clean) code
size(200);
import graph3;
currentprojection=perspective(5,4,4);
real R=3;
real a=1;
triple f(pair t) {
triple z;
z= ((R+a*cos(t.y))*cos(t.x),(R+a*cos(t.y))*sin(t.x),a*sin(t.y));
return z;
}
bool active(pair pos) {return (R+a*cos(pos.y))*cos(pos.x)<=3.2;}
bool active2(pair pos) {return (R+a*cos(pos.y))*cos(pos.x)>=-3.2;}
surface s=surface(f,(0,0),(2pi,2pi),200,200, active);
surface ss=shift((2*R+2*a,0,0))*surface(f,(0,0),(2pi,2pi),200,200,active2);
draw(s,green,render(compression=Low,merge=true));
draw(ss,green,render(compression=Low,merge=true));
triple f11 (pair z)
{
real y,ty,x;
if (z.y>=0)
{
x=z.x+1/2*(z.x+R)^2/a;
y=sqrt(z.y)*sqrt((R+a)^2-x*x);
}
else
{
x=z.x+1/2*(z.x+R)^2/a;
y=-sqrt(abs(z.y))*sqrt((R+a)^2-x*x);
}
return (-z.x,y,sqrt(a^2-(R-sqrt((x)^2+y^2))^2));
}
triple f12(pair z)
{
triple z=f11(z);
return (z.x,z.y,-z.z);
}
triple f21(pair z)
{
triple z=f11(z);
return (-z.x,z.y,z.z);
}
triple f22(pair z)
{
triple z=f11(z);
return (-z.x,z.y,-z.z);
}
int N=100;
surface st=surface(f11,(-R-a,-1),(-R,1),N,N);
surface sst=surface(f12,(-R-a,-1),(-R,1),N,N);
surface sstb=shift((2*R+2*a,0,0))*surface(f21,(-R-a,-1),(-R,1),N,N);
surface sstbm=shift((2*R+2*a,0,0))*surface(f22,(-R-a,-1),(-R,1),N,N);
surface pont1=surface(st,sst,sstb,sstbm);
draw(pont1,green);//blue);
and the result
