# 3D helix torus with hidden lines

How can I generate the above image using TikZ (possibly other methods will work if they produce better looking results)?
I would like to have a nicely shaded (not so dark image) without gridlines and possibly enough opacity on the torus to see the hidden part of the helix.

The image given is very close to what I want though and I would just try to make it visually a little better.

-
Just in case it might be useful, here's the Maple worksheet that I used to generate that image. – Scott H. Sep 7 '12 at 7:31
@ScottH. Thanks – Jubao Sep 8 '12 at 18:15

Here's a Sketch/TikZ approach.

Running sketch on this file:

def helix {
def n_segs 600
sweep [draw=orange] { n_segs, rotate(24*360 / n_segs, (1.5,0,0), [0,0,1]), rotate(1*360/n_segs, (0,0,0), [0,1,0]) } (2.01,0,0)
}

def torus {
def n_segs 60
sweep [draw=none, fill=cyan, fill opacity=0.75] {n_segs, rotate(360/n_segs, (0,0,0), [0,1,0])}
sweep {n_segs, rotate(360/n_segs, (1.5,0,0), [0,0,1])}
(2,0,0)
}

put { view((10,4,2)) } {{helix} {torus}}

global { language tikz }


generates a .tex file which can be compiled using pdflatex.

The helix winding around the helix winding around a torus

can be generated using

def helix {
def n_segs 10000
sweep [draw=orange] {
n_segs,
rotate(1000*360 / n_segs, (2,0,0), [0,1,0]),
rotate(24*360 / n_segs, (1.5,0,0), [0,0,1]),
rotate(1*360/n_segs, (0,0,0), [0,1,0])
} (2.04,0,0)
}

def torus {
def n_segs 50
sweep [draw=none, fill=cyan, fill opacity=0.75] {n_segs, rotate(360/n_segs, (0,0,0), [0,1,0])}
sweep {n_segs, rotate(360/n_segs, (1.5,0,0), [0,0,1])}
(1.9,0,0)
}

put { view((10,4,2)) } {{torus} {helix}}

global { language tikz }

-
How big is the TikZ file? :-) – Joseph Wright Sep 10 '12 at 15:56
@JosephWright: Hehe, I tried posting it, but the system wouldn't let me. It's about 250kB, so not much fun to edit by hand. The resulting PDF is 47kB. – Jake Sep 10 '12 at 16:04
I'm impressed: draws faster than the PSTrick one! – Joseph Wright Sep 10 '12 at 16:13
OK. Thanks for your effort. I am still waiting for the pure TikZ. If it is impossible then I will assign you the bounty. – kiss my armpit Sep 10 '12 at 16:13
Wow. The helix winding around the other winding around the torus is incredible. This is most likely impossible in a 'normal' graphics program such as Illustrator or Inkscape and quite hard even in Blender or other 3D programs. – Alexander Sep 10 '12 at 16:49

The example shows the function

 x(u,v)=(R1 + (R0 +RL*sin(u))*sin(k*v))*cos(v)-RL*cos(u)*sin(v)
y(u,v)=(R1 + (R0 +RL*sin(u))*sin(k*v))*sin(v)+RL*cos(u)*cos(v)
z(u,v)=(R0 + RL*sin(u))*cos(k*v)


with the parameter setting shown in the example. RL: radius of the coil line; R1: Torus outer; R0: Torus inner radius; k:number of coils

run it with xelatex or latex>dvips>ps2pdf (takes some time to run!)

\documentclass{minimal}
\usepackage{pst-solides3d}
\pagestyle{empty}
\begin{document}

