function plotcovar(u,v,C,m1,m2)
%get from http://www.visiondummy.com/2014/04/draw-error-ellipse-representing-covariance-matrix/
% Get the 95% confidence interval error ellipse
chisquare_val = 2.4477;
%INFO: pour une covariance=identité, trace un cercle de RAYON chisquare_val
%pour une covariance diagonale, les termes de variance sont les carrés des écarts type
%l'ellipse tracée un aura ses DEMI-axes de longeur ecart type * chisquare_val



%hold on
%plotcovar(10,20,[2,4;3,1],'r+','b-');
%plotcovar(10,20,[1,0;0,1],'r+','b-');

plot(u,v,m1)

%bvdp, added step to inforce that C is symetric 
%if this precaution is not taken, this kind of matrix can be provided 
%C=[0.999999999999951 2.73114864057789e-14;-1.50990331349021e-14 0.999999999999964]
%and it will result in complex eigen value, which won't be treatable by the
%atan2 function used later...
Cm=(C(2,1)+C(1,2))/2;
C(2,1)=Cm;
C(1,2)=Cm;


% Calculate the eigenvectors and eigenvalues
covariance =C;
[eigenvec, eigenval ] = eig(covariance);
[eigenvec, eigenval ] = eig(covariance,'nobalance');

% Get the index of the largest eigenvector
[largest_eigenvec_ind_c, r] = find(eigenval == max(max(eigenval)));
largest_eigenvec = eigenvec(:, largest_eigenvec_ind_c);

% Get the largest eigenvalue
largest_eigenval = max(max(eigenval));

% Get the smallest eigenvector and eigenvalue
if(largest_eigenvec_ind_c == 1)
    smallest_eigenval = max(eigenval(:,2));
    smallest_eigenvec = eigenvec(:,2);
else
    smallest_eigenval = max(eigenval(:,1));
    smallest_eigenvec = eigenvec(1,:);
end

% Calculate the angle scriptshellbetween the x-axis and the largest eigenvector
angle = atan2(largest_eigenvec(2), largest_eigenvec(1));

% This angle is between -pi and pi.
% Let's shift it such that the angle is between 0 and 2pi
if(angle < 0)
    angle = angle + 2*pi;
end

% Get the coordinates of the data mean
avg = [u;v]; %mean(data);

theta_grid = linspace(0,2*pi);
phi = angle;
X0=avg(1);
Y0=avg(2);
a=chisquare_val*sqrt(largest_eigenval);
b=chisquare_val*sqrt(smallest_eigenval);

% the ellipse in x and y coordinates 
ellipse_x_r  = a*cos( theta_grid );
ellipse_y_r  = b*sin( theta_grid );

%Define a rotation matrix
R = [ cos(phi) sin(phi); -sin(phi) cos(phi) ];

%let's rotate the ellipse to some angle phi
r_ellipse = [ellipse_x_r;ellipse_y_r]' * R;

% Draw the error ellipse
plot(r_ellipse(:,1) + X0,r_ellipse(:,2) + Y0,m2)
%hold on;

% Plot the original data
%plot(data(:,1), data(:,2), '.');
%mindata = min(min(data));
%maxdata = max(max(data));
%Xlim([mindata-3, maxdata+3]);
%Ylim([mindata-3, maxdata+3]);
%hold on;

% Plot the eigenvectors
%quiver(X0, Y0, largest_eigenvec(1)*sqrt(largest_eigenval), largest_eigenvec(2)*sqrt(largest_eigenval), '-m', 'LineWidth',2);
%quiver(X0, Y0, smallest_eigenvec(1)*sqrt(smallest_eigenval), smallest_eigenvec(2)*sqrt(smallest_eigenval), '-g', 'LineWidth',2);
%hold on;

% Set the axis labels
%hXLabel = xlabel('x');
%hYLabel = ylabel('y');