Want to share your content on R-bloggers? click here if you have a blog, or here if you don't.

I wanted to plot an algebraic isosurface with POV-Ray but the expression of the polynomial defining the isosurface was very long (the polynomial had degree 12). Moreover there was a square root in the coefficients ($$\sqrt{3}$$) as well as $$\cos t$$ and $$\sin t$$, where $$t$$ is a parameter I wanted to vary in order to make an animation. So I needed a tool able to expand a polynomial with some literal values in the coefficients. This is not possible with the Ryacas package.

I finally found this tool: the caracas package. It allows to use the Python library SymPy in R. I didn’t carefully read its documentation yet, I don’t know whether it has other features. But this feature is a great one.

Here is a small example:

library(caracas)
def_sym(x, y, z, a, b) # symbolic values
poly <- sympy_func(
x^2 + a*x^2 + 2/3*y + b*y + x*z + a*x*z, "Poly", domain = "QQ[a,b]"
)
as.character(poly)

This gives:

"Poly((a + 1)*x^2 + (a + 1)*x*z + (b + 2/3)*y, x, y, z, domain='QQ[a,b]')"

That is great. Here QQ[a,b] is the field $$\mathbb{Q}[a,b]$$. I lost a significant part of my knowledge in mathematics but I think this is a field. It doesn’t matter. Roughly speaking, this is the set of rational numbers to which we add the two elements $$a$$ and $$b$$. So there are treated as constants, as if they were some numbers.

To get a coefficient, for example the one of $$xz = x^1y^0z^1$$:

sympy <- get_sympy()
sympy$Poly$nth(poly$pyobj, 1L, 0L, 1L) This gives: a + 1 Everything needed for writing the POV-Ray code was there. I wrote a small script in addition to generate this code. I show it below with the above small example: library(caracas) library(partitions) # to get the compositions of an integer, # representing the degrees with a given total def_sym(x, y, z, a, b) poly <- sympy_func( x^2 + a*x^2 + 2/3*y + b*y + x*z + a*x*z, "Poly", domain = "QQ[a,b]" ) sympy <- get_sympy() f <- function(comp){ xyz <- sprintf("xyz(%s): ", toString(comp)) coef <- sympy$Poly$nth(poly$pyobj, comp[1L], comp[2L], comp[3L])
if(coef == 0) return(NULL)
paste0(xyz, coef, ",")
}
for(deg in 0L:2L){
comps <- compositions(deg, 3L)
povray <- apply(comps, 2L, f, simplify = FALSE)
cat(
unlist(povray), sep = "\n", file = "povray.txt", append = deg > 0L
)
}

And here is the povray.txt file generated by this script:

xyz(0, 1, 0): b + 2/3,
xyz(2, 0, 0): a + 1,
xyz(1, 0, 1): a + 1,

One just has to remove the trailing comma, and this the desired POV-Ray code.

I won’t leave you without showing the animation:

Credit to ‘ICN5D’ for the isosurface.