\begin{pspicture}[solidmemory](-6,-4)(6,4)
\psset{viewpoint=30 0 15 rtp2xyz,Decran=30,lightsrc=viewpoint}
\psSolid[object=tore,r1=5,r0=1,ngrid=36 36,
fillcolor=blue!30,action=none,name=Torus]%
%\axesIIID(4.5,4.5,0)(5,5,4)
\codejps{/R1 5 def /RL 0.05 def /R0 1.1 def /k 25 def}%
\defFunction[algebraic]{helix}(u,v)
{(R1 + (R0 +RL*sin(u))*sin(k*v))*cos(v)-RL*cos(u)*sin(v)}
{(R1 + (R0 +RL*sin(u))*sin(k*v))*sin(v)+RL*cos(u)*cos(v)}
{(R0 + RL*sin(u))*cos(k*v)}
\psSolid[object=surfaceparametree,
base=0 6.2831853 0 6.2831853,
linecolor=blue,linewidth=0.01,fillcolor=yellow,
ngrid=0.8 0.01,function=helix,action=none,name=Helix]%
\psSolid[object=fusion,base=Torus Helix,grid=false]
%\gridIIID[Zmin=-3,Zmax=3,showAxes=false](-2,2)(-2,2)
\end{pspicture}

\begin{pspicture}[solidmemory](-6,-6)(6,6)
\psset{viewpoint=30 0 90 rtp2xyz,Decran=30,lightsrc=viewpoint}
\psSolid[object=tore,r1=5,r0=1,ngrid=36 36,
fillcolor=blue!30,action=none,name=Torus]%
%\axesIIID(4.5,4.5,0)(5,5,4)
\codejps{/R1 5 def /RL 0.05 def /R0 1.1 def /k 25 def}%
\defFunction[algebraic]{helix}(u,v)
{(R1 + (R0 +RL*sin(u))*sin(k*v))*cos(v)-RL*cos(u)*sin(v)}
{(R1 + (R0 +RL*sin(u))*sin(k*v))*sin(v)+RL*cos(u)*cos(v)}
{(R0 + RL*sin(u))*cos(k*v)}
\psSolid[object=surfaceparametree,
base=0 6.2831853 0 6.2831853,
linecolor=blue,linewidth=0.01,fillcolor=yellow,
ngrid=0.8 0.01,function=helix,action=none,name=Helix]%
\psSolid[object=fusion,base=Torus Helix,grid=false]
%\gridIIID[Zmin=-3,Zmax=3,showAxes=false](-2,2)(-2,2)
\end{pspicture}

\end{document}


an animation is here: http://tug.org/PSTricks/main.cgi?file=Animation/gif/gif

With \psSolid[object=fusion,base=Torus Helix,grid=false,opacity=0.5]  (setting transparency) and a thinner helix (decrease /RL) one gets

and just for fun with

\listfiles
\documentclass{minimal}
\usepackage{pst-solides3d}
\begin{document}
\begin{pspicture}[solidmemory](-6.5,-3.5)(6.5,3)
\psset{viewpoint=30 0 15 rtp2xyz,Decran=30,lightsrc=viewpoint}
\psSolid[object=tore,r1=5,r0=1,ngrid=36 36,tablez=0 0.05 1 {} for,
zcolor= 1 .5 .5 .5 .5 1,action=none,name=Torus]
\pstVerb{/R1 5 def /R0 1.2 def /k 20 def /RL 0.15 def /kRL 40 def}%
\defFunction[algebraic]{helix}(t)
{(R1+R0*cos(k*t))*sin(t)+RL*sin(kRL*k*t)}
{(R1+R0*cos(k*t))*cos(t)+RL*cos(kRL*k*t)}
{R0*sin(k*t)+RL*sin(kRL*k*t)}
\psSolid[object=courbe,
resolution=7800,
fillcolor=black,incolor=black,
r=0,
range=0 6.2831853,
function=helix,action=none,name=Helix]%
\psSolid[object=fusion,base=Torus Helix,grid]
\end{pspicture}
\end{document}


