function [MAC_data,EMA_FEA] = MacW2_mod(fig,Phi_l,NodesPhi_l,W2_l,Phi_r,NodesPhi_r,W2_r,k)

%MacW2 calculates the MAC index of correlation between 2 sets of 
%eigenvectors weighted with the corresponding eigenvalue mismatch
%   INPUT:
%       fig        = figure number (0 per deactivate plotting)
%
%       Phi_l      = first set (left) of eigenvectors (n_l x r_l)
%       NodesPhi_l = correspondence table of the first set of eigenvectors
%                    (n_l x 2)
%       W2_l       = first set of eigenvalues (r_l x 1)
%
%       Phi_r      = second set (right) of eigenvectors (n_r x r_r)
%       NodesPhi_r = correspondence table of the second set of eigenvectors
%                    (n_r x 2)
%       W2_r       = second set of eigenvalues (r_r x 1)
%       k          = weight of the error on the eigenvalues with respect
%                    to the error on the eigenvectors
%
%   OUTPUT:
%       MAC_data   = values of MAC matrix (r_l x r_r)
%
%[MAC_data] = MacW2(fig,Phi_l,NodesPhi_l,W2_l,Phi_r,NodesPhi_r,W2_r,k)
%
%   New algorithm using eigenvalues error
%   MAC = MAC.*exp(-k*(abs(2*(W2_l-W2_r)/(W2_l+W2_r)))^0.5)
%   (old form MAC = MAC.*abs((W2_l-W2_r)/(W2_l+W2_r)) )

%% Selection of corresponding dofs
% W2_l=abs(W2_l).*sign(real(W2_l));
% W2_r=abs(W2_r).*sign(real(W2_r));
% W2_l=imag(W2_l);
% W2_r=imag(W2_r);
%Initialisation of the vector row. It will contain, for a generic i
%position, the row of the second set of eigenvectors that contains the i-th
%DOF of the first set
row = zeros(size(Phi_l,1),1);
for count = 1 : size(Phi_l,1)
    %Since only a single value is ammitted using this form, there will be
    %an error if more than a single value is returned by find. This, 
    %however, would mean that there are duplicated nodes or DOF, that is 
    %impossible (or an error)
    row(count) = find((NodesPhi_r(:,1)+0.1*NodesPhi_r(:,2))==(NodesPhi_l(count,1)+0.1*NodesPhi_l(count,2)));
end
%Creation of the identity matrix, with the dimension of the second set of
%eigenvectors (that is usualy the largest one)
I = speye(size(Phi_r,1));
%Permutation matrix, obtained from the identical one considering only the
%DOFs, in the second set, that are present in the first one
P = I(row,:);
%The second set of eigenvectors is updated, considering only the nodes that
%are present in the first one
Phi_r_red = Phi_r;
Phi_r = P*Phi_r;

