Better to do a mesh:
function parametricMesh(f, umin, umax, vmin, vmax, nu, nv)
u_ = LinRange(umin, umax, nu)
v_ = LinRange(vmin, vmax, nv)
vertices = Array{Float64}(undef, nu*nv, 3)
triangles = Array{Int64}(undef, 2*(nu-1)*(nv-1), 3)
k = 1
for i in 1:nv
v = v_[i]
for j in 1:nu
vertices[k,:] = f(u_[j], v)
k = k+1
end
end
k = 1
for i in 1:(nv-1)
for j in 1:(nu-1)
a = (i-1) * nu + j
b = (i-1) * nu + j + 1
c = i*nu + j + 1
d = i*nu + j
triangles[2*(k-1)+1,:] = [b, a, d]
triangles[2*k,:] = [c, b, d]
k = k+1
end
end
return (vertices = vertices, triangles = triangles)
end
function r(phi; a, b, m, n1, n2, n3)
return 1 / (abs(cos(m*phi/4)/a)^n2 + (abs(sin(m*phi/4)/b)^n3))^(1/n1)
end
p1 = (a=1, b=1, m=4, n1=10, n2=10, n3=10)
p2 = p1
function f(phi, theta)
x = r(theta; a=p1.a, b=p1.b, m=p1.m, n1=p1.n1, n2=p1.n2, n3=p1.n3) * cos(theta) *
r(phi; a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * cos(phi)
y = r(theta; a=p1.a, b=p1.b, m=p1.m, n1=p1.n1, n2=p1.n2, n3=p1.n3) * sin(theta) *
r(phi; a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * cos(phi)
z = r(phi; a=p2.a, b=p2.b, m=p2.m, n1=p2.n1, n2=p2.n2, n3=p2.n3) * sin(phi)
return [x, y, z]
end
ssmesh = parametricMesh(f, -pi/2, pi/2, -pi, pi, 50, 50)
using GLMakie
using GeometryBasics
vertices = map(Point3f, eachrow(ssmesh.vertices))
faces = map(TriangleFace, eachrow(ssmesh.triangles))
nrmls = normals(vertices, faces)
gb_mesh = GeometryBasics.Mesh(meta(vertices; normals = nrmls), faces)
mesh(gb_mesh, color = "yellow")