Python script for Costa's minimal surface

# from https://mathstodon.xyz/@tpfto/115282798013947883

import mpmath as mpm
import numpy as np
import plotly.graph_objs as go

# prepare a polar mesh
rh = np.linspace(0.05, 0.95, 70)
th = np.linspace(-np.pi, np.pi, 90)
rG, tG = np.meshgrid(rh, th)
zz = rG * np.exp(1j * tG)

# Jacobi epsilon function, https://dlmf.nist.gov/22.16.E30
def ellipjep(zz, m):
    km = mpm.ellipk(m)
    em = mpm.ellipe(m)
    qm = mpm.qfrom(m=m)
    fac = mpm.mp.pi / (2 * km)
    uu = fac * zz
    zn = fac * mpm.jtheta(4, uu, qm, 1) / mpm.jtheta(4, uu, qm)
    res = zz * (em / km) + zn
    return res

# vectorize the elliptic functions
jaccn = np.vectorize((lambda u, m: complex(mpm.ellipfun("cn", u, m))))
jacdn = np.vectorize((lambda u, m: complex(mpm.ellipfun("dn", u, m))))
jacsn = np.vectorize((lambda u, m: complex(mpm.ellipfun("sn", u, m))))
jacep = np.vectorize((lambda u, m: complex(ellipjep(u, m))))

# period
cc = float(2 * mpm.ellipk(0.5))

# precompute elliptic functions over mesh
st = jacsn(cc * zz, 0.5)
ct = jaccn(cc * zz, 0.5)
dt = jacdn(cc * zz, 0.5)
et = jacep(cc * zz, 0.5)

# Costa minimal surface, https://doi.org/10.1007/BF02584707
x = np.real(
    (np.pi / 2 + (cc / 2) ** 2) * zz
    - cc / 2 * (et + ct * dt / st)
    - (np.pi * st * (1 / ct + ct)) / (4 * cc * dt)
)
y = np.imag(
    (np.pi / 2 - (cc / 2) ** 2) * zz
    + cc / 2 * (et + ct * dt / st)
    - (np.pi * st * (1 / ct + ct)) / (4 * cc * dt)
)
z = np.sqrt(np.pi / 8) * np.log(np.abs(ct**2))

# 'birdie' from https://tpfto.wordpress.com/2019/04/17/faking-a-birds-colors-and-miscellanea/
birdie = [
    "#3E26A8", "#4536D4", "#484CEF", "#4562FC", "#347AFC", "#2C90F0", "#21A3E3", "#0CB3D3", "#12BEB9",
    "#34C79C", "#5CCC76", "#94C94A", "#C8C129", "#F1BA37", "#FEC933", "#F5E327", "#F9FB14"
]

# and now, the plot...
surface = go.Surface(x=x, y=y, z=z, colorscale=birdie, showscale=False)
data = [surface]

psetxy = dict(
    gridcolor="#EEE8D5",
    zerolinecolor="#EEE8D5",
    showbackground=True,
    backgroundcolor="#839496",
    range=[-4, 4]
)
psetz = dict(
    gridcolor="#EEE8D5",
    zerolinecolor="#EEE8D5",
    showbackground=True,
    backgroundcolor="#839496",
    range=[-2.5, 2.5]
)
layout = go.Layout(
    title="Costa's minimal surface",
    scene=dict(
        aspectmode="manual",
        aspectratio=dict(x=1, y=1, z=0.625),
        xaxis=psetxy,
        yaxis=psetxy,
        zaxis=psetz
    ),
)

fig = go.Figure(data=data, layout=layout)
fig.show()