function draw_elipsoid(C)


[U,L] = eig(C);

% For N standard deviations spread of data, the radii of the eliipsoid will
% be given by N*SQRT(eigenvalues).

N = 1; % choose your own N
radii = N*(diag(L));

% generate data for "unrotated" ellipsoid
[xc,yc,zc] = ellipsoid(0,0,0,radii(1),radii(2),radii(3));

% rotate data with orientation matrix U and center M
a = kron(U(:,1),xc);
b = kron(U(:,2),yc);
c = kron(U(:,3),zc);

data = a+b+c; n = size(data,2);
x = data(1:n,:); y = data(n+1:2*n,:); z = data(2*n+1:end,:);

% now plot the rotated ellipse
figure;
sc = surf(x,y,z);
set(sc,'FaceColor',[1 0 0],'FaceAlpha',0.5, 'EdgeAlpha', 0,'edgecolor','none')
%shading interp

xlabel('X');
ylabel('Y');
zlabel('Z');

xlim([-0.002,0.002])
ylim([-0.002,0.002])
zlim([-0.002,0.002])


