NAME
    super_surfaces

FUNCTION
    super_surfaces(n, a, b, m, n1, n2, n3)

NOTES
    Generate super surfaces based on parametric equations of
    
        x = r(v) cos(v) r(u) cos(u)
        y = r(v) sin(v) r(u) cos(u)
        z = r(u) sin(u)
        r(v) = [|a cos(m v / 4)|^n1 + |b sin(m v / 4)|^n2]^(-n3)
        
    Refer to Paul Bourke's home page at http://astronomy.swin.edu.au/~pbourke/surfaces/supershape3d/
    
    Example:
        require("plot_simple")
        plot = plot_simple.new()
        plot:add_static(zeGrf.new("light"))
        require("super_surfaces")
        xyz, nor = super_surfaces(36, 1, 1, 1, 1, 1, 1)
        xyz:scale(100, 100, 100)
        shape = zeGrf.new("polygon")
        shape:set{type = "quads", vertex = xyz,
                  vertex_normal = nor, color = {0, .7, .7, 1}}
        plot:add(shape)
        plot:animate()
NPUTS
    n - number of segments in 0 to 2pi
    a, b, m, n1, n2, n3 - shape parameters

OUTPUTS
    Two zeVertex objects containing coordinates and normals.

SOURCE

require("surface_generator")

function super_surfaces(n, a, b, m, n1, n2, n3)
    assert(n > 16)
    assert(a ~= 0)
    assert(b ~= 0)
    assert(m > 0)
    assert(n1 > 0)
    assert(n2 > 0)
    assert(n3 > 0)

    local function radius(v)
        return 1 / math.pow(math.pow(math.abs(a*math.cos(m*v/4)), n1) + 
                            math.pow(math.abs(b*math.sin(m*v/4)), n2), n3)
    end
    
    local function sfunc(u, v)
        local ru, rv = radius(u), radius(v)
        return rv*math.cos(v)*ru*math.cos(u),
               rv*math.sin(v)*ru*math.cos(u),
               ru*math.sin(u)
    end

    local function nfunc(u, v)
        local d = 1.e-5
        local x, y, z = sfunc(u, v)
        local xu, yu, zu, xv, yv, zv
        if math.abs(u-math.pi/2) > d then
            xu, yu, zu = sfunc(u+d, v)
            xu, yu, zu = xu-x, yu-y, zu-z
        else
            xu, yu, zu = sfunc(u-d, v)
            xu, yu, zu = x-xu, y-yu, z-zu
        end
        if math.abs(v-math.pi) > d then
            xv, yv, zv = sfunc(u, v+d)
            xv, yv, zv = xv-x, yv-y, zv-z
        else
            xv, yv, zv = sfunc(u, v-d)
            xv, yv, zv = x-xv, y-yv, z-zv
        end
        return zeMake.normal2(0, 0, 0, xu, yu, zu, xv, yv, zv)
    end

    local U, V = zeUtl.new("double", "double")
    U:range(-math.pi/2, math.pi/n, n+1)
    V:range(-math.pi, math.pi/n, 2*n+1)
    
    local xyz = surface_generator(U, V, sfunc)
    require("check_normal")
    local nor = check_normal(xyz, surface_generator(U, V, nfunc))

    return xyz, nor
end