Supertoroid.py

import bpy
import numpy as np

resU = 100
resV = 100
n1 = 1.2
n2 = 1.5
r0 = 2
r1 = 1

def C(theta, degree):
    return sgn(np.cos(theta))*np.power(np.abs(np.cos(theta)), degree)

def S(theta, degree):
    return sgn(np.sin(theta))*np.power(np.abs(np.sin(theta)), degree)

def sgn(x):
    if x < 0:
        return -1
    elif x == 0:
        return 0
    else:
        return 1

def getVertices(r0, r1, n1, n2, resU=100, resV=100):
    verts = []
    U = np.linspace(0, 2*np.pi, resU, False)
    V = np.linspace(0, 2*np.pi, resV, False)
    for u in U:
        for v in V:
            x = C(u, n1)*(r0+r1*C(v, n2))
            y = S(u, n1)*(r0+r1*C(v, n2))        
            z = r1*S(v, n2)
            verts.append([x, y, z])
    return verts

def getFaces(resU=100, resV=100):
    faces = []
    for i in range(resU):
        for j in range(resV):
            faces.append([i+j*resV, i+(j+1)%resV*resV, (i+1)%resU+(j+1)%resV*resV, (i+1)%resU+j*resV])
    return faces

msh = bpy.data.meshes.new(name="toroidmesh")
#msh = bpy.data.meshes.new("cubemesh")
msh.from_pydata(getVertices(r0, r1, n1, n2, resU, resV), [], getFaces(resU, resV))
msh.update()
obj = bpy.data.objects.new(name="Supertoroid", object_data=msh)
#obj = bpy.data.objects.new("cube", msh)
scene = bpy.context.scene
#scene.objects.link(obj)
scene.collection.objects.link(obj)