% Computes the effective masses for a given COMSOL eigenfrequency analysis.
% 
% Before running this script, you need to run an eigenfrequency analysis 
% in COMSOL, save the COMSOL file and connect your COMSOL program to  
% MATLAB. You don't need to export your fem solution to MATLAB.
% 
% Run the script as any MATLAB script. The script opens a dialog box for 
% selecting the COMSOL file containing the eigenfrequency solution computed 
% by COMSOL. 
% 
% Effective masses are calculated and saved to a Microsoft Excel file. The 
% script will ask you were to save the file and opens it automatically when
% finished.
% 
% Copyright (c) 2010 Lars Tiedemann

clear all;

% The Excel *.xls result file
xls_file_temp = [tempname '.xls'];
xls_save = 'Save file...';
xls_label_frequency = {'f'};
xls_label_sum = {'Total'};
xls_sheet_prefix = '';
xls_sheet_M_total = 'Total Masses';
xls_sheet_M_eff = 'Effective Masses';
xls_sheet_M_frac = 'Mass Fraction';
xls_sheet_M_cum = 'Cumulative Mass Fraction';

% Get the COMSOL *.mph file
[mph_file, mph_path, mph_filter] = uigetfile('*.mph');
flload([mph_path mph_file]);

% Get the eigenvector matrix (nxm)
u = fem.sol.u;
u = real(u(:,:));
% Number of DOF
n = size(u, 1);
% Number of eigenmodes
m = size(u, 2);
% Number of independent variables, e.g. u,v,w 
dims = [];
for i = 1: size(fem.appl, 2)
    dims = [dims fem.appl{1,i}.dim(1,:)];
end
dims = unique(dims);
k = size(dims, 2);

% Get eigenfrequencies (jx1)
eigenfreq = fem.sol.lambda;
eigenfreq(:,:) = abs(imag(eigenfreq(:,:)))/(2*pi);
eigenfreq = eigenfreq';

% Get the mass matrix (nxn)
M = assemble(fem,'out',{'E'},'eigname','lambda');
%Stiffness matrix
%K=assemble(fem,'out',{'K'},'eigname','lambda');

% Rigid body influence vector (nxk)
r = repmat(eye(k), [n/k 1]);

% Total system masses (1xk)
M_total = sum(M*r, 1);

% Modal mass matrix (mxm)
Mr = u'*M*u;

% Compute participation vector
L = u' * M * r;

% Replace the eigenvector by the mass normalized eigenvector (nxm)
%u_n = u*Mr^(-1/2);
% Check correct normalization
%identity_matrix = u_n'*M*u_n;
%fem.sol = femsol(u_n, 'lambda', fem.sol.lambda);
% Save file with normalized eigenvectors
%flsave [mph_path mph_file] fem

% Modal participation factor
Lambda = zeros([m k]);
% Effective mass matrix
M_eff = zeros([m k]);
% Mass fraction
M_frac = zeros([m k]); 
% Cumulative mass fraction
M_cum = zeros([m k]); 

% loop dimensions
for j = 1:k
    % loop modes
    for i = 1:m
        Lambda(i, j) = L(i, j)   /Mr(i, i);
        M_eff(i, j) = L(i, j)^2 /Mr(i, i);
        M_cum(i, j) = sum(M_eff(:, j));
    end
end
% Sum of effective masses
M_sum = sum(M_eff, 1);
for j = 1: k
    M_cum(:, j) = M_cum(:, j) / M_sum(1, j);
    M_frac(:, j) = M_eff(:, j) / M_sum(1, j);
end

% Write output to Excel file
xls_sheet = [xls_sheet_prefix xls_sheet_M_total];
xlswrite(xls_file_temp, dims, xls_sheet, 'A1');
xlswrite(xls_file_temp, M_total, xls_sheet, 'A2');

xls_sheet = [xls_sheet_prefix xls_sheet_M_eff];
xlswrite(xls_file_temp, xls_label_frequency, xls_sheet, 'A1');
xlswrite(xls_file_temp, dims, xls_sheet, 'B1');
xlswrite(xls_file_temp, eigenfreq, xls_sheet, 'A2');
xlswrite(xls_file_temp, M_eff, xls_sheet, 'B2');
xlswrite(xls_file_temp, xls_label_sum, xls_sheet, ['A' int2str(m+2)]);
xlswrite(xls_file_temp, M_sum, xls_sheet, ['B' int2str(m+2)]);

xls_sheet = [xls_sheet_prefix xls_sheet_M_frac];
xlswrite(xls_file_temp, xls_label_frequency, xls_sheet, 'A1');
xlswrite(xls_file_temp, dims, xls_sheet, 'B1');
xlswrite(xls_file_temp, eigenfreq, xls_sheet, 'A2');
xlswrite(xls_file_temp, M_frac, xls_sheet, 'B2');

xls_sheet = [xls_sheet_prefix xls_sheet_M_cum];
xlswrite(xls_file_temp, xls_label_frequency, xls_sheet, 'A1');
xlswrite(xls_file_temp, dims, xls_sheet, 'B1');
xlswrite(xls_file_temp, eigenfreq, xls_sheet, 'A2');
xlswrite(xls_file_temp, M_cum, xls_sheet, 'B2');

% Save excel file and open
xls_file = strrep([mph_path mph_file], '.mph', '.xls');
[xls_file, xls_path, xls_filter] = uiputfile('*.xls', xls_save, xls_file);
movefile(xls_file_temp, [xls_path xls_file]);
winopen([xls_path xls_file]);