-
Is there any way to fix the artifacts where the yellow helix seems to be missing in the first picture? I would also like to use a parameteric curve instead of a solid(I tried to modify the code but always get errors). I'd like the yellow helix to look more like a grid line on the torus rather than a separate object. – Jubao Sep 8 '12 at 18:18
\pstVerb{/R1 5 def /R0 1.2 def /k 20 def /RL 0.1 def /kRL 40 def}% – Herbert Sep 11 '12 at 14:55
then you have not the current pst-solides3d.tex|pro files. Compare with the ones from texnik.dante.de – Herbert Sep 19 '12 at 13:52
I am also having trouble obtaining the image above. I update the pst-solides3d.tex / pro files from the website you mention. Should I also update the rubans files ? – Mathusalem Jun 3 '13 at 11:12
see my edited answer with the current file list, which can be read at the end of the logfile. – Herbert Jun 3 '13 at 11:21

Here's a work-in-progress. It's missing the most important thing: the visibility of the helix. At first I thought that the points where it becomes (in)visible are evenly spaced, but they are not quite. I think I'll have to do some vector algebra to (hopefully) find a solution. The torus is not a 3D object, it's made from many almost transparent circular rings.

## Code

\documentclass[tikz,border=5mm]{standalone}

\begin{document}

\newcommand{\xangle}{-30}
\newcommand{\yangle}{210}
\newcommand{\zangle}{90}

\newcommand{\xlength}{1}
\newcommand{\ylength}{1}
\newcommand{\zlength}{1}

\pgfmathsetmacro{\xx}{\xlength*cos(\xangle)}
\pgfmathsetmacro{\xy}{\xlength*sin(\xangle)}
\pgfmathsetmacro{\yx}{\ylength*cos(\yangle)}
\pgfmathsetmacro{\yy}{\ylength*sin(\yangle)}
\pgfmathsetmacro{\zx}{\zlength*cos(\zangle)}
\pgfmathsetmacro{\zy}{\zlength*sin(\zangle)}

\begin{tikzpicture}
[   x={(\xx cm,\xy cm)},
y={(\yx cm,\yy cm)},
z={(\zx cm,\zy cm)},
]

\pgfmathsetmacro{\RO}{1.1}
\pgfmathsetmacro{\RI}{5}
\pgfmathsetmacro{\K}{10}

\pgfmathsetmacro{\doubleRO}{2*\RO}

\foreach \h in {0,0.01,...,\doubleRO}
{   \pgfmathsetmacro{\pm}{sqrt(\h*(\doubleRO-\h))}
\fill[opacity=0.007,blue,even odd rule] (0,0,\h-\RO) circle (\RI+\pm) (0,0,\h-\RO) circle (\RI-\pm);
}

\xdef\initialx{\RI}
\xdef\initialy{0}
\xdef\initialz{\RO}

\foreach \v in {0.1,0.2,...,360.1}
{   \pgfmathsetmacro{\newx}{(\RI + \RO*sin(\K*\v))*cos(\v)}
\pgfmathsetmacro{\newy}{(\RI + \RO*sin(\K*\v))*sin(\v)}
\pgfmathsetmacro{\newz}{\RO*cos(\K*\v)}
\pgfmathsetmacro{\mycolor}{cos(\v)*50+50}
\draw[red!\mycolor!green,thick] (\initialx,\initialy,\initialz) -- (\newx,\newy,\newz);
\xdef\initialx{\newx}
\xdef\initialy{\newy}
\xdef\initialz{\newz}

}

\end{tikzpicture}

\end{document}


## (insufficient) Result

-
if you want hidden lines and surfaces then you have to build small line or polygon segments of all objects, then build the direction vector and sort all calculated vectors. And in the end you have to fill the canvas with the line and polygon segments depending to the sorted list of direction vectors. – Herbert Sep 11 '12 at 13:17
see my deleted answer for the same limitations using pgfplots – cmhughes Sep 11 '12 at 18:58
The torus is nice and you got the windings. Maybe put a solid black torus behind the blue one so it looks shaded rather than glowing – Jubao Sep 16 '12 at 19:34

Updated: Includes workarounds for previous difficulties with compiling complex paths.

