In my attempt to learn Python and its mathematical packages, I have translated a library for generating Povray plots in Euler Math Toolbox into Python. The module is not polished yet, but if you are interested, ask in the comments for a Beta version.

Have a look at this tutorial for Euler Math Toolbox to see what can be done with this interface to Povray. The code for the knot in EMT is here.

For a start, let us generate a simple plot with an axis. I am going to explain below what’s happening here.

from povray import *
setengine("pvengine")

import numpy as np

def f(x,y):
    return x**y

povstart(zoom=4,angle=rad(210),viewcenter=[0.5,0.5,0.5])
t=np.arange(0,1,0.01)
fg = graph(f,t,t,look=look(green))
writeln(intersection([fg,cylinder([0,0,0],[0,0,1],1)]))
for i in range(1,4):
    writeln(axis(0,1,axis=i))
povend()

The Povray engine I am using is a version of Povray for Windows (see here). The name of the executable is „povengine.exe“. It is in the search Path of Windows. For Mac M1, I found a „homebrew“ version of a command line tool named „povray“. It is not fully supported. (These problems are one reason I abandoned Macs for math.) In any case, you can set the correct program to start Povray with „setengine“.

The function „povstart“ opens an output file, by default „povout.pov“ and writes the necessary headers and includes to this file. It does also set the default view to the scene. You can change the light here too.

We then create x,y,z-matrices containing the coordinates of the surface inside the function „graph“. It uses „meshgrid“ internally and assumes that f can be vectorized.

The „graph“ function writes the scene to the file and embeds it in an object of Povray, returning the object. The scene consists of 4 triangles for each rectangle in the grid. The center of gravity is used as one corner in each of the four triangles. Povray needs a list of corners and a list of triangles, and „graph“ writes both lists to the file.

The intersection() function shows the way that this module works: It does not write to the file, but returns a Povray string representing the object, in this case the intersection of two other objects, both represented as strings. The cylinder() function is one of them, the other our graph. Here is, how the Povray strings typically look like.

cylinder { <0,0,0>, <0,0,1>, 1
    texture { pigment { color rgbf <0.470588,0.470588,0.470588> }  } 
    finish { ambient 0.2 }
}

That’s one string with lines separated by „\n“.

The alternative to this approach was to generate Povray objects and use these objects. But I found the string approach more appealing. However, for the surface, the string would be huge. Thus, the object is written to the file directly.

The axis() function creates an axis, as follows. The texture can be changed using a „look“ parameter for each object.

union {
  cylinder { <-0.1,0,0>,<1.1,0,0>,0.02 }
  cone { 
    <1.22,0,0>,0
    <1.1,0,0>,0.08
  }
  texture { pigment { color rgbf <0.470588,0.470588,0.470588> } }
}

Now, let us try the knot.

from sympy import *
import numpy as np

t = symbols('t')
fx = 4*(1+sin(3*t)/3)*sin(2*t)
fy = 4*(1+sin(3*t)/3)*cos(2*t)
fz = 2*cos(3*t)

vt = np.linspace(-np.pi,np.pi,1000)
vs = np.linspace(-np.pi,np.pi,160)

f = [fx,fy,fz]
ffunc = lambdify(t,f)
df = [diff(f[k],t) for k in range(3)]
dffunc = lambdify(t,df)

w1 = [df[1],-df[0],0]
w1func = lambdify(t,w1)

def crossproduct (v,w):
    return np.array([v[1]*w[2]-v[2]*w[1],-v[0]*w[2]+v[2]*w[0],v[0]*w[1]-v[1]*w[0]])
def norm(v):
    return np.sqrt(v[0]*v[0]+v[1]*v[1]+v[2]*v[2])

def w1n(t):
    h = w1func(t)
    return h/norm(h)
def w2n(t):
    h = crossproduct(dffunc(t),w1func(t))
    return h/norm(h)

def S(t,s):
    return ffunc(t)+np.cos(s)*w1n(t)+np.sin(s)*w2n(t)
x = np.zeros((len(vt),len(vs)))
y = np.zeros((len(vt),len(vs)))
z = np.zeros((len(vt),len(vs)))
for i in range(len(vt)):
    for j in range(len(vs)):
        u = S(vt[i],vs[j])
        x[i,j] = u[0]
        y[i,j] = u[1]
        z[i,j] = u[2]

import povray as pv

pv.povstart(distance=15)
surf = pv.xyz_surface(x,y,z,look=pv.look(pv.gray,phong=1))
pv.writeln(surf)
pv.povend(engine="pvengine")

This needs a lot of explaining. The center curve of the knot is a curve in space parameterized by „fx“, „fy“ and „fz“. I did some tests first, to make sure that this is indeed working, using a 3D plot in Matplotlib. The result was the following. It is a different style from the 3D plots you get in Euler Math Toolbox.

Then I build a tube around that curve. To achieve that, the expression „df“ is a vector representing the derivative of „f“, i.e., the vector pointing in the direction of the curve at each point. Then „w1“ is a vector perpendicular to „df“.

To get another perpendicular vector „w2“, I use the crossproduct of „df“ and „w1“. Both, „w1“ and „w2“ are then normalized to have length 1. However, this is already done numerically. To parameterize the tube, I use

\(S(t,s) = f(t) + \cos(s)w_1(t) + \sin(s)w_2(t)\)

This function is mapped to the all s and t, using a double loop.