R source code

library("threejs")

N     = 20000
theta = runif(N)*2*pi
phi   = runif(N)*2*pi
R     = 1.5
r     = 1.0

x = (R + r*cos(theta))*cos(phi)
y = (R + r*cos(theta))*sin(phi)
z = r*sin(theta)

d = 6
h = 6
t = 2*runif(N) - 1
w = t^2*sqrt(1-t^2)
x1 = d*cos(theta)*sin(phi)*w
y1 = d*sin(theta)*sin(phi)*w

i = order(phi)
j = order(t)
col = c(rainbow(length(phi))[order(i)],
        rainbow(length(t), start=0, end=2/6)[order(j)])

M = cbind(x=c(x, x1), y=c(y, y1), z=c(z, h*t))
scatterplot3js(M, size=0.1, color=col, bg="black", pch=".")