% The definition of the superquadrics are taken from the paper
%       Gordon Kindlmann. Superquadric Tensor Glyphs.
%       In Proceedings IEEE TVCG/EG Symposium on Visualization 2004,
%         pages 147-154, May 2004, http://www.cs.utah.edu/~gk/papers/vissym04/
% Formulas (3) and (5) in 
% http://www.cs.utah.edu/~gk/papers/vissym04/vissym04kindlmann.pdf
%
% In the paper the definition range for phi is [0,pi] and for theta [0,2pi].
% The present implementation in VTK uses different ranges: phi in [-pi/2,pi], theta in [-pi,pi]
% The formulas have been adapted accordingly.

function retval=superquadricZ(alpha,beta,scaling)
    n=50;
    phimin=-pi/2;
    phimax=pi/2;
    thetamin=-pi;
    thetamax=pi;
    dphi=(phimax-phimin)/n;
    dtheta=(thetamax-thetamin)/n;
    [i,j] = meshgrid(1:n+1,1:n+1);
    phi   = phimin + (i-1) * dphi;
    theta = thetamin + (j-1) * dtheta;
    x = scaling(1) .* sign(cos(theta)) .* abs(cos(theta)).^alpha .* sign(cos(phi)) .* abs(cos(phi)).^beta;
    y = scaling(2) .* sign(sin(theta)) .* abs(sin(theta)).^alpha .* sign(cos(phi)) .* abs(cos(phi)).^beta;
    z = scaling(3) .* sign(sin(phi))   .* abs(sin(phi)).^beta;

    mesh(x,y,z);
    xlabel("x");
    ylabel("y");
    zlabel("z");
    endfunction;

function retval=superquadricX(alpha,beta,scaling)
    n=50;
    phimin=-pi/2;
    phimax=pi/2;
    thetamin=-pi;
    thetamax=pi;
    dphi=(phimax-phimin)/n;
    dtheta=(thetamax-thetamin)/n;
    [i,j] = meshgrid(1:n+1,1:n+1);
    phi   = phimin + (i-1) * dphi;
    theta = thetamin + (j-1) * dtheta;
    z =  scaling(1) .* sign(cos(theta)) .* abs(cos(theta)).^alpha .* sign(cos(phi)) .* abs(cos(phi)).^beta;
    y = -scaling(2) .* sign(sin(theta)) .* abs(sin(theta)).^alpha .* sign(cos(phi)) .* abs(cos(phi)).^beta;
    x =  scaling(3) .* sign(sin(phi))   .* abs(sin(phi)).^beta;

    mesh(x,y,z);
    xlabel("x");
    ylabel("y");
    zlabel("z");
    endfunction;

function retval=superquadricY(alpha,beta,scaling)
    n=50;
    phimin=-pi/2;
    phimax=pi/2;
    thetamin=-pi;
    thetamax=pi;
    dphi=(phimax-phimin)/n;
    dtheta=(thetamax-thetamin)/n;
    [i,j] = meshgrid(1:n+1,1:n+1);
    phi   = phimin + (i-1) * dphi;
    theta = thetamin + (j-1) * dtheta;
    x = -scaling(1) .* sign(cos(theta)) .* abs(cos(theta)).^alpha .* sign(cos(phi)) .* abs(cos(phi)).^beta;
    z =  scaling(2) .* sign(sin(theta)) .* abs(sin(theta)).^alpha .* sign(cos(phi)) .* abs(cos(phi)).^beta;
    y =  scaling(3) .* sign(sin(phi))   .* abs(sin(phi)).^beta;

    mesh(x,y,z);
    xlabel("x");
    ylabel("y");
    zlabel("z");
    endfunction;

alpha=0.25; beta=4; scale=[1 1 1]; superquadricX(alpha,beta,scale);
fprintf (stderr, "press any key to continue...\n"); pause;
alpha=0.25; beta=4; scale=[1 1 1]; superquadricY(alpha,beta,scale);
fprintf (stderr, "press any key to continue...\n"); pause;
alpha=0.25; beta=4; scale=[1 1 1]; superquadricZ(alpha,beta,scale);

