NAME
    super_ellipsoid

FUNCTION
    super_ellipsoid(slices, rx, ry, rz, p1, p2)

NOTES
    Creates 2D super ellipse line.

    Based on the equations
    
        x = rx * cos(v)^p1 * cos(u)^p2
        y = ry * cos(v)^p1 * sin(u)^p2
        z = rz * sin(v)^p1
        
    of Paul Bourke' home page (http://astronomy.swin.edu.au/~pbourke/surfaces/mobius/).
    Where 0 <= u, v <= pi / 2 and 0 < p1, p2 < infinite.
    
    The normal is
        nx = cos(v)^(2-p1) * cos(u)^(2-p2) / rx
        ny = cos(v)^(2-p1) * sin(u)^(2-p2) / ry
        nz = sin(v)^(2-p1) / rz

    Example: 
 NPUTS
    slices - slices between 0 and  pi/2
    rx     - size parameter in x-direction
    ry     - size parameter in y-direction
    rz     - size parameter in y-direction
    p1     - shape parameter
    p2     - shape parameter

OUTPUTS
    Two zeArray objects containing coordinates and normals.

SOURCE

require("register")

function super_ellipsoid(slices, rx, ry, rz, p1, p2)
    assert(slices >= 16)
    assert(rx > 1)
    assert(ry > 1)
    assert(rz > 1)
    assert(p1 > 0)
    assert(p2 > 0)
    
    local function vertex_xyz(p1, p2, rx, ry, rz, ucos, vcos, usin, vsin)
        vcos = vcos^p1
        ucos = ucos^p2
        vsin = vsin^p1
        usin = usin^p2
        return rx * vcos * ucos,
               ry * vcos * usin,
               rz * vsin
    end

    local function normal_xyz(p1, p2, rx, ry, rz, ucos, vcos, usin, vsin)
        vcos = vcos^(2-p1)
        ucos = ucos^(2-p2)
        vsin = vsin^(2-p1)
        usin = usin^(2-p2)
        return vcos * ucos / rx,
               vcos * usin / ry,
               vsin /rz
    end

    local function uv_vecs(slices, p)
        local d = 0.5 * math.pi / slices
        local cos, sin = zeUtl.new("double", "double")
        cos:range(0, d, slice+1)
        zeMath.cos(cos)
        zeMath.pow(cos, p)
        cos:setele(0, 0, 1)
        cos:setele(slices, 0, 0)
        sin:range(0, d, slice+1)
        zeMath.sin(sin)
        zeMath.pow(sin, p)
        sin:setele(0, 0, 0)
        sin:setele(slices, 0, 1)
        return cos, sin
    end
    
    local vec, cos1, cos2, sin1, sin1
         = zeUtl.new("double", "double", "double", "double", "double")
    cos1:range(0, d, slices+1)
    zeMath.cos(cos1)
    cos1:copy(cos2)
    zeMath(cos1,p

    local vertex, normal = zeGrf.new("vertex", "vertex")

    
    local x1, y1, z1 = rx, 0, 0
    local x2, y2, z2, x3, y3, z3, x4, y4, z4
    
    local nx1, ny1, nz1 = 1, 0, 0
    local nx2, ny2, nz2, nx3, ny3, nz3, nx4, ny4, nz4

    local vec, cos1, cos2, sin1, sin1
         = zeUtl.new("double", "double", "double", "double", "double")
    cos1:range(0, d, slices+1)
    zeMath.cos(cos1)
    cos1:copy(cos2)
    zeMath(cos1,p
            
    for i = 1, slices do
        vcos2, vsin2 = math.cos(i * d), math.sin(i * d)
        x4, y4, z4 = vertex_xyz(p1, p2, rx, ry, rz, ucos, vcos, usin, vsin)
        for j = 1, slices do
            uc2, us2 = math.cos(j * d), math.sin(j * d)
            x2 = rx * vc2^p1 * uc2^p2
            y2 = ry * vc2^p1 * us2^p2
            vertex:add(x1, y1, z1)
            vertex:add(x1, y1, z1)
        end
    end
    
    xyz:resize(slices, 3)
    nor:resize(slices, 3)

        x = rx * cos(theta)^p1 * cos(phi)^p2
        y = ry * cos(theta)^p1 * sin(phi)^p2
        z = rz * sin(theta)^p1

    cos:range(0, d, slices)
    zeMath.cos(cos)
    sin:range(0, d, slices)
    zeMath.sin(sin)
    
    cos:copy(vec)
    zeMath.pow(vec, p)
    vec:mul(rx)
    vec:setele(0, 0, rx)
    vec:setele(slices-1, 0, 0)
    xyz:setarr(0, vec)

    sin:copy(vec)
    zeMath.pow(vec, p)
    vec:mul(ry)
    vec:setele(0, 0, 0)
    vec:setele(slices-1, 0, ry)
    xyz:setarr(1, vec)

    sin:copy(vec)
    zeMath.pow(vec, p-1)
    vec:mul(ry)
    vec:mul(cos)
    vec:setele(0, 0, 1)
    vec:setele(slices-1, 0, 0)
    nor:setarr(0, vec)
    
    cos:copy(vec)
    zeMath.pow(vec, p-1)
    vec:mul(rx)
    vec:mul(sin)
    vec:setele(0, 0, 0)
    vec:setele(slices-1, 0, 1)
    nor:setarr(1, vec)
    
    zeMake.normalize(nor)
        
    return xyz, nor
end