How to actually use the Weierstrass-Enneper parameterization to draw a minimal surface?

Question

I'm interested in drawing (with Mathematica for example) the generalized Scherk saddle tower with threefold symmetry, a shape that I find very attractive. In an article (see here) I found the following Weierstrass-Enneper parameterization:

\begin{eqnarray} x&=&\Re \int_{z_0}^{z_1}\frac{1}{2}f(z)(1-g^2(z))\;\mathrm{d}z\\ y&=&\Re \int_{z_0}^{z_1}\frac{\mathrm{i}}{2}f(z)(1+g^2(z))\;\mathrm{d}z\\ z&=&\Re \int_{z_0}^{z_1}f(z)g(z)\;\mathrm{d}z,\\ \end{eqnarray} with $$ f(z)=\frac{1}{z^6+1}\quad \text{and} \quad g(z)=z^2. $$

Unfortunately, I have no idea how to use these formulas to actually draw the surface. To begin with, which integration contour should I use? What do $z_0$ and $z_1$ stand for? And once I know the values of these integrals (either analytically or numerically), how do I obtain the $(x(u,v),y(u,v),z(u,v))$ parameterization that I can feed to Mathematica?

Any help would be greatly appreciated!

Answer 1

Here's an approach which, while naive, can be improved considerably. I'll simply define the functions $f(w)$ and $g(w)$ and then the functions $x(w)$, $y(w)$, and $z(w)$ in terms of the integrals, making sure the values at zero are zero. I changed the complex variable $z$ to $w$ to avoid conflict with the Cartesian coordinate $z$. The input parameters will be the polar coordinates $r$ and $t$ with domain specified to describe the unit disk.

Clear[f, g, x1, y1, z1];
f[w_] = 1/(w^6 + 1);
g[w_] = w^2;
x1[w_] = Integrate[f[w] (1 - g[w]^2)/2, w];
(* Make value zero at zero *)
x1[w_] = x1[w] - x1[0];
y1[w_] = Integrate[I*f[w] (1 + g[w]^2)/2, w];
z1[w_] = Integrate[f[w] g[w], w];
ParametricPlot3D[
 Re[{x1[r*Exp[I*t]], y1[r*Exp[I*t]], z1[r*Exp[I*t]]}],
 {r, 0, 1}, {t, 0, 2 Pi}]

Well, not so bad really. We clearly have some issues with branch cuts and we get only part of the surface, which is really all to be expected. The easiest thing to do is restrict our input domain to a region that avoids the branch cuts and then use the resulting piece, together with the known symmetry of the desired result to assemble the final figure. Here's a piece that, I think, should be a reasonable building block.

basePic = ParametricPlot3D[
  Re[{x1[r*Exp[I*t]], y1[r*Exp[I*t]], z1[r*Exp[I*t]]}],
  {r, 0, 1}, {t, Pi/3, Pi}, PlotRange -> Automatic,
  PlotPoints -> 100, MaxRecursion -> 8]

We can place three copies of this, symmetrically rotated about the $z$-axis, reflect the resulting figure about a plane normal to the $z$-axis, and finally translate that a few times to obtain a nice image. Here is the result.

And here is the code, which I deferred since it is rather, well, codey. It's really all just a matter of directly manipulating the graphics primitives that form the basePic.