Here's an Asymptote approach that allows nth order helixes. I show examples of a first-order helix (which wraps once around the torus, as in the original question), a second-order helix, and a third-order helix.

Here's the code (configured for a first-order helix):

settings.outformat = "png";
settings.render = 16;
settings.prc = false;
real unit = 2cm;
unitsize(unit);

import graph3;

void drawsafe(path3 longpath, pen p, int maxlength = 400) {
int length = length(longpath);
if (length <= maxlength) draw(longpath, p);
else {
int divider = floor(length/2);
drawsafe(subpath(longpath, 0, divider), p=p, maxlength=maxlength);
drawsafe(subpath(longpath, divider, length), p=p, maxlength=maxlength);
}
}

struct helix {
path3 center;
path3 helix;
int numloops;
int pointsperloop = 12;
/* t should range from 0 to 1*/
triple centerpoint(real t) {
return point(center, t*length(center));
}
triple helixpoint(real t) {
return point(helix, t*length(helix));
}
triple helixdirection(real t) {
return dir(helix, t*length(helix));
}
/* the vector from the center point to the point on the helix */
triple displacement(real t) {
return helixpoint(t) - centerpoint(t);
}
bool iscyclic() {
return cyclic(helix);
}
}

path3 operator cast(helix h) {
return h.helix;
}

helix helixcircle(triple c = O, real r = 1, triple normal = Z) {
helix toreturn;
toreturn.center = c;
toreturn.helix = Circle(c=O, r=r, normal=normal, n=toreturn.pointsperloop);
toreturn.numloops = 1;
}

helix toreturn;
toreturn.numloops = numloops;
from toreturn unravel pointsperloop;
toreturn.center = center.helix;
int n = numloops * pointsperloop;
triple[] newhelix;
for (int i = 0; i <= n; ++i) {
real theta = (i % pointsperloop) * 2pi / pointsperloop;
real t = i / n;
triple ihat = unit(center.displacement(t));
triple khat = center.helixdirection(t);
triple jhat = cross(khat, ihat);
triple newpoint = center.helixpoint(t) + radius*(cos(theta)*ihat + sin(theta)*jhat);
newhelix.push(newpoint);
}
toreturn.helix = graph(newhelix, operator ..);
}

int loopfactor = 20;
helix wrap(helix input, int order, int initialloops = 10, real initialradius = 0.6, int loopfactor=loopfactor) {
helix toreturn = input;
int loops = initialloops;
for (int i = 1; i <= order; ++i) {
loops *= loopfactor;
}
}

currentprojection = perspective(12,0,6);

helix circle = helixcircle(r=2, c=O, normal=Z);

/* The variable part of the code starts here. */
int order = 1;    // This line varies.
real safefactor = 1;
for (int i = 1; i < order; ++i)

helix todraw = wrap(circle, order=order, initialradius = helixradius);    // This line varies (optional loopfactor parameter).

surface torus = surface(Circle(c=2X, r=0.99*saferadius, normal=-Y, n=32), c=O, axis=Z, n=32);
material toruspen = material(diffusepen=gray, ambientpen=white);
draw(torus, toruspen);

drawsafe(todraw, p=0.5purple+linewidth(1pt));  // This line varies (linewidth only).


The output:

For a second-order helix, change the last portion of the above code to the following:

/* The variable part of the code starts here. */
int order = 2;    // This line varies.
real safefactor = 1;
for (int i = 1; i < order; ++i)

helix todraw = wrap(circle, order=order, initialradius = helixradius, loopfactor=40);    // This line varies (optional loopfactor parameter).

surface torus = surface(Circle(c=2X, r=0.99*saferadius, normal=-Y, n=32), c=O, axis=Z, n=32);
material toruspen = material(diffusepen=gray, ambientpen=white);
draw(torus, toruspen);

drawsafe(todraw, p=0.5purple+linewidth(0.6pt));  // This line varies (linewidth only).


The output:

For a third-order helix, change the "variable" portion of the code to the following:

