3D helix torus with hidden lines

Question

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.

Answer 1

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 }

Answer 2

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}

Answer 3

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;
  return toreturn;
}

helix helixAbout(helix center, int numloops, real radius) {
  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 ..);
  return toreturn;
}

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

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 helixradius = 0.5;
real safefactor = 1;
for (int i = 1; i < order; ++i)
  safefactor -= radiusfactor^i;
real saferadius = helixradius * safefactor;

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 helixradius = 0.5;
real safefactor = 1;
for (int i = 1; i < order; ++i)
  safefactor -= radiusfactor^i;
real saferadius = helixradius * safefactor;

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 helixradius = 0.5;
real safefactor = 1;
for (int i = 1; i < order; ++i)
  safefactor -= radiusfactor^i;
real saferadius = helixradius * safefactor;

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;
  return toreturn;
}

helix helixAbout(helix center, int numloops, real radius) {
  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 ..);
  return toreturn;
}

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

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 helixradius = 0.5;
real safefactor = 1;
for (int i = 1; i < order; ++i)
  safefactor -= radiusfactor^i;
real saferadius = helixradius * safefactor;

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).

Answer 4

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);

//inner radius
real R=2;
//outer radius
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}

Answer 5

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