%% MAC evaluation
%MAC_data is initialised. It will contain the MAC index value for each 
%possible combination of eigenvectors, between the first and the second set
MAC_data = zeros(size(W2_l,1),size(W2_r,1));
%ErrW2 is initialised. It will contain the error about the eigenvalues.
%Once the eigenvectors are coupled, the correspondent eigenvalues are
%compared, calculating the difference between them
ErrW2 = zeros(size(W2_l,1),size(W2_r,1));
for count1 = 1 : size(W2_l,1)
    for count2 = 1 : size(W2_r,1)
        if norm(Phi_l(:,count1)) > 0 && norm(Phi_r(:,count2)) > 0
            %ErrW2 is equal to 100 when the two eigenvalues are identical,
            %while is equal to 0 when their difference is close to Inf
            %ErrW2(count1,count2) = 100*(1-abs((W2_r(count2,1)-W2_l(count1,1))/(W2_r(count2,1)+W2_l(count1,1))));
            ErrW2(count1,count2) = 100*exp(-k*(abs(2*(W2_r(count2,1)-W2_l(count1,1)) ...
                /(W2_r(count2,1)+W2_l(count1,1))))^0.5);
            %MAC_data is equal to 100 when the two eigenvectors are
            %identical and 0 when they are totally different
            %MAC_data(count1,count2) = 100*((abs(Phi_l(:,count1).'*Phi_r(:,count2)))^2) ...
            %    /((Phi_l(:,count1).'*Phi_l(:,count1))*(Phi_r(:,count2).'*Phi_r(:,count2)));
%             MAC_data(count1,count2) = 100*((abs(Phi_l(:,count1)'*Phi_r(:,count2))+abs(Phi_l(:,count1).'*Phi_r(:,count2)))^2 ...
%                 /((Phi_l(:,count1)'*Phi_l(:,count1)+abs(Phi_l(:,count1).'*Phi_l(:,count1)))*(Phi_r(:,count2)'*Phi_r(:,count2)+abs(Phi_r(:,count2).'*Phi_r(:,count2)))));
       
           MAC_data(count1,count2) = 100*((abs(Phi_l(:,count1)'*Phi_r(:,count2))+abs(Phi_l(:,count1).'*Phi_r(:,count2)))^2 ...
                /((Phi_l(:,count1)'*Phi_l(:,count1)+abs(Phi_l(:,count1).'*Phi_l(:,count1)))*(Phi_r(:,count2)'*Phi_r(:,count2)+abs(Phi_r(:,count2).'*Phi_r(:,count2)))));
        
        
        end
    end
end
%MAC_data = ErrW2;
temp = MAC_data;
%MAC_data is defined using both the parameters: MAC_data, related to the
%eigenvectors and ErrW2, related to the eigenvalues. Since both have values
%between 0 and 100, it is necessary to multiply by 0.01, so that the new
%MAC_data is a number <100
MAC_data = 0.01*MAC_data.*ErrW2;

