--- title: "Drawing a torus with Julia" author: "Stéphane Laurent" date: '2022-12-15' tags: julia, geometry, graphics output: md_document: variant: markdown preserve_yaml: true html_document: highlight: kate keep_md: no highlighter: pandoc-solarized --- I already showed how to draw a torus whose equator passes through three given points with [three.js](https://laustep.github.io/stlahblog/posts/threejsTorus.html), [Asymptote](https://gist.github.com/stla/f461877b9426df53d2db2d18db93a2e6), [R](https://laustep.github.io/stlahblog/posts/rglTorus.html), [POV-Ray](https://laustep.github.io/stlahblog/posts/povrayTorus.html), and [Haskell](https://laustep.github.io/stlahblog/posts/HaskellTorus.html). [Here](https://github.com/stla/PyTorusThreePoints) 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: ```julia 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: ```julia 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() ``` ![](./figures/julia_tori.png){width="50%"}