# Draw a smooth surface

I would like to draw a smooth surface like https://mathematica.stackexchange.com/questions/20912/how-to-draw-a-higher-genus-surface but using Asymptote, TikZ, Metapost, PStricks, or other TeX friendly software. The following comes out pretty rough. If I increase nx it crashes (a frequent problem for me using Asymptote). settings.outformat = "pdf";
settings.render = 16;
size(8cm);
real R=1.2;
real r=.5;
real a=sqrt(2);
real Rsqr=R^2;
real rsqr=r^2;
real asqr=a^2;
real diffsqr=Rsqr-rsqr;
real sqrdiffsqr=diffsqr^2;
real twsumsqr=2*(rsqr+Rsqr);
real twdiffsqr=2*diffsqr;

import contour3;
currentprojection=perspective(4,6,8);
real f(real x, real y, real z)
{
real ysqr=y^2;
real zsqr=z^2;
real p=(-r - R + x)^2;
real q=(r + R + x)^2;
return
-asqr
+
(
sqrdiffsqr
- twsumsqr*(p + ysqr)
+ twdiffsqr*zsqr
+ (p+ysqr+zsqr)^2
)
*
(
sqrdiffsqr
- twsumsqr*(q + ysqr)
+ twdiffsqr*zsqr
+ (q+ysqr+zsqr)^2
);
}

material opaq=material(diffusepen=gray(0.5), emissivepen=gray(0.4), specularpen=black);

draw(
surface(
contour3(
f,
(-2*(r + R), -(r + R),-r - a),
(2*(r + R), (r + R), r + a),
nx=60)
),
surfacepen=opaq
);

