Las ecuaciones paramétricas de la superficie de Costa son:

Costa[z_] := Module[{p0, p1, p2, q1, q2, q3},
    p0 = WeierstrassP[z, {189.0729999999999, 0}];
    p1 = WeierstrassPPrime[z, {189.0729999999999, 0}];
    p2 = WeierstrassZeta[Conjugate[z], {189.0729999999999, 0}];
    {q1, q2, q3} = Append[cr[(6.28319*p0)/p1 + 3.1415*z - p2],
    1.253309999999999*Log[Abs[(p0 - 6.875189999999999)/(p0 + 6.875189999999999)]]];
    {{q1, q2, q3}, {-q1, q2, q3}, {-q1, -q2, q3}, {q1, -q2, q3}}]

donde z es un número complejo y

cr[z_] := {Re[z], Im[z]}