e1 = CBF(1/2).elliptic_p(CBF(I))
ii = CBF(I)

def fx(u, v):
    z = CBF(u) + ii*CBF(v)
    z0 = CBF(z).elliptic_zeta(ii)
    z1 = CBF(z - 1/2).elliptic_zeta(ii)
    zi = CBF(z - ii/2).elliptic_zeta(ii)
    return float(1/2*(-z0+pi*CBF(u)+pi**2/4/e1+pi/(2*e1)*(z1-zi)).real_part())

def fy(u, v):
    z = CBF(u) + ii*CBF(v)
    z0 = CBF(z).elliptic_zeta(ii)
    z1 = CBF(z-1/2).elliptic_zeta(ii)
    zi = CBF(z-ii/2).elliptic_zeta(ii)
    return float(1/2*(-ii*z0+pi*CBF(v)+pi**2/4/e1-pi/(2*e1)*(ii*z1-ii*zi)).real_part())

def fz(u, v):
    ell = (CBF(u)+ii*CBF(v)).elliptic_p(ii)
    return float(sqrt(CBF(2*pi))/4 * abs((ell-e1) / (ell+e1)).log())

P = parametric_plot3d((fx, fy, fz), (u, 0, 1), (v, 0, 1))