1;
clear all
n=50;
etamax=pi/2;
etamin=-pi/2;
wmax=pi;
wmin=-pi;
deta=(etamax-etamin)/n;
dw=(wmax-wmin)/n;
a=[3/4 3/4 1];
epsilon=[1.8 1.8]
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
i=1;
endfor
xlabel ('x')
ylabel ('y')
zlabel ('z')
title(num2str(epsilon))


#gset terminal x11;
mesh(x,y,z)
print("spq.ps","-dps;-mono")

#gset terminal x11;