// Assumes Lmin, Lmax, Lmid, Emin, Emax, Emid are all defined
// given L and Es 
scf(0);
clf();
L = [Lmin, Lmid, Lmax];
a = 1 ./ sqrt(L);
E = [Emin, Emid, Emax];
//x  = a*cos(t)*sin(s);
//y  = b*sin(t)*sin(s);
//z  = c*cos(s)
theta = 0:%pi/50:2*%pi;
phi = 0:%pi/50:%pi;
n = size(theta,2);
m = size(phi,2);
lat = zeros(3,n);
long = zeros(3,m);
// latitude
// percent up or down
ulist=-0.5:0.5:0.5;
for j = 1:3,
u = ulist(j);
up = u;
left = sqrt(1-up^2);
for i = 1:n,
	t = theta(i);
	lat(:,i) = left*(E(:,2)*a(1,2)*cos(t) + E(:,3)*a(1,3)*sin(t))+up*a(1,1)*E(:,1);
end;
param3d1(lat(1,:)',lat(2,:)',list(lat(3,:)',[3]));
end;
// equator
for i = 1:n,
	t = theta(i);
	lat(:,i) = (E(:,2)*a(1,2)*cos(t) + E(:,3)*a(1,3)*sin(t));
end;
param3d(lat(1,:),lat(2,:),lat(3,:));
// longitude
tlist = 0:%pi/2:2*%pi;
for i = 1:4
t = tlist(i);
for j = 1:m;
	s = phi(j);
	long(:,j) = (E(:,2)*a(1,2)*cos(t) + E(:,3)*a(1,3)*sin(t))*sin(s)+E(:,1)*a(1,1)*cos(s);
end;
param3d1(long(1,:)',long(2,:)',list(long(3,:)',[5]));
end;
xx=[0 0 0;Emax(1,1)/sqrt(Lmax) Emid(1,1)/sqrt(Lmid) Emin(1,1)/sqrt(Lmin)];
yy=[0 0 0;Emax(2,1)/sqrt(Lmax) Emid(2,1)/sqrt(Lmid) Emin(2,1)/sqrt(Lmin)];
zz=[0 0 0;Emax(3,1)/sqrt(Lmax) Emid(3,1)/sqrt(Lmid) Emin(3,1)/sqrt(Lmin)];
param3d1(xx,yy,list(zz,[2,2,2]));
xtitle('Ellipsoid eigendirections in blue, equator in black');
graph=gca();
graph.isoview="on";