% Exercise 13: Coeficient of Influence
clear all;
clc
%
% Calculate the modal parameters of the Cantilever beam
%
% System data
%
L = 0.6;                % length of the beam
b = 40.e-3;             % width of the beam
h = 4.e-3;              % thickness of the beam
E = 7e10;               % Elasticity modulus of the steel material
I = b*h^3/12;           % Second order moment of area
Aa = b*h;                % Area
ni = 0.33;              % Poisson coeficient
ro = 2700;              % Densidity   


N=10;                   % Number of elements of the beam
k=1;                    % POint of the Response 
s=N;                    % Point of the excitation

NN=5;                   % Number of modes used in the trucate matrix
ksi = 2e-3;             % modal damping
fmax=500;               % maximal frequency


le = L/N;               % Length of the element
me = b*h*le*ro;         % Mass of the elements                 

M = zeros(N);           % Mass Matrix
A = zeros(N);           % Fexibility matrix

for ii=1:N
    M(ii,ii)=me;
    a(ii)=ii*le;
    x(ii)=ii*le;
end
M(N,N)= me/2;


for ii=1:N
    for jj=1:N
        if x(jj) <= a(ii)	% Nós antes da força
        
            A(ii,jj)=1/(E*I)*(-x(jj)^3/6 + a(ii)*x(jj)^2/2);

        else		% Nós depois da força

            A(ii,jj)=1/(E*I)*(1/6*(x(jj)-a(ii))^3 - x(jj)^3/6 + a(ii)*x(jj)^2/2);

        end
    end
end
    
    
% Inversao de A 

K = inv(A);             % Stiffness matrix

% Eigenvalue problems

[F,LL]=eig(K,M);         % F: eigenvector matrix (modal matrix)
                        % L: eigenvalue matrix (spectral matrix                       
                        
% Sorting the eigenvalue e eigenvector                       
[lam,II]=sort(diag(LL));
for ii=1:length(LL)
    FF(:,ii)=F(:,II(ii));
end

% Ortogonality test
F=FF;
Mr = F.'*M*F;
Kr = F.'*K*F;



%% Mr(1:5,1:5);
% Kr(1:5,1:5);

% Ortonomalization the matrix F

for ii=1:N
    F(:,ii) = F(:,ii)/ Mr(ii,ii)^0.5;
end

Mr = F.'*M*F;
Kr = F.'*K*F;


fn=sqrt(lam)/2/pi;          % undamped natural frequency

figure(1)
plot(x,F(:,1),x,F(:,2),x,F(:,3),x,F(:,4))
title('VIGA - SEM MASSA')


for j=1:fmax
    f=j;
    alfa=0.;
    for jj=1:NN
        alfa=alfa+(F(k,jj)*F(s,jj))/(-(2*pi*f)^2+(2*pi*fn(jj))^2+i*2*ksi*(2*pi*fn(jj))*(2*pi*f)); 
    end
    alfaks(j)=alfa;
    freq(j)=f;
end

Inerks=-(2*pi*freq).^2.*alfaks;

figure(6)
plot(freq,20*log10(abs(alfaks)))
title('VIGA - SEM MASSA')

figure(7)
plot(freq,20*log10(abs(Inerks)))
title('VIGA - SEM MASSA')
