function [area_elipse,center_distances,coeff] = Error_elipse3D(Scores,Classes,Clr)
    [cUnq,~,cInd] = unique(Classes);
    numP = size(cUnq,1);
    
    area_elipse = zeros(1,numP);
    lim = xlim;
    for k = 1:numP
        % 2D area calculation
        % indices of points in this group
        idx = (Classes == k);
    
        % substract mean
        Mu = mean(Scores(idx,:));
        X0 = bsxfun(@minus, Scores(idx,:), Mu);
    
        STD = 3;                     % 2 standard deviations
        conf = 2*normcdf(STD)-1;     % covers around 95% of population
        scale = chi2inv(conf,2);     % inverse chi-squared with dof=#dimensions
        
        Cov = cov(X0) * scale;
        [V, D] = eig(Cov);
        [D, order] = sort(diag(D), 'descend');
        D = diag(D);
        V = V(:, order);
    
        t = linspace(0,2*pi,100);
        e = [cos(t) ; sin(t)];        % unit circle
        VV = V*sqrt(D);               % scale eigenvectors
        e = bsxfun(@plus, VV*e, Mu'); % project circle back to orig space

        % center ellipse error
        center_classes(k,1:2) = Mu';

        % plot cov and major/minor axes
        fill = patch(e(1,:), e(2,:), Clr(k), 'EdgeColor', 'none');
        area_elipse(k) = (polyarea(e(1,:), e(2,:)))/(lim(2)-lim(1));
        alpha(fill,0.2)
        text(center_classes(k,1)*1.1,center_classes(k,2)*1.1,num2str(area_elipse(k),3),'Color','r','FontWeight','bold')
        scatter(center_classes(k,1),center_classes(k,2),'h','filled','r')


        % 3D area calculation
%         % indices of points in this group
%         idx = (Classes == k);
%     
%         % substract mean
%         Mu = mean(Scores(idx,:));
%         X0 = bsxfun(@minus, Scores(idx,:), Mu);
%     
%         STD = 2;                     % 2 standard deviations
%         conf = 2*normcdf(STD)-1;     % covers around 95% of population
%         scale = chi2inv(conf,2);     % inverse chi-squared with dof=#dimensions
%         
%         Cov = cov(X0) * scale;
% 
%         [U,L] = eig(Cov);
%         
%         % 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*sqrt(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,:)+Mu(1); 
%         y = data(n+1:2*n,:)+Mu(2); 
%         z = data(2*n+1:end,:)+Mu(3);
% 
%         % center ellipse error
%         center_classes(k,1:3) = Mu';
%      
%         fill = surf(x,y,z,'FaceAlpha',0.1);
%         fill.EdgeColor = convertCharsToStrings(Clr(k));
%         fill.EdgeAlpha = 0.5;
%         fill.MeshStyle = 'column';
%         fill.FaceColor = convertCharsToStrings(Clr(k)); 
% 
%         % plot cov and major/minor axes
%         area_elipse(k) = trapz(x,trapz(y,z));
%         text(center_classes(k,1)*1.1,center_classes(k,2)*1.1,num2str(area_elipse(k),3),'Color','r','FontWeight','bold')
%         scatter(center_classes(k,1),center_classes(k,2),'h','filled','r')    

    end    
    
    center_distances = pdist(center_classes,'euclidean');
    
    % This algorithm only works for models with three classes
    for n = 1:numP
        obj_line = line([center_classes(1,1) center_classes(n+1,1)],[center_classes(1,2) center_classes(n+1,2)],'color',[0 0 0],'linestyle','--');
        
        P1 = center_classes(1,1:2);
        P2 = center_classes(n+1,1:2);
        mid_pt= (P1(:) + P2(:)).'/2;

        text(mid_pt(1,1),mid_pt(1,2),num2str(center_distances(n),3),'FontWeight','bold')

        coeff(n) = center_distances(n)/(area_elipse(1,1) + area_elipse(1,n+1));

        if n == (numP-1)
            obj_line = line([center_classes(n,1) center_classes(n+1,1)],[center_classes(n,2) center_classes(n+1,2)],'color',[0 0 0],'linestyle','--');
            
            P1 = center_classes(n,1:2);
            P2 = center_classes(n+1,1:2);
            mid_pt= (P1(:) + P2(:)).'/2;
    
            text(mid_pt(1,1),mid_pt(1,2),num2str(center_distances(n+1),3),'FontWeight','bold')
            coeff(n+1) = center_distances(n+1)/(area_elipse(1,n+1) + area_elipse(1,n+1));
            return
        end
    end  
end


