NAME
    gumdrop_torus

FUNCTION
    gumdrop_torus(step)

NOTES
    Finds the iso-surface of the equation:

        4 * [x^4 + (y^2 + z^2)^2] +
        17 * x^2 * (y^2 + z^2) -
        20 * (x^2 + y^2 + z^2) + 17 = 0 
        
    Refer to Paul Bourke' home page (http://astronomy.swin.edu.au/~pbourke/surfaces/gumdrop/).
    
    Example:
        require("plot_simple")
        plot = plot_simple.new()
        plot:add_static(zeGrf.new("light"))
        require("gumdrop_torus")
        xyz, nor = gumdrop_torus(0.05)
        xyz:scale(80, 80, 80)
        nor:scale(-1, -1, -1)
        shape = zeGrf.new("polygon")
        shape:set{type = "triangles", vertex = xyz,
                  vertex_normal = nor, color = {0, .5, .5, 1}}
        plot:add(shape)
        require("custom_materials")
        mat = custom_materials.new()
        mat:pewter()
        plot:add_static(mat.material)
        plot:animate()

INPUTS
    step  - calculation step, should be much smaller than 1

OUTPUTS
    Two zeVertex objects containng trianle vertices and their normal.

SOURCE

require("iso_surface")

function gumdrop_torus(step)
    assert(step > 0)
    assert(step < 1)

    local function sfunc(x, y, z)
        local x2, y2, z2 = x*x, y*y, z*z
        return 4*(x2*x2 + (y2 + z2)*(y2 + z2)) +
               17*x2*(y2 + z2) -
               20*(x2 + y2 + z2) + 17
    end
    
    local function nfunc(x, y, z)
        local x2, y2, z2 = x*x, y*y, z*z
        local nx = 16*x*x2 + 34*x*(y2 + z2) - 40*x
        local ny = 16*y*(y2 + z2) + 34*x2*y - 40*y
        local nz = 16*z*(y2 + z2) + 34*x2*z - 40*z
        local d = math.sqrt(nx*nx + ny*ny + nz*nz)
        return -nx/d, -ny/d, -nz/d
    end
    
    local xarr, yarr, zarr = zeUtl.new("double", "double", "double")
    local n = 2 / step
    xarr:range(-2, step, 2*n+1)
    yarr:range(-2, step, 2*n+1)
    zarr:range(-2, step, 2*n+1)
    
    local xyz, nor = iso_surface(0, xarr, yarr, zarr, sfunc)
    n = xyz:size()-1
    nor:clear()
    for k = 0, n do
        nor:add(nfunc(xyz:get(k)))
    end
    
    return xyz, nor
end