primitives = First[Normal[basePic]];
rotate[prim_, t_] := prim /. {
    Line[pts_] :> Line[RotationMatrix[t, {0, 0, 1}].# & /@ pts],
    Polygon[pts_] :> Polygon[RotationMatrix[t, {0, 0, 1}].# & /@ pts]
    };
primitives = 
  Table[Rotate[primitives, t, {0, 0, 1}], {t, {0, 2 Pi/3, 4 Pi/3}}];
primitives = {primitives, primitives /. {
     Line[pts_] :> 
      Line[ReflectionMatrix[{0, 0, 1}].(# - {0, 0, Re[z1[Exp[2 I]]]}) +
          {0, 0, Re[z1[Exp[2 I]]]} & /@ pts],
     Polygon[pts_, VertexNormals -> vn_] :> 
      Polygon[ReflectionMatrix[{0, 0, 
             1}].(# - {0, 0, Re[z1[Exp[2 I]]]}) +
          {0, 0, Re[z1[Exp[2 I]]]} & /@ pts,
       VertexNormals -> (-{#[[1]], #[[2]], #[[3]]} & /@ vn)]}
   };
primitives = {primitives, primitives /. {
     Line[pts_] :> Line[# - {0, 0, 4 Re[z1[Exp[2 I]]]} & /@ pts],
     Polygon[pts_, VertexNormals -> vn_] :> 
      Polygon[# - {0, 0, 4 Re[z1[Exp[2 I]]]} & /@ pts, 
       VertexNormals -> vn]},
   primitives /. {
     Line[pts_] :> Line[# + {0, 0, 4 Re[z1[Exp[2 I]]]} & /@ pts],
     Polygon[pts_, VertexNormals -> vn_] :> 
      Polygon[# + {0, 0, 4 Re[z1[Exp[2 I]]]} & /@ pts, 
       VertexNormals -> vn]}
   };
Graphics3D[primitives, 
   PlotRange -> {{-0.7, 0.7}, {-0.7, 0.7}, {-1.2, 1.8}},
   Axes -> True, ViewPoint -> {0.95, -3.2, 0.54}]

Answer 2

$\newcommand{\Cpx}{\mathbf{C}}$Not exactly an answer, but a collection of miscellanea too extensive for a comment.

Caveat: I've never done any of this myself; YMMV. :)

You can fix $z_{0}$ arbitrarily: a different choice amounts to an additive constant (i.e., a translation of the surface). The complex quantity $z_{1} = u + iv$ comprises the surface parameters. For brevity, let's denote the parametrization by $\Phi$, so that (loosely) $$ \Phi(z_{1}) = \Phi(u, v) = \bigl(x(u, v), y(u, v), z(u, v)\bigr) = (x, y, z). $$ The proximate difficulty for plotting is that the component functions aren't closed formulas (e.g., explicit elementary functions), but definite integrals. If you're able to calculate the integrals as elementary functions (here your integrands are rational, so in principle you're in business, though the partial fractions may be a mess...), then you simply instruct your plotting software accordingly.

If you're less lucky, you'd presumably ask the software to find numerical values for the integrals (as functions of the "upper limit" $z_{1} = u + iv$).

Now we come to the issue of which integration contour to use. For starters, let's write $\Phi_{\gamma}(z_{1})$ to denote the value obtained by integrating over a contour $\gamma$. To understand conceptually the effect a different choice of $\gamma$ may have, let's pull out some "heavy machinery", covering spaces.

The respective integrands in the parametrization are defined in some open subset $U$ of $\Cpx$ (here the complement of the sixth roots of $-1$). The surface defined by $\Phi$ is a covering space of $U$. The nature of the covering is determined by the real parts of the residues of the integrands around each pole.

For example, suppose you pick $z_{0} = 0$ (which naively seems like a good choice), and you want to find $\Phi(z)$ for some $z$. Let $\gamma_{0}$ and $\gamma_{1}$ be arbitrary contours from $0$ to $z$ in $U$. We need to understand how $\Phi_{\gamma_{0}}(z)$ differs from $\Phi_{\gamma_{1}}(z)$. The key idea is, the difference $$ \Phi_{\gamma_{0}}(z) - \Phi_{\gamma_{1}}(z) $$ is an integral around a loop $\gamma$ based at $0$ (follow $\gamma_{0}$ forward, then $\gamma_{1}$ backward), and as such is the sum of (real parts of) residues weighted by the winding numbers of $\gamma$ about each pole. (In the specific example at hand, the components of $\Phi$ contain $\log$ terms, which give rise to the periods.)

As a practical matter, it's probably easiest to calculate the real part of the residue of $\Phi$ about each pole. These vectors are periods; the surface parametrized by $\Phi$ is invariant under translation by its periods. (Addition of a period is precisely what happens when a contour winds around a pole.)

Now pick a dense, simply-connected subset of $U$ (e.g., remove six radial rays from each pole to infinity), plot the portion of the surface obtained by using radial straight-line contours, and translate the resulting surface by the periods to get as large a piece as desired.

Here are some links of possible interest. (I have no affiliation with any of them.)

3D-XploreMath, and the gallery of minimal surfaces.

Scherk saddle tower at Shapeways (a commercial site offering 3D printed models for sale).