/* The variable part of the code starts here. */
int order = 3;    // This line varies.
real safefactor = 1;
for (int i = 1; i < order; ++i)

helix todraw = wrap(circle, order=order, initialradius = helixradius);    // This line varies (optional loopfactor parameter).

surface torus = surface(Circle(c=2X, r=0.99*saferadius, normal=-Y, n=32), c=O, axis=Z, n=32);
material toruspen = material(diffusepen=gray, ambientpen=white);
draw(torus, toruspen);

drawsafe(todraw, p=0.5purple+linewidth(0.2pt));  // This line varies (linewidth only).


The output:

One final comment: the rendering mechanism in Asymptote for translucent surfaces does not do a good job with more than one layer of translucent surface (i.e., when a translucent surface is obscuring another translucent surface (or piece of one)). Here's an example of the difficulty:

This is a bit less noticeable in mrc's answer, but it's still there, at least for now.

However, this difficulty can be mitigated (when only one surface is involved and it has no self-intersections) by sorting the constituent patches so that they are drawn in order of distance from the camera:

settings.outformat = "png";
settings.render = 16;
settings.prc = false;
real unit = 2cm;
unitsize(unit);

import graph3;

void drawsafe(path3 longpath, pen p, int maxlength = 400) {
int length = length(longpath);
if (length <= maxlength) draw(longpath, p);
else {
int divider = floor(length/2);
drawsafe(subpath(longpath, 0, divider), p=p, maxlength=maxlength);
drawsafe(subpath(longpath, divider, length), p=p, maxlength=maxlength);
}
}

void sort(surface s) {
projection P = currentprojection;
//The following code is copied from three_surface.asy
// Sort patches by mean distance from camera
triple camera=P.camera;
if(P.infinity) {
triple m=min(s);
triple M=max(s);
camera=P.target+camerafactor*(abs(M-m)+abs(m-P.target))*unit(P.vector());
}

real[][] depth=new real[s.s.length][];
for(int i=0; i < depth.length; ++i)
depth[i]=new real[] {abs(camera-s.s[i].cornermean()),i};

depth=sort(depth);
//end of copied code

int[] permutation = sequence(new int(int i) {return (int)depth[i][4];}, depth.length);

int[][] inversionTool = new int[permutation.length][5];
for (int i = 0; i < permutation.length; ++i)
inversionTool[i] = new int[] {permutation[i], i};
inversionTool = sort(inversionTool);
int inverse(int i) {return inversionTool[i][6];};

patch[] sortedS = new patch[depth.length];
for (int i = 0; i < sortedS.length; ++i) {
sortedS[i] = s.s[permutation[i]];
}
s.s = sortedS;

for (int[] currentrow : s.index)
for (int i = 0; i < currentrow.length; ++i)
currentrow[i] = inverse(currentrow[i]);
}

struct helix {
path3 center;
path3 helix;
int numloops;
int pointsperloop = 12;
/* t should range from 0 to 1*/
triple centerpoint(real t) {
return point(center, t*length(center));
}
triple helixpoint(real t) {
return point(helix, t*length(helix));
}
triple helixdirection(real t) {
return dir(helix, t*length(helix));
}
/* the vector from the center point to the point on the helix */
triple displacement(real t) {
return helixpoint(t) - centerpoint(t);
}
bool iscyclic() {
return cyclic(helix);
}
}

path3 operator cast(helix h) {
return h.helix;
}

helix helixcircle(triple c = O, real r = 1, triple normal = Z) {
helix toreturn;
toreturn.center = c;
toreturn.helix = Circle(c=O, r=r, normal=normal, n=toreturn.pointsperloop);
toreturn.numloops = 1;
}

