function daenmenB(aa,T,Q)
  n=15;
  p=0:pi/n:2*pi;
  t=(0:pi/n:pi)';
  ct=cos(t);st=sin(t);cp=cos(p);sp=sin(p);
  a1=1/sqrt(aa(1));a2=1/sqrt(aa(2));a3=1/sqrt(aa(3));
  xxx=a1*st*cp;
  yyy=a2*st*sp;
  yoko=ones(1,2*n+1);
  zzz=a3*ct*yoko;
  for i=1:n+1
    xxxx=T*Q*[xxx(i,:);yyy(i,:);zzz(i,:);yoko];
    x(i,:)=xxxx(1,:);y(i,:)=xxxx(2,:);z(i,:)=xxxx(3,:);
  end
  mesh(x,y,z)
  axis equal
  grid on
end

ax=1;ay=3;az=3;hxy=0;hyz=-1;hzx=0;bx=1;by=-1;bz=3;c=2; % ellipsoid
A=[ax hxy hzx;hxy ay hyz;hzx hyz az];
b=[bx;by;bz];
B=[A b;b' c];
[P AA]=eig(A);
aa=diag(AA);
x0=A\b;
T=[eye(3) -x0;zeros(1,3) 1];
Q=[P zeros(3,1);zeros(1,3) 1];
BBB=Q'*T'*B*T*Q;
aa=-aa/BBB(4,4);
daenmenB(aa,T,Q)