clear all;

% Parametros modais de uma viga engastada livre.
%
%           Coeficients de Infuencia
%
%
% Dados
%
L = 0.6;              % Longitude da viga
b = 40.e-3;              % Largura da viga
h = 4.e-3;               % Espessura da viga
E = 7.1e10;              % Modulo de Yung do aço
I = b*h^3/12;           % Momento de segunda ordem de area
A = b*h;                % Area da seçao transversal
ni = 0.33;              % Coeficiente de Poisson
ro = 2770;              % Densidade do aço     


N=20;
k=13;
s=2;

le = L/N;              % Longitude de cada elemento
me = b*h*le*ro;         % Massa de cada elemento                 % numero de coordenadas generalizadas

M = zeros(N);
A = zeros(N);

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

for i=1:N
    for j=1:N
        if x(j) <= a(i)
            A(i,j)=1/(E*I)*(-x(j)^3/6 + a(i)*x(j)^2/2);
        else
            A(i,j)=1/(E*I)*(1/6*(x(j)-a(i))^3 - x(j)^3/6 + a(i)*x(j)^2/2);
        end
    end
end
    
% Inversao de A 

K = inv(A);

% Problema de autovalores

[F,L]=eig(K,M);         % F matriz de autovetores nao ortonormalizados (matriz modal)
                        % L matriz espectral

                        
[lam,I]=sort(diag(L));
for ii=1:length(L)
    FF(:,ii)=F(:,I(ii));
end

% Teste de ortogonalidade
F=FF;
Mr = F.'*M*F;
Kr = F.'*K*F;

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

% Ortonomalizacao de F atraves da matriz de massa

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

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

fn=sqrt(lam)/2/pi;

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

ksi = 2e-3;

NN=4;

for j=1:1000
    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(4)
plot(freq,20*log10(abs(alfaks)))
title('VIGA - SEM MASSA')

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