How to compute the normals to Costa's minimal surface? I am trying to draw Costa's minimal surface in high resolution using the PovRay raytracer. For this I need to compute points on the surface as well as the normals. It is relatively easy to compute the points with Mathematica, because it has Weierstrass functions built in:
g = WeierstrassInvariants[{1/2, I/2}];
e = WeierstrassP[1/2, g];

costa[x_, y_] :=
 With[
  {
   wz1 = WeierstrassZeta[x + I y, g],
   wz23 = 
    WeierstrassZeta[x + I y - 1/2, g] - 
     WeierstrassZeta[x + I y - I/2, g],
   wp = WeierstrassP[x + I y, g]
   },
  {
   (Pi x + Pi^2/(4 e) - Re[wz1] + Pi/(2 e)*Re[wz23])/2,
   (Pi y + Pi^2/(4 e) + Im[wz1] + Pi/(2 e)*Im[wz23])/2,
   (Sqrt[2 Pi]/4)*Log[Abs[(wp - e)/(wp + e)]]
   }
  ]

To draw Costa's surface you can do ParametricPlot3D[costa[x, y], {x, 0, 1}, {y, 0, 1}]:

Let me write costa[x,y] as $C(x,y)$. Presumably, to compute the normal at $C(x_0,y_0)$ we take the cross product $(\partial C/\partial x) \times (\partial C/\partial y)$ at $(x_0, y_0)$. This is not so easy to do with Mathematica because it does not know how to compute the derivatives of $\wp$ and $\zeta$.
How can I get explicit formulas for the normal? Do I really have to do the cross product of partial derivatives, or is there an easier way?
 A: In "A complete embedded minimal surface in R3 with genus one and three ends"
(J. Differential Geom. 21 (1985)) Hoffman and Meeks proved that the Costa Surface is embedded. They used (as Costa hisself did) the Weierstrass representation for this surface in terms of a function $g$ and a 1-form $fdz$ (see formula 4 in this paper). The function g is explicitly given (see theorem 1) as 
$g=\frac{2\sqrt{2\pi }\mathfrak P(0.5)}{\mathfrak P'}$
where $\mathfrak P$ and $\mathfrak P'$ are the Weierstrass $\mathfrak P$-function and its derivative of torus given by the square lattice $\mathbb Z +\mathbb Z i$. 
In general the function $g$ in the Weierstrass representation has the geometric meaning of being the stereographic projection (here you need to be careful to get the appropriate stereographic projection, i.e. the right formula) of the Gauss map of the corresponding minimal surface.
Of course mathematica knows $\mathfrak P'$ as well as $\mathfrak P.$
A: I went to the trouble of computing the normals with Mathematica (thanks to my colleague Pavle Saksida for teaching me how to compute derivatives). You can see the result at my GitHub costa-surface repository. The result is not pretty. Just for the fun of it, and to torture MathJax, here is the $x$-component of the normal at a point parametrized by $z \in \mathbb{C}$, as computed and simplified by Mathematica ($g$ is as computed above in the original question):
$$\frac{\sqrt{\pi}}{8\Re(\wp(1/2,g))\sqrt{2}}\cdot\left(\left(\pi\Im(\wp(z-i/2,g)-\wp(z-1/2,g))-2\Re(\wp(1/2,g))\Im(\wp(z,g))\right)\cdot \left(\Im(\wp'(z,g))\cdot\left(\frac{\Re(\wp(1/2,u)-\wp(z,g))}{\Im(\wp(1/2,u)-\wp(z,g))^2+\Re(\wp(1/2,g)-\wp(z,g))^2}+\frac{\Re(\wp(1/2,g)+\wp(z,g))}{\Im(\wp(1/2,u)+\wp(z,g))^2+\Re(\wp(1/2,g)+\wp(z,g))^2}\right)+\Re(\wp'(z,g))\cdot\left(\frac{\Im(-\wp(1/2,u)+\wp(z,g))}{\Im(\wp(1/2,u)-\wp(z,g))^2+\Re(\wp(1/2,g)-\wp(z,g))^2}-\frac{\Im(\wp(1/2,g)+\wp(z,g))}{\Im(\wp(1/2,u)+\wp(z,g))^2+\Re(\wp(1/2,g)+\wp(z,g))^2}\right)\right)-\left(2\Re(\wp(1/2,g))(\pi-\Re(\wp(z,g)))+\pi\Re(\wp(z-i/2,g)-\wp(z-1/2,g))\right) \cdot \left(\frac{\Im(-\wp(1/2,u)+\wp(z,g))\Im(\wp'(z,g))}{\Im(\wp(1/2,g)-
\wp(z,g))^2+\Re(\wp(1/2,g)-\wp(z,g))^2}-\frac{\Im(\wp(1/2,g)+\wp(z,g))\Im(\wp'(z,g))}{\Im(\wp(1/2,g)+\wp(z,g))^2+\Re(\wp(1/2,g)+\wp(z,g))^2}+\Re(\wp'(z,g))\cdot\left(\frac{\Re(-\wp(1/2,u)+\wp(z,g))}{\Im(\wp(1/2,u)-\wp(z,g))^2+\Re(\wp(1/2,g)-\wp(z,g))^2}-\frac{\Re(\wp(1/2,g)+\wp(z,g))}{\Im(\wp(1/2,u)+\wp(z,g))^2+\Re(\wp(1/2,g)+\wp(z,g))^2}\right)\right)\right)
$$
So I wonder how Weierstraß would have computed this.
