NAME
    bifolia_surface

FUNCTION
    bifolia_surface(step)

NOTES
    Finds the iso-surface of the equation:
        [x^2 + y^2 + z^2]^2 - 3*y*[x^2 + z^2] = 0
        
    Refer to Paul Bourke' home page (http://astronomy.swin.edu.au/~pbourke/surfaces/gumdrop/)
    or http://www.uib.no/People/nfytn/surfaces.htm.
    
    Example:
        require("plot_simple")
        plot = plot_simple.new()
        plot:add_static(zeGrf.new("light"))
        require("bifolia_surface")
        xyz, nor = bifolia_surface(0.05)
        xyz:scale(150, 150, 150)
        nor:scale(-1, -1, -1)
        shape = zeGrf.new("polygon")
        shape:set{type = "triangles", vertex = xyz,
                  vertex_normal = nor, color = {0, .7, .7, 1}}
        plot:add(shape)
        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 bifolia_surface(step)
    assert(step > 0)
    assert(step < 1)

    local function sfunc(x, y, z)
        local x2, y2, z2 = x*x, y*y, z*z
        return (x2 + y2 + z2)^2 - 3*y*(x2 + z2)
    end
    
    local function nfunc(x, y, z)
        local x2, y2, z2 = x*x, y*y, z*z
        local nx = 4*x*(x2 + y2 + z2) - 6*x*y
        local ny = 4*y*(x2 + y2 + z2) - 3*(x2 + z2)
        local nz = 4*z*(x2 + y2 + z2) - 6*z*y
        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 = 1 / step
    xarr:range(-1, step, 2*n+1)
    yarr:range(-0.01, step, n+1)
    zarr:range(-1, 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