Costa's Minimal Surface

Weierstraß data

A minimal surface can be described often conveniently in terms of its Weierstraß data.
For Costa's surface, these consist of a meromorphic function G (the Gauß map) and a meromorphic height differential dh on the square torus. Given these data, a surface parametrization is obtained by the formula

F [z]=(1/2) Re(∫_0^z)(({{G - (1/G)}, {i (G + (1/G))}, {2}}))dh

The use of elliptic functions can be avoided if we allow for multivalued functions instead:

G[
          z_] := ρ (Sqrt[ (1/(z ((z^2) - 1)))]);
dh[z_] := (1/(-1 + (z^2)));

Here the Gauß map G is multivalued, its natural holomorphic continuation is defined on a square torus. For the moment, we restrict the domain of definition of G and dh to the first quadrant.

Without going into to much detail, we point out that there are four points where the Gauß map is singular, namely at 0,1,-1 and ∞. The point 0 corresponds to the central symmetry point of the surface where the coordinate axes meet. The points 1 and -1 correspond to the catenoidal ends, while ∞ corresponds to the middle planar end. The behaviour of a minimal surface at the ends is completely determined by the singularities of G and dh at these points.

We introduce the 1-forms

ϕ1[z_] :=  G[z] dh[z];
ϕ2[z_] :=  (1/G)[z] dh[z];

These can be integrated in terms of hypergeometric functions:

Φ1[w_] := 
        2 ρ I (Sqrt[w]) Hypergeometric2F1[(1/4), (3/2), (5/4), 
            (w^2)];
Φ2[
          w_] := -(2/(3  ρ)) I (w^(3/2)) Hypergeometric2F1[(1/2), (3/4), (7/4), (w^2)];

To avoid problems with branch cuts later, we will use instead:

Φ1[
        w_] := -((2 (Sqrt[w]) (Sqrt[1 - (w^2)]) Hypergeometric2F1[(1/4), (3/2), 
                (5/4), (w^2)])/(Sqrt[-1 + (w^2)]))
Φ2[
        w_] := -((2 (w^(3/2)) (Sqrt[-1 + (w^2)]) Hypergeometric2F1[
                (3/4), (1/2), (7/4), 
                (w^2)])/(3 (Sqrt[1 - (w^2)])))

To do the integration with Mathematica, it is advisable to rewrite first

ϕ1[z_] := ρ 
           (1/((-1 + (z^2)) (Sqrt[z (-1 + (z^2))])));
ϕ2[z_] := (1/ρ) 
        (z/(Sqrt[z (-1 + (z^2))]))

Then indeed

(∫_0^w) ϕ1[z] \[DifferentialD]z

(∫_0^w) ϕ2[z] \[DifferentialD]z

The Lopez-Ros Parameter

The so called Lopez-Ros parameter ρ is given by

ρ = Gamma[(3/4)]/((Sqrt[2]) Gamma[(5/4)]);

This choice of the parameter ensures that the image of a cycle encircling the points 0 and 1 encircling the points is a closed curve on the surface. To see this, we would like to integrate ϕ1 and ϕ2 from 0 to 1. However, this is not possible for ϕ1, as the singularity at 1 is not integrable. To circumvent this difficulty, we replace ϕ1 by a 1-form ϕ1b with the same periods:

ϕ1b[
        z_] := ρ (( 4 (z^(3/2)))/((Sqrt[(z^2) - 1]) (((1 + (z^2)))^2)))

This is true (up to some careful sign checking) because the difference of ϕ1b and ϕ1 is the deivative of a function which is meromorphic on the torus:

(ϕ1b[z] + ϕ1[z])/
          D[-((2 (Sqrt[z]) ρ)/((Sqrt[-1 + (z^2)]) (1 + 
                        (z^2)))), z] // PowerExpand // 
      Simplify

Then the definition of ρ is checked by

(∫_0^1) ϕ1b[
              z] \[DifferentialD]z == (∫_0^1) ϕ2[
              z] \[DifferentialD]z // FullSimplify

So we obtain finally the parametrization

ω1[
          z_]  := -(1/2) (Φ1[z] - Φ2[
                z]);
ω2[
          z_]  := (I/2) (Φ1[z] + Φ2[
                z] );
ω3[z_]  := 
        (1/2) Log[-1 + z] - (1/2) Log[1 + z];

