clear all
close all
clc

%% Script computing the MAC matrix from TestLab .txt file
% hoainam.huynh@umons.ac.be (December 2018 (UVic))

%%%%%%%%%%%%%%
% User input %
%%%%%%%%%%%%%%

% Data
TestLabFileName = 'TestLab_Mode_'; % name template Mode_x_n : ex: Mode_x_1
EasyDynFileName = 'EasyDyn_3D_LMS_Mode_';
nbrmode = 4;

%%%%%%%%%%%
% Program %
%%%%%%%%%%%

% TestLab mode shapes
% Editing files to read them 
for n=1:nbrmode
        % Replace comma by point in numbers
        FileName2 = strcat(TestLabFileName,num2str(n));
        wholefile = fileread(strcat(FileName2,'.txt'));  % read in entire file
        newfiledata1 = strrep(wholefile,',','.');        % replace , by .
        fileid = fopen(strcat(FileName2, '.txt'),'w');   % open file to write
        fprintf(fileid,'%s',newfiledata1);               % print to file
        fclose(fileid);
end
% Load modes
datamode1 = load(strcat(TestLabFileName,'1.txt')); % load mode 1
[n_row, n_col] = size(datamode1);
TestLabModeShape = zeros(n_row,nbrmode);
TestLabModeShape(:,1) = complex(datamode1(:,1),datamode1(:,2));
for n=2:nbrmode
      FileName2 = strcat(TestLabFileName,num2str(n));
      datamode2 = load(strcat(FileName2,'.txt')); % load mode n
      TestLabModeShape(:,n) = complex(datamode2(:,1),datamode2(:,2));
end

% EasyDyn modes shapes (but computed in TestLab
%, in EasyDyn time domain simulation of vibration)
% Editing files to read them 
for n=1:nbrmode
        % Replace comma by point in numbers
        FileName2 = strcat(EasyDynFileName,num2str(n));
        wholefile = fileread(strcat(FileName2,'.txt'));  % read in entire file
        newfiledata1 = strrep(wholefile,',','.');        % replace , by .
        fileid = fopen(strcat(FileName2, '.txt'),'w');   % open file to write
        fprintf(fileid,'%s',newfiledata1);               % print to file
        fclose(fileid);
end
% Load modes
datamode1 = load(strcat(EasyDynFileName,'1.txt')); % load mode 1
[n_row, n_col] = size(datamode1);
EasyDynModeShape = zeros(n_row,nbrmode);
EasyDynModeShape(:,1) = complex(datamode1(:,1),datamode1(:,2));
for n=2:nbrmode
      FileName2 = strcat(EasyDynFileName,num2str(n));
      datamode2 = load(strcat(FileName2,'.txt')); % load mode n
      EasyDynModeShape(:,n) = complex(datamode2(:,1),datamode2(:,2));
end

% MAC Matrix TestLab
MAC_matrix=zeros(nbrmode,nbrmode);
for mode_row=1:nbrmode
    Psi1 = TestLabModeShape(:,mode_row);
    for mode_col=1:nbrmode
        Psi2 = EasyDynModeShape(:,mode_col);
        MAC_matrix(mode_row,mode_col)=abs(sum(Psi1.*conj(Psi2)))^2/((sum(Psi1.*conj(Psi1)))*(sum(Psi2.*conj(Psi2))));
    end
    
end

% Plot of the MAC matrix
[l c] = size(MAC_matrix) ;
MAC_matrix = [MAC_matrix zeros(l,1)] ;
MAC_matrix = [MAC_matrix ; zeros(1,c+1)] ;
[l c] = size(MAC_matrix) ;
figure ;
colormap(cool) ;
pcolor(MAC_matrix) ;
colorbar ;
title('MAC matrix: TestLab vs. EasyDyn') ;
ylabel('LMS TestLab') ;
xlabel('EasyDyn') ;
