# Drawing a torus with Julia

Posted on December 15, 2022 by Stéphane Laurent

I already showed how to draw a torus whose equator passes through three given points with three.js, Asymptote, R, POV-Ray, and Haskell. Here I have a method for Python, using PyVista, but it is a bit useless because once you know how to draw a circular path passing through three points you can “tubify” it with PyVista, that is to say you can add a radius to the path. Similarly, you can do that with the R package rgl by using the cylinder3d function.

Now I want to show how to do that with Julia. I’m using the Meshes package because I find it easy to use.

Here is the code:

using LinearAlgebra
using Meshes

# plane passing by points p1, p2, p3 #
function plane3pts(p1, p2, p3)
normal = LinearAlgebra.cross(p2 - p1, p3 - p1)
if LinearAlgebra.norm_sqr(normal) == 0
error("The three points are colinear.")
end
offset = LinearAlgebra.dot(Meshes.coordinates(p1), normal)
return (normal = normal, offset = offset)
end

# center, radius and normal of the circle passing by three points #
function circleCenterAndRadius(p1, p2, p3)
v12 = p2 - p1
v13 = p3 - p1
p12 = Meshes.coordinates(p1 + v12/2)
p13 = Meshes.coordinates(p1 + v13/2)
normal, offset = plane3pts(p1, p2, p3)
A = transpose(hcat(normal, v12, v13))
b = [
offset,
LinearAlgebra.dot(p12, v12),
LinearAlgebra.dot(p13, v13)
]
center = Meshes.Point(inv(A) * b)
r = LinearAlgebra.norm2(p1 - center)
return (center = center, radius = r, normal = normal)
end

# key transformation #
function keyTransform(p1, p2, p3)
center, radius, normal = circleCenterAndRadius(p1, p2, p3)
normal = normal / LinearAlgebra.norm2(normal)
s = sqrt(normal[1]^2 + normal[2]^2)
if s == 0
return (matrix = I, center = center, radius = radius)
end
u = [normal[2] / s, -normal[1] / s, 0]
v = LinearAlgebra.cross(normal, u)
R = hcat(u, v, normal)
return (rotMatrix = R, center = center, radius = radius)
end

"""
torusMesh(R, r; nu = 50, nv = 30)

Mesh of a torus with major radius R and minor radius r,
with the z-axis as its axis of rotation and the xy-plane as its
equatorial plane.

# Arguments
- R: major radius
- r: minor radius
- nu, nv: numbers of subdivisions

"""
function torusMesh(R, r; nu = 50, nv = 30)
kxy = R^2 - r^2
kz = sqrt(kxy) * r
s = sqrt(kxy) / r
u_ = LinRange(-s*pi, s*pi, nu)
v_ = LinRange(-pi, pi, nv)
points = [Meshes.Point(
[kxy * sin(u/s), kxy * cos(u/s), kz * sin(v)] / (R - r*cos(v))
) for u in u_ for v in v_
]
topo = Meshes.GridTopology((nv, nu), (true, true))
return Meshes.SimpleMesh(points, topo)
end

"""
torusMesh(r, p1, p2, p3; nu = 50, nv = 30)

Mesh of a torus whose equator passes through three given points.

# Arguments
- r: minor radius
- p1, p2, p3: the three points
- nu, nv: numbers of subdivisions

"""
function torusMesh(r, p1, p2, p3; nu = 50, nv = 30)
M, center, R = keyTransform(p1, p2, p3)
mesh = torusMesh(R, r; nu = nu, nv = nv)
vertices = [
center + Meshes.Vec3(M * Meshes.coordinates(v))
for v in Meshes.vertices(mesh)
]
return Meshes.SimpleMesh(vertices, Meshes.topology(mesh))
end

And an example:

using MeshViz

a = 1 / sqrt(3)
p1 = [ a, -a, -a]
p2 = [ a,  a,  a]
p3 = [-a, -a,  a]
p4 = [-a,  a, -a]

mesh1 = torusMesh(0.1, p1, p2, p3)
mesh2 = torusMesh(0.1, p1, p2, p4)
mesh3 = torusMesh(0.1, p1, p3, p4)
mesh4 = torusMesh(0.1, p2, p3, p4)

function draw()
MeshViz.viz!(mesh2; color = :red)
MeshViz.viz!(mesh3; color = :blue)
MeshViz.viz!(mesh4; color = :green)
end

MeshViz.viz(mesh1; color=:yellow)
draw()