F[z_]  :=   Re[{ω1[z], ω2[z], ω3[
            z]}]

Parametrizing the Surface

For drawing Costa's surface, we will exploit first its symmetries. First, the vertical coordinate planes are symmetry planes, and the 180°rotation about the y=x axes becomes also a reflective surface symmetry. These together decompose the surfaces into eigth congruent pieces, one of which is represented by the first quadrant.

Within this quadrant, three points are of special importance: The origin will be a corner of the surface patch, 1 will be mapped to the planar end, and ∞ to one of the catenoid ends. To get a parametrization of the surface patch adapted to this geometry, we need to foliate the first quadrant by coordinate lines which look like polar coordinates near the ends. This can be done using

p[w_]  := (Sqrt[1.  + E^( w)])

This map maps an infinite parallel strip to the first quadrant so that the ends of the strip are mapped to 1 and ∞.

<< Graphics`ComplexMap`

CartesianMap[p, {-2, 3}, {0, Pi}, 
        PlotRange -> All];

Finally, we need to set up the coordinate lines of the strip so that the image grid looks pretty. This is best done manually:

xlist = {-3.5, -3. , -2.5, -2. , -1.5, -1. , - .5, -0.3, -0.2, - .1, - .05, - .03, 
          0. ,  .03,  .05,  .1,  .2, 0.4, 0.6,  .8, 1. , 1.2, 1.5, 2. , 2.5, 
          3. , 3.5, 4. , 5. , 6. };

The negative values are mapped to the planar end, the positive to the catenoid end.

xlen = Length[xlist];
X[s_] := 
      If[s < 0.5  ∨ s > xlen + 0.4, -1, 
        xlist[[Round[s]]]]

This defines a step function, taking the wanted x-values near integral points. Similarly for the horizontal lines:

ylist = 
        Pi*{0, 0.08, 
            0.16,  .24,  .32,  .4,  .48,  .56,  .64,  .72,  .8,  .88,  .95,  .98,  .99,  .999,  .9999, 1};

ylen = Length[ylist];
Y[s_] := 
      If[s < 0.5  ∨ s > ylen + 0.4, -1, 
        ylist[[Round[s]]]]

Now we produce the fundamental piece:

gr1 = ParametricPlot3D[
            Evaluate[F[p[X[x] + I Y[y]]]],  {x, 1, xlen}, {y, 1, ylen},  
            PlotPoints -> {xlen,  ylen},  PlotRange -> All, 
            DisplayFunction -> Identity][[1]];

Arranging the Surface

Next we produce a single graphics object with  eight copies of the computed patch. W do this by applying reflections to the patch. As some reflections reverse the orientation, we correct this by reversing the order of the points in the respective polygons:

revpoly[Polygon[{v1_, v2_, v3_, v4_}]] := 
      Polygon[{v4, v3, v2, v1}]

First the reflection about the  y=x line:

gr2 = 
        Flatten[{gr1,  
            gr1 /.  {x_?NumericQ,  y_,  z_} :> {y, 
                  x, -z}}];

Then the four copies obtained by reflecting about tbe coordinate planes:

\(\(gr3 = 
        Flatten[{gr2, \ 
            Map[revpoly, 
              gr2 /. \ {x_?NumericQ, \ y_, \ z_} :> {\(-x\), y, 
                    z}], 
		\ gr2 /. \ {x_?NumericQ, \ y_, \ 
                  z_} :> {\(-x\), \(-y\), z}, 
		Map[revpoly, 
              gr2 /. \ {x_?NumericQ, \ y_, \ z_} :> {x, \(-y\), 
                    z}]}];\)\)

If you don't want to see the grid edges:

gr4 = Prepend[gr3, EdgeForm[]];

Show[Graphics3D[gr4], Boxed -> False, Axes -> False, 
        ViewPoint -> {1.5,  -3.0,  0.9}];

Finally, here is an GLExplorer version:

Needs["GLExplorer`GLRenderer`"]

GLShow[Graphics3D[gr3], LockLights -> True, PolygonEdges -> False, 
        Boxed -> False, Axes -> False, Tessellation -> True, 
        PlotRange -> All, ShadeModel -> Smooth, 
        DepthCompression -> False, 
        ImageSize -> {256, 256}];