I am trying to draw [Costa's minimal surface](http://mathworld.wolfram.com/CostaMinimalSurface.html) 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}]`:

![enter image description here](https://i.stack.imgur.com/WGW51.png)

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?