File:Algebraic surfaces.png

using GLMakie
using Meshing, GeometryBasics

set_theme!(theme_black())

isoval = 100
algo = MarchingCubes(iso=isoval, insidepositive=false)

function show_isosurface(f,h,ξ)
  s = [h(x,y,z) for x in ξ, y in ξ, z in ξ] .+ isoval

    mc = GeometryBasics.Mesh(s, algo)
    mesh(f,normal_mesh(mc), 
    color=(:dodgerblue, .7), 
    diffuse = Vec3f0(0.8),
    specular = Vec3f0(1.1),
    shininess = 30f0,
    backlight = 5f0,
    axis = (show_axis = false, ))
end

f = Figure(size=(900, 900))
# isosurfaces for h(x,y,z)=0

## Torus 
ξ = -40:.5:40;
c = 20
a = 15
torus(x,y,z) = ((hypot(x,y)-c)^2+z^2-a^2)

show_isosurface(f[1,1],torus,ξ);

## Entzensberger star
Xm = 1
ξ = -Xm:.01:Xm;
estar(x,y,z) = 100*(x^2*y^2 + y^2*z^2 + x^2*z^2) - (1 - x^2 - y^2 - z^2)^3

show_isosurface(f[1,2],estar,ξ);


## Tetrahedron
Xm = 10
ξ = -Xm:.05:Xm;
tetra(x,y,z) = x^4 + 2*x^2*y^2 + 2*x^2*z^2 + y^4 + 2*y^2*z^2 + z^4 + 8*x*y*z - 10*x^2 - 10*y^2 - 10*z^2 + 20

show_isosurface(f[2,1],tetra,ξ);


## Decocube 
Xm = 2.5
ξ = -Xm:.025:Xm;
b = 1
c = 2.2
t = 1.2

deco(x,y,z) = ((x^2 + y^2 - c^2)^2 + (z - 1)^2*(z + 1)^2)*((y^2 + z^2 - c^2)^2 +
(x - 1)^2*(x + 1)^2)*((z^2 + x^2 - c^2)^2 + (y - 1)^2*(y + 1)^2) -
t*(1 + b*(x^2 + y^2 + z^2))

show_isosurface(f[2,2],deco,ξ);