• Could you update your question to show the parameretisation separately from the code? – cmhughes Mar 5 '15 at 16:33
• @cmhughes: The surface is not parametrized--it's plotted implicitly. – Charles Staats Mar 6 '15 at 1:59
• I am hoping to draw a smooth surface, with no corners, like a smooth neck between two tori. So I am not quite as satisfied as I hoped I would be with the union of two overlapping tori, because it has corners. – Benjamin McKay Mar 8 '15 at 14:20
• The original pdf is a little better than the png file you see above, but not as much as I would like. – Benjamin McKay Mar 8 '15 at 14:26

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 • Amazing. A giant improvement on my attempt. – Benjamin McKay Mar 11 '15 at 12:00
• I can see some tiny "stretch marks" at the closest points to the viewer on the "neck" in the first solution. Luckily they get smoothed out in the second solution. Very nice! – Benjamin McKay Mar 11 '15 at 13:17
• The "stretch marks" are due to the not well-balanced points. To give an idea it is similar when you take the square root on equidistribued point. – O.G. Mar 11 '15 at 13:19
• Out of curiosity (since I haven't tried it myself): Is the output in vector format or not? – Werner Mar 11 '15 at 14:43
• No, the output is bitmap here, you can choose the resolution. With asy you can have interactive OpenGL picture or with the -prc output an interactive PRC picture in the pdf document (only supported by AdobeReader, not by evince, okular...). – O.G. Mar 11 '15 at 15:12

Update The smoothcontour3 module has been incorporated into Asymptote version 2.33 (released 11 May 2015), so if you have that version or later, you should not need to download the module.

## Introducing smoothcontour3.asy

I have finally gotten this module for plotting smooth implicitly defined surfaces in a form that I am willing to share with others. It's not perfect, but I think the results are usually pretty good: settings.outformat = "png";;
settings.render = 16;
size(8cm);

import smoothcontour3;

// Erdos lemniscate of order n:
real erdos(pair z, int n) { return abs(z^n-1)^2 - 1; }

real h = 0.12;

real lemn3(real x, real y) { return erdos((x,y), 3); }

real f(real x, real y, real z) {
return lemn3(x,y)^2 + (16*abs((x,y))^4 + 1) * (z^2 - h^2);
}

draw(
implicitsurface(f,
a=(-3,-3,-3),
b=(3,3,3),
overlapedges=true
),
surfacepen=material(diffusepen=gray(0.5),
emissivepen=gray(0.4),
specularpen=gray(0.1)
)
);


Here's what happens when you apply it to the OP's original code (slightly modified; I extended the bounds by epsilon and added a tiny bit of reflection): settings.outformat="png";
settings.render=16;

import smoothcontour3;
size(8cm);

real R=1.2;
real r=.5;
real a=sqrt(2);
real Rsqr=R^2;
real rsqr=r^2;
real asqr=a^2;
real diffsqr=Rsqr-rsqr;
real sqrdiffsqr=diffsqr^2;
real twsumsqr=2*(rsqr+Rsqr);
real twdiffsqr=2*diffsqr;

currentprojection=perspective(4,6,8);
real f(real x, real y, real z)
{
real ysqr=y^2;
real zsqr=z^2;
real p=(-r - R + x)^2;
real q=(r + R + x)^2;
return
-asqr
+
(
sqrdiffsqr
- twsumsqr*(p + ysqr)
+ twdiffsqr*zsqr
+ (p+ysqr+zsqr)^2
)
*
(
sqrdiffsqr
- twsumsqr*(q + ysqr)
+ twdiffsqr*zsqr
+ (q+ysqr+zsqr)^2
);
}

material softgray = material(diffusepen=gray(0.5),
emissivepen=gray(0.4),
specularpen=gray(0.1));

triple aa = (-2*(r + R), -(r + R),-r - a) - 1e-2*(1,1,1);
triple bb = (2*(r + R), (r + R), r + a) + 1e-2*(1,1,1);

draw(implicitsurface(f, aa, bb, overlapedges=true),
surfacepen=softgray);


There is a vector graphics version. The necessity of using meshpen is somewhat unfortunate, but translucency becomes feasible and it gives you some idea how the algorithm works. (Note that the algorithm is adaptive, unlike any of Asymptote's built-in plotting algorithms.) settings.outformat="pdf";
settings.render=0;

import smoothcontour3;
size(8cm);

real R=1.2;
real r=.5;
real a=sqrt(2);
real Rsqr=R^2;
real rsqr=r^2;
real asqr=a^2;
real diffsqr=Rsqr-rsqr;
real sqrdiffsqr=diffsqr^2;
real twsumsqr=2*(rsqr+Rsqr);
real twdiffsqr=2*diffsqr;

currentprojection=perspective(4,6,8);
real f(real x, real y, real z)
{
real ysqr=y^2;
real zsqr=z^2;
real p=(-r - R + x)^2;
real q=(r + R + x)^2;
return
-asqr
+
(
sqrdiffsqr
- twsumsqr*(p + ysqr)
+ twdiffsqr*zsqr
+ (p+ysqr+zsqr)^2
)
*
(
sqrdiffsqr
- twsumsqr*(q + ysqr)
+ twdiffsqr*zsqr
+ (q+ysqr+zsqr)^2
);
}

material softgray = material(diffusepen=gray(0.5) + opacity(0.9),
emissivepen=gray(0.4),
specularpen=gray(0.1));

triple aa = (-2*(r + R), -(r + R),-r - a) - 1e-2*(1,1,1);
triple bb = (2*(r + R), (r + R), r + a) + 1e-2*(1,1,1);

draw(implicitsurface(f, aa, bb, n=20, nx=30),
surfacepen=softgray,
meshpen=blue+0.2pt);


Here's a description of the implicitsurface function, since the code currently has no external documentation (except this answer, I suppose):

The implicit surface has three required parameters: a function f (of type either real(triple) or real(real, real, real)) and two triples a and b. It returns a plot of the zero locus of the function within the rectangular solid with opposite corners at a and b.

[Note: if f is overloaded with both signatures, the version with signature real(triple) will be used.]

Optional parameters:

• int n - the number of initial segments in each of the x, y, z directions.
• bool overlapedges - if true, the patches of the surface are slightly enlarged to compensate for an artifact in which the viewer can see through the boundary between patches. (Some of this may actually be a result of edges not lining up perfectly, but I'm fairly sure a lot of it arises purely as a rendering artifact.)
• int nx - override n in the x direction
• int ny - override n in the y direction
• int nz - override n in the z direction
• int maxdepth - the maximum depth to which the algorithm will subdivide in an effort to find patches that closely approximate the true surface. (The algorithm is adaptive).

The algorithm is designed to deal with a differentiable function with a smooth zero locus. It can often (but not always) cope with singularities in the zero locus, at the expense of slowing down the computation dramatically. On one occasion I have seen it cope with a continuous, piecewise differentiable function, although it took a while to compute. Whatever you do, do not pass it a piecewise constant function—it will take forever and return nothing good.

Original answer: Implicitly defined surfaces always come out rough using Asymptote. (I've actually started a module to do it better, but there's no guarantee I'll ever get it to work.) So far as I know, none of the other TeX-friendly software you list even attempts to plot implicitly defined surfaces. So your best bet for TeX-friendly software is probably the SageTeX package.

• Amazing. Thanks. I hope you find some time for mathematics. Computer programming is a dangerous temptation for a mathematician. – Benjamin McKay May 2 '15 at 9:33
• A little too dangerous in my case--I'm a software engineer now. :) – Charles Staats May 2 '15 at 15:50
• Wonderful work ! – O.G. May 3 '15 at 19:32

Another PSTricks solution:

\documentclass[pstricks]{standalone}
\usepackage{pst-solides3d}
\begin{document}

\begin{pspicture}(-8,-4)(8,4)
\psset{viewpoint=50 10 40 rtp2xyz,Decran=50,lightsrc=viewpoint,solidmemory,
object=tore,r1=3,r0=1,ngrid=36 60,fillcolor=[rgb]{1 0.5 0}}
\psSolid[name=toreA,action=none](0,3.5,0)
\psSolid[name=toreB,action=none](0,-3.5,0)
\psSolid[object=fusion,base=toreA toreB,grid=false]
\end{pspicture}

\end{document}


Use latex->dvips->ps2pdf or use xelatex but that will take very long time PSTricks can not smooth the surfaces. One can decrease the size of the grid [ngrid] to approach it, but the computation may become too long.

However, PSTrink allows to merge the objects (for instance two tores - similar or not). In the first example, the computation is done by Ghostscript. To reduce the computation time, as in the second example, one can save the files ( uncomment the line % the first time); the saving takes place during the viewing with Ghostscript.  \documentclass{article}
\usepackage{pst-solides3d}
\begin{document}
\begin{center}
\begin{pspicture}(-5,-4)(5,4)
\psset[pst-solides3d]{viewpoint=50 0 60 rtp2xyz,Decran=50,lightsrc=50 20 60 rtp2xyz}
\psset{solidmemory}
%\pstVerb{/n1 36 def /n2 60 def}%
%\psSolid[object=tore,
%         r1=3,r0=1,
%         ngrid=n1 n2,
%         action=none,
%         plansepare={[0 1 0 0]},
%         name=coupeplaneA](0,3.5,0)
%         action=writesolid,file=toreA,
%         rm=0]
%\psSolid[object=tore,
%         r1=3,r0=1,
%         ngrid=n1 n2,
%         action=none,
%         plansepare={[0 1 0 0]},
%         name=coupeplaneB](0,-3.5,0)
%         action=writesolid,file=toreB,
%         rm=0]
\psSolid[object=datfile,file=toreA,action=none,name=toreA1,fillcolor=yellow!50]%
\psSolid[object=datfile,file=toreB,action=none,name=toreB1,fillcolor=yellow!50]%
\psSolid[object=fusion,base=toreA1 toreB1,grid=true,linecolor={[rgb]{0 0.5 0}}]
\end{pspicture}
\end{center}
\begin{center}
\begin{pspicture}(-5,-4)(5,4)
\psset[pst-solides3d]{viewpoint=20 20 60 rtp2xyz,Decran=20}%
\psset{solidmemory}
\psSolid[object=datfile,file=toreA,action=none,name=toreA1,fillcolor={[rgb]{1 0.5 0}}]%
\psSolid[object=datfile,file=toreB,action=none,name=toreB1,fillcolor={[rgb]{1 0.5 0}}]%
\psSolid[object=fusion,base=toreA1 toreB1,grid,
intersectiontype=0,
intersectionplan={ -1 0.05 1 {/iP exch def [0 0 1 iP] } for},
intersectioncolor=(black)]
\end{pspicture}
\end{center}
\end{document}


I'm not sure what you mean by TeX friendly software. You can use Sage through LaTeX using the sagetex package but, unfortunately, I get output similar to your Asymptote code. If I work with Sage outside of LaTeX, like you're doing with Asymptote, then the result is better. Go to a Sage Cell Server and copy/paste the following code in:

var('u,v');
a,b = 1.2,.5 # outer and inner radii
x = (a + b*cos(u))*cos(v)
y = (a + b*cos(u))*sin(v)
z = b*sin(u)

var('t')
T = parametric_plot3d([x,y,z], (u,0,2*pi),(v,0,2*pi), opacity=1,aspect_ratio=[1,1,1],color="yellow",plot_points=[100,100],frame=False)
M = T + T.translate(2,1.5,0)
M.show(frame=False)


and then press "Evaluate" to get this: The double torus here is made by taking a plot of the torus and attaching it to a translation of the torus, so you can fiddle with it to your heart's content. Below the image, to the left, there's a link to download the image as a PNG file (which can be used when compiling with pdflatex). The image is clearer but you don't have the same control over the picture as you do if Sage is installed on your computer. Using Sage directly lets you right click on the image (see below) and do all sorts of things with it: rotate, zoom in, view it as a 3D anaglyph and so on. By rotating and zooming you can adjust the image to your satisfaction. This image doesn't have the problems that you encounter with Asymptote. If you've got red-cyan glasses for 3-D you can get the 3-D effect, too. You can play with the numbers in the aspect_ratio to make the double torus flatter or fatter. • If you want a smooth parametric plot (or union thereof), Asymptote can do that at least as well as Sage. But that surface is not smooth--you can see the corners. – Charles Staats Mar 6 '15 at 7:27
• Your posted answer was, "implicitly defined surfaces always come out rough" which seems to indicate the user had to live with the really rough picture above. The Sage code seemed like a noticeable improvement in smoothness, if not smooth as you note. The other answers posted seem to have some roughness as well. If Asymptote can do "at least as well" as the Sage code, why not provide the code that could help reduce the roughness problem that Benjamin McKay notes? – DJP Mar 6 '15 at 16:14
• First, in my mind, the key aspect of the question was that it was asking for an implicitly plotted surface. I've thought a lot in the past about how to plot a smooth genus-2 surface, and so far come up blank. Also, because of my training as a theoretical mathematician, I consider these pictures a good deal less smooth than the OP's: in the OP the roughness is an artifact, whereas a union of two tori is intrinsically non-smooth. – Charles Staats Mar 7 '15 at 1:25

It's fairly straight-forward to draw a torus in Sketchup. The advantage is that you can orbit through 3D space a find the best viewpoint, before printing the output to PDF (a quick pdfcrop will trim out any boundary whitespace). Replicating it is just as easy.

Double torus in Sketchup: PDF view: Triple torus in Sketchup: PDF view: 