function retval=superellipse(epsilon,a)
  n=50;
  etamax=pi/2;
  etamin=-pi/2;
  wmax=pi;
  wmin=-pi;
  deta=(etamax-etamin)/n;
  dw=(wmax-wmin)/n;
  k=0;
  l=0;
  for i=1:n+1
    eta(i)=etamin+(i-1)*deta;
    for j=1:n+1
      w(j)=wmin+(j-1)*dw;
      x(i,j)=a(1)*sign(cos(eta(i)))*abs(cos(eta(i)))^epsilon(1)*sign(cos(w(j)))*abs(cos(w(j)))^epsilon(2);
      y(i,j)=a(2)*sign(cos(eta(i)))*abs(cos(eta(i)))^epsilon(1)*sign(sin(w(j)))*abs(sin(w(j)))^epsilon(2);
      z(i,j)=a(3)*sign(sin(eta(i)))*abs(sin(eta(i)))^epsilon(1);
    endfor;
  endfor;
  mesh(x,y,z);
  return;
endfunction;