import numpy as np
import matplotlib.pyplot as plt

u = np.linspace(0.001, 1, 180)
v = np.linspace(-np.pi, np.pi, 240)

U, V = np.meshgrid(u, v)

def gu(c, w):
    return np.exp(-((U - c) / w) ** 2)

def ab(c, w):
    d = np.angle(np.exp(1j * (V - c)))
    return np.exp(-(d / w) ** 2)

# Basic body/head radius
rho = (
    0.58
    + 0.34 * gu(0.30, 0.19)      # belly
    + 0.58 * gu(0.70, 0.18)      # head
    - 0.16 * gu(0.50, 0.10)      # narrower neck
    - 0.20 * gu(0.96, 0.08)      # top shaping
)

# Ear bumps
L = gu(0.94, 0.075) * ab(0, 0.23)
R = gu(0.94, 0.075) * ab(np.pi, 0.23)

# Arm bumps
LA = gu(0.48, 0.095) * ab(0, 0.30)
RA = gu(0.48, 0.095) * ab(np.pi, 0.30)

# Foot bumps
LF = gu(0.10, 0.07) * ab(0, 0.38)
RF = gu(0.10, 0.07) * ab(np.pi, 0.38)

# Parametrized surface
X = (
    rho * np.cos(V)
    + 0.55 * (L - R)
    + 0.38 * (LA - RA)
    + 0.32 * (LF - RF)
)

Y = (
    0.78 * rho * np.sin(V)
    - 0.20 * (LF + RF)
)

Z = (
    -1.65
    + 3.35 * U
    + 1.15 * (L + R)
    - 0.12 * (LA + RA)
    - 0.13 * (LF + RF)
)

# Plot
fig = plt.figure(figsize=(8, 9))
ax = fig.add_subplot(111, projection="3d")

ax.plot_surface(
    X, Y, Z,
    rstride=2,
    cstride=2,
    linewidth=0
)

ax.set_xlabel("x")
ax.set_ylabel("y")
ax.set_zlabel("z")
ax.set_title("Single Parametrized Full-Body Surface")

ax.view_init(elev=12, azim=-67)

plt.show()