%% Plotting MAC
if fig>0
    figure(fig)
    set(gcf,'Color',[1 1 1]);
    set(gcf,'Name','MAC');
    [count1,count2] = meshgrid(1:size(W2_l,1),1:size(W2_r,1));
    scatterbar3(count1,count2,MAC_data.',0.8)
    clear count1 count2
    colormap(flipud(hot))
    colorbar
    view(60,70)
    %title(['\bfModal Assurance Criterion'],'FontName','TimesNewRoman','FontSize',16)
    xlabel('FEA_{left}','FontSize',12)
    ylabel('FEA_{right}','FontSize',12)
    zlim([0 100]);
    temp = daspect;
    daspect([min(temp(1,1:2)) min(temp(1,1:2)) temp(1,3)]);
    grid on
end

%% W2 correlation (TBI)
%Dimensions of the input data
[n_Phi_ex m_Phi_ex] = size(Phi_l);
[n_NodesPhi_ex m_NodesPhi_ex] = size(NodesPhi_l);
[n_W2_l m_W2_l] = size(W2_l);
[n_Phi_an m_Phi_an] = size(Phi_r);
[n_NodesPhi_an m_NodesPhi_an] = size(NodesPhi_r);
[n_W2_r m_W2_r] = size(W2_r);


%Initialisation of the vector row. It will contain, for a generic i
%position, the row of the second set of eigenvectors that contains the i-th
%DOF of the first set
row = zeros(n_Phi_ex,1);
for count = 1 : n_Phi_ex
    %Since only a single value is ammitted using this form, there will be
    %an error if more than a single value is returned by find. This, 
    %however, would mean that there are duplicated nodes or DOF, that is 
    %impossible (or an error)
    row(count) = find((NodesPhi_r(:,1)+0.1*NodesPhi_r(:,2))==(NodesPhi_l(count,1)+0.1*NodesPhi_l(count,2)));
end
%Creation of the identity matrix, with the dimension of the second set of
%eigenvectors (that is usualy the largest one)
I = speye(n_Phi_an);
%Permutation matrix, obtained from the identical one considering only the
%DOFs, in the second set, that are present in the first one
P = I(row,:);
%The second set of eigenvectors is updated, considering only the nodes that
%are present in the first one
Phi_an_red = Phi_r;
Phi_an = P*Phi_r;

%Creation of a vector, from 1 to the number of considered eigenvectors
EMA_FEA = (1:n_W2_l).';
%Permutation matrix
P = zeros(n_W2_r);
%MAC realized by choosing, for each eigenvectors of the first set, the best
%match of the second set
for count = 1 : n_W2_l
    %pos contains the index to be used to obtain, from MAC_data, a vector
    %containing its sorted values. Using pos, its last element directs to
    %the maximum value of the MAC index for the considered eigenvector
    [~, pos] = sort(MAC_data(count,:));
    %count1 will be a progressive index, while count2 will be the trigger
    %used to exit the while cycle when the optimum condition is reached
    count1 = 0;
    count2 = 0;
    %The following cycle continues until a "free spot" is reached. In
    %particular, the cycle is repeated until the eigenvector from the first
    %set find a free eigenvector of the second set. Since pos is used, that
    %one will be the one with the maximum MAC index, between the ones not
    %already "taken"
    while count2==0
        if max(P(:,pos(n_W2_r-count1)))==0
            if (n_W2_r-count1-1)~=0
                if MAC_data(count,pos(n_W2_r-count1))==MAC_data(count,pos(n_W2_r-count1-1))
                    if max(P(:,pos(n_W2_r-count1)))==0 & max(P(:,pos(n_W2_r-count1-1)))==0
                        count1 = count1+1;
                    else
                        EMA_FEA(count,2) = pos(n_W2_r-count1);
                        P(count,pos(n_W2_r-count1)) = 1;
                        count2 = 1;
                    end
                else
                    EMA_FEA(count,2) = pos(n_W2_r-count1);
                    P(count,pos(n_W2_r-count1)) = 1;
                    count2 = 1;
                end
            else
                EMA_FEA(count,2) = pos(n_W2_r-count1);
                P(count,pos(n_W2_r-count1)) = 1;
                count2 = 1;
            end
        else
            count1 = count1+1;
        end
    end
end
%If the number of eigenvectors of the second set is larger of the ones of
%the first set, the not assigned eigenvectors are listed one after the
%other. P matrix is modified so that these last eigenvectors are not lost
if n_W2_r > n_W2_l
    count1 = n_W2_l+1;
    for count3 = 1 : n_W2_r
        if max(P(1:(count1-1),count3))==0
            P(count1,count3) = 1;
            count1 = count1+1;
        end
    end
end
%Eigenvectors and eigenvalues are updated, so that they are in the desired
%order (to have the highest MAC index, with preference to the first
%eigenvectors)
Phi_an_red = Phi_an_red*P';
W2_r = P*W2_r;

function scatterbar3(X,Y,Z,width)
%SCATTERBAR3   3-D scatter bar graph.
[r,c]=size(Z);
for j=1:r,
    for k=1:c,
        if ~isnan(Z(j,k))
            drawbar(X(j,k),Y(j,k),Z(j,k),width/2)
        end
    end
end

zlim=[min(Z(:)) max(Z(:))];
if zlim(1)>0,zlim(1)=0;end
if zlim(2)<0,zlim(2)=0;end
axis([min(X(:))-width max(X(:))+width min(Y(:))-width max(Y(:))+width zlim])
caxis([min(Z(:)) max(Z(:))])

function drawbar(x,y,z,width)
h(1)=patch([-width -width width width]+x,[-width width width -width]+y,[0 0 0 0],'b');
h(2)=patch(width.*[-1 -1 1 1]+x,width.*[-1 -1 -1 -1]+y,z.*[0 1 1 0],'b');
h(3)=patch(width.*[-1 -1 -1 -1]+x,width.*[-1 -1 1 1]+y,z.*[0 1 1 0],'b');
h(4)=patch([-width -width width width]+x,[-width width width -width]+y,[z z z z],'b');
h(5)=patch(width.*[-1 -1 1 1]+x,width.*[1 1 1 1]+y,z.*[0 1 1 0],'b');
h(6)=patch(width.*[1 1 1 1]+x,width.*[-1 -1 1 1]+y,z.*[0 1 1 0],'b');
set(h,'facecolor','flat','FaceVertexCData',z)