helix toreturn;
toreturn.numloops = numloops;
from toreturn unravel pointsperloop;
toreturn.center = center.helix;
int n = numloops * pointsperloop;
triple[] newhelix;
for (int i = 0; i <= n; ++i) {
real theta = (i % pointsperloop) * 2pi / pointsperloop;
real t = i / n;
triple ihat = unit(center.displacement(t));
triple khat = center.helixdirection(t);
triple jhat = cross(khat, ihat);
triple newpoint = center.helixpoint(t) + radius*(cos(theta)*ihat + sin(theta)*jhat);
newhelix.push(newpoint);
}
toreturn.helix = graph(newhelix, operator ..);
}

int loopfactor = 20;
helix wrap(helix input, int order, int initialloops = 10, real initialradius = 0.6, int loopfactor=loopfactor) {
helix toreturn = input;
int loops = initialloops;
for (int i = 1; i <= order; ++i) {
loops *= loopfactor;
}
}

currentprojection = perspective(12,0,6);

helix circle = helixcircle(r=2, c=O, normal=Z);

/* The variable part of the code starts here. */
int order = 1;    // This line varies.
real safefactor = 1;
for (int i = 1; i < order; ++i)

helix todraw = wrap(circle, order=order, initialradius = helixradius);    // This line varies (optional loopfactor parameter).

surface torus = surface(Circle(c=2X, r=0.99*saferadius, normal=-Y, n=32), c=O, axis=Z, n=32);
sort(torus);
material toruspen = material(diffusepen=gray + opacity(0.5), ambientpen=white);
draw(torus, toruspen);

drawsafe(todraw, p=0.4magenta+linewidth(1pt));  // This line varies (linewidth only).

-
Thanks for mentioning that. I was wondering about that---it's incredibly noticeable for opacity(0.9), but I chose opacity(0.3) to minimize the weirdness. – mrc Dec 14 '13 at 1:33
Also: I'm trying to get better at asymptote. Would you mind explaining the idea behind your approach a little bit? I was going to just come up with the parametric equation for the 2nd order curve (following math.stackexchange.com/questions/143897/wrapping-curves) but that got way too complicated . . . – mrc Dec 14 '13 at 1:40
@mrc: As discussed in the question you link to, the key to wrapping a curve by a helix is to have a moving frame. In my approach, I assume that the curve to be wrapped (say, **v**(t)) was itself constructed as a helix (the "center" **c**(t)). I also assume that when both curves are reparametrized so that t ranges from 0 to 1, the vector from **c**(t) to **v**(t) is orthogonal to **v'**(t). Thus, this vector provides the **i**(t) vector (called ihat in the code), **v'**(t) provides **k**(t), and **j**(t) is obtained as a cross product. Once we have a continuous moving frame,... – Charles Staats Dec 14 '13 at 3:46
... the new helix is parametrized as h(t) = i cos(t) + j sin(t). The only reason I use path3s for the intermediate helixes rather than functions of type triple(real) is so that I can use the function dir(path3, real), which returns the unit tangent vector to the path at the specified t-value. – Charles Staats Dec 14 '13 at 3:52
@mrc: You might be interested in my newest addition to this answer. – Charles Staats Feb 15 '14 at 1:08

Here's an approach using asymptote. It's pretty basic, but if you open it in Adobe Reader you can rotate the torus!

\documentclass[12pt]{article}
\usepackage{asymptote}

\begin{document}

\begin{center}
\begin{asy}[width=0.5\textwidth]
import graph3;

size(200,0);
currentprojection=orthographic(4,0,2);

real R=2;
real a=0.75;

//surface:
triple f(pair t) {
return ((R+a*cos(t.y))*cos(t.x),(R+a*cos(t.y))*sin(t.x),a*sin(t.y));
}

//path:
real x(real t) {return cos(t/15)*(R + a*cos(t));}
real y(real t) {return sin(t/15)*(R + a*cos(t));}
real z(real t) {return a*sin(t);}

pen p=blue+opacity(0.33);
// make surface and path
surface s=surface(f,(0,0),(2pi,2pi),8,8,Spline);
path3 q=graph(x,y,z,0,30*pi,operator ..);

// draw surface and path
draw(s,p);
draw(q);

\end{asy}
\end{center}
\end{document}


-