%% Segunda Prova

clear all; close all; clc

%% Dados de entrada

tic;

L1=200E-3;
L2=300E-3;
L3=700E-3;
L4=200E-3;
L5=100E-3;
L=[L1 L2 L3 L4 L5];
LK=1100E-3;
D1= 80E-3;
D2=110E-3;
D3=150E-3;
D4=100E-3;
D5= 80E-3;
DD=300E-3;
I1=pi*D1^4/64;
I2=pi*D2^4/64;
I3=pi*D3^4/64;
I4=pi*D4^4/64;
I5=pi*D5^4/64;
A1=pi*(D1/2)^2;
A2=pi*(D2/2)^2;
A3=pi*(D3/2)^2;
A4=pi*(D4/2)^2;
A5=pi*(D5/2)^2;
H=50E-3;
E=201e9;
dens=7860;
M=pi*(DD/2)^2*H*dens;
K0=2E6;
K1=10E20;
SEC=5;
NMODOS=10;
N=3000;
X=linspace(0,sum(L),N+1);

%BETA=[4.730041 7.853205 10.995608 14.137165];
for i=1:NMODOS
    BETA(i)=3.1364*i+1.5881;
end

%% Modos de vibrar do modelo contínuo

C=1;
MODOS(:,1)=C*ones(1,length(X));
MODOS(:,2)=linspace(-C,C,length(X));

for j=3:2+length(BETA)
    ALPHA=(sin(BETA(j-2))-sinh(BETA(j-2)))/(cosh(BETA(j-2))-cos(BETA(j-2)));
    AUX1=BETA(j-2).*X/(max(X));
    AUX2=C*(sin(AUX1)+sinh(AUX1)+ALPHA*(cos(AUX1)+cosh(AUX1)));
    MODOS(:,j)=AUX2;
end

%% Normalização dos modos

for j=1:2+length(BETA)
    soma=0;
    for i=1:length(MODOS(:,j))
        soma=soma+MODOS(i,j)^2;
    end
    MODOS(:,j)=MODOS(:,j)./sqrt(soma);
end

%% Separaçao para cada secção

aux1=0;
aux2=0;
aux3=0;
aux4=0;
aux5=0;

for i=1:length(X)
    if X(i)<=L1;
        aux1=aux1+1;
        X1(aux1)=X(i);
        MSECAO1(aux1,:)=MODOS(i,:);
    end
    if L1<X(i)&& X(i)<=L1+L2
        aux2=aux2+1;
        X2(aux2)=X(i);
        MSECAO2(aux2,:)=MODOS(i,:);
    end
    if L1+L2<X(i)&& X(i)<=L1+L2+L3
        aux3=aux3+1;
        X3(aux3)=X(i);
        MSECAO3(aux3,:)=MODOS(i,:);
    end
    if L1+L2+L3<X(i)&& X(i)<=L1+L2+L3+L4
        aux4=aux4+1;
        X4(aux4)=X(i);
        MSECAO4(aux4,:)=MODOS(i,:);
    end
    if L1+L2+L3+L4<X(i)
        aux5=aux5+1;
        X5(aux5)=X(i);
        MSECAO5(aux5,:)=MODOS(i,:);
    end
end

% Montagem das derivadas

% Derivada primeira

for i=1:length(X)-1
    DER(i,:)=(MODOS(i+1,:)-MODOS(i,:))/(X(i+1)-X(i));
end
for i=1:length(X1)-1
    DER1(i,:)=(MSECAO1(i+1,:)-MSECAO1(i,:))/(X1(i+1)-X1(i));
end
for i=1:length(X2)-1
    DER2(i,:)=(MSECAO2(i+1,:)-MSECAO2(i,:))/(X2(i+1)-X2(i));
end
for i=1:length(X3)-1
    DER3(i,:)=(MSECAO3(i+1,:)-MSECAO3(i,:))/(X3(i+1)-X3(i));
end
for i=1:length(X4)-1
    DER4(i,:)=(MSECAO4(i+1,:)-MSECAO4(i,:))/(X4(i+1)-X4(i));
end
for i=1:length(X5)-1
    DER5(i,:)=(MSECAO5(i+1,:)-MSECAO5(i,:))/(X5(i+1)-X5(i));
end

% Derivada segunda

for i=1:length(X)-2
    DERDER(i,:)=(DER(i+1,:)-DER(i,:))/(X(i+1)-X(i));
end
for i=1:length(X1)-2
    DERDER1(i,:)=(DER1(i+1,:)-DER1(i,:))/(X1(i+1)-X1(i));
end
for i=1:length(X2)-2
    DERDER2(i,:)=(DER2(i+1,:)-DER2(i,:))/(X2(i+1)-X2(i));
end
for i=1:length(X3)-2
    DERDER3(i,:)=(DER3(i+1,:)-DER3(i,:))/(X3(i+1)-X3(i));
end
for i=1:length(X4)-2
    DERDER4(i,:)=(DER4(i+1,:)-DER4(i,:))/(X4(i+1)-X4(i));
end
for i=1:length(X5)-2
    DERDER5(i,:)=(DER5(i+1,:)-DER5(i,:))/(X5(i+1)-X5(i));
end

% Montagem de matriz de integrais

for i=1:2+length(BETA)
    for j=1:2+length(BETA)
        soma=0;
        for k=1:length(MSECAO1(:,1))-1
            soma=soma+(X1(k+1)-X1(k))*(MSECAO1(k+1,i)*MSECAO1(k+1,j)+MSECAO1(k,i)*MSECAO1(k,j))/2;
        end
        INT1MODO(i,j)=soma;
        soma=0;
        for k=1:length(MSECAO2(:,1))-1
            soma=soma+(X2(k+1)-X2(k))*(MSECAO2(k+1,i)*MSECAO2(k+1,j)+MSECAO2(k,i)*MSECAO2(k,j))/2;
        end
        INT2MODO(i,j)=soma;
        soma=0;
        for k=1:length(MSECAO3(:,1))-1
            soma=soma+(X3(k+1)-X3(k))*(MSECAO3(k+1,i)*MSECAO3(k+1,j)+MSECAO3(k,i)*MSECAO3(k,j))/2;
        end
        INT3MODO(i,j)=soma;
        soma=0;
        for k=1:length(MSECAO4(:,1))-1
            soma=soma+(X4(k+1)-X4(k))*(MSECAO4(k+1,i)*MSECAO4(k+1,j)+MSECAO4(k,i)*MSECAO4(k,j))/2;
        end
        INT4MODO(i,j)=soma;
        soma=0;
        for k=1:length(MSECAO5(:,1))-1
            soma=soma+(X5(k+1)-X5(k))*(MSECAO5(k+1,i)*MSECAO5(k+1,j)+MSECAO5(k,i)*MSECAO5(k,j))/2;
        end
        INT5MODO(i,j)=soma;
        soma=0;
        for k=1:length(DERDER1(:,1))-1
            soma=soma+(X1(k+1)-X1(k))*(DERDER1(k+1,i).*DERDER1(k+1,j)+DER1(k,i).*DERDER1(k,j))/2;
        end
        INT1DER(i,j)=soma;
        soma=0;
        for k=1:length(DERDER2(:,1))-1
            soma=soma+(X2(k+1)-X2(k))*(DERDER2(k+1,i).*DERDER2(k+1,j)+DERDER2(k,i).*DERDER2(k,j))/2;
        end
        INT2DER(i,j)=soma;
        soma=0;
        for k=1:length(DERDER3(:,1))-1
            soma=soma+(X3(k+1)-X3(k))*(DERDER3(k+1,i).*DERDER3(k+1,j)+DERDER3(k,i).*DERDER3(k,j))/2;
        end
        INT3DER(i,j)=soma;
        soma=0;
        for k=1:length(DERDER4(:,1))-1
            soma=soma+(X4(k+1)-X4(k))*(DERDER4(k+1,i).*DERDER4(k+1,j)+DERDER4(k,i).*DERDER4(k,j))/2;
        end
        INT4DER(i,j)=soma;
        soma=0;
        for k=1:length(DERDER5(:,1))-1
            soma=soma+(X5(k+1)-X5(k))*(DERDER5(k+1,i).*DERDER5(k+1,j)+DERDER5(k,i).*DERDER5(k,j))/2;
        end
        INT5DER(i,j)=soma;
    end
end

% Coeficientes da matriz

% Posição das Molas

MODOK11 = MODOS(1,:);

xk=0;
aux=0;

while xk<LK
    aux=aux+1;
    xk=X(aux);
end

MODOK0 = MODOS(aux,:);

xk=0;
aux=0;

while xk<sum(L)
    aux=aux+1;
    xk=X(aux);
end

MODOK12 = MODOS(aux,:);

% Posição da Massa

xm=0;
aux=0;

while xm<L1+L2
    aux=aux+1;
    xm=X(aux);
end

MODOM = MODOS(aux,:);

% Massa e Rigidez

for i=1:2+length(BETA)
    for j=1:2+length(BETA)
        mij(i,j)=dens*(A1*INT1MODO(i,j)+A2*INT2MODO(i,j)+A3*INT3MODO(i,j)+A4*INT4MODO(i,j)+A5*INT5MODO(i,j))+M*MODOM(i)*MODOM(j);
        kij(i,j)=E*(I1*INT1DER(i,j)+I2*INT2DER(i,j)+I3*INT3DER(i,j)+I4*INT4DER(i,j)+I5*INT5DER(i,j))+K1*(MODOK11(i)*MODOK11(j)+MODOK12(i)*MODOK12(j))+K0*MODOK0(i)*MODOK0(j);
    end
end

%% Obtendo as frequências naturais

[MODOS_A,FRAD]=eig(kij,mij);
[u,v]=sort(diag(FRAD));

MASSA=MODOS_A'*mij*MODOS_A;

for i=1:length(u)
    MODOS_A(:,i)=MODOS_A(:,i)/sqrt(MASSA(i,i));
end

% for j=1:2+length(BETA)
%     soma=0;
%     for i=1:length(u)
%         soma=soma+MODOS_A(i,j)^2;
%     end
%     MODOS_A(:,j)=MODOS_A(:,j)./sqrt(soma);
% end

for i=1:length(u)
    MODOS_A_O(:,i)=MODOS_A(:,v(i));
end

FHZ=sqrt(u)/(2*pi);

F1=FHZ(1)
F2=FHZ(2)
F3=FHZ(3)
F4=FHZ(4)
F5=FHZ(5)

MODOS_REAIS=MODOS*MODOS_A_O;

figure(1)
subplot(1,2,1)
plot(X,MODOS(:,1),'b',X,MODOS(:,2),'r',X,MODOS(:,3),'k',X,MODOS(:,4),'m',X,MODOS(:,5),'g')
title('Modos do modelo continuo')
grid
legend('1 Modo','2 Modo','3 Modo','4 Modo','5 Modo')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS(:,1:5))) 1.1*max(max(MODOS(:,1:5)))])

subplot(1,2,2)
plot(X,MODOS_REAIS(:,1),'b',X,MODOS_REAIS(:,2),'r',X,MODOS_REAIS(:,3),'k',X,MODOS_REAIS(:,4),'m',X,MODOS_REAIS(:,5),'g')
title('Modos reais do modelo')
grid
legend('1 Modo','2 Modo','3 Modo','4 Modo','5 Modo')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])

%% Leitura dos valores obtidos no ANSYS

load PRNSOL1MODOY.txt
load PRNSOL2MODOY.txt
load PRNSOL3MODOY.txt
load PRNSOL4MODOY.txt
load PRNSOL5MODOY.txt
load PRNSOL6MODOY.txt
load NLIST.txt

xnosaux=NLIST(1:73,2);
[a,b]=sort(xnosaux);

for i=1:length(a)
    xnos(i)=NLIST(b(i),2);
    modo1(i)=PRNSOL1MODOY(b(i),2);
    modo2(i)=PRNSOL2MODOY(b(i),2);
    modo3(i)=PRNSOL3MODOY(b(i),2);
    modo4(i)=PRNSOL4MODOY(b(i),2);
    modo5(i)=PRNSOL5MODOY(b(i),2);
    modo6(i)=PRNSOL6MODOY(b(i),2);
end
    
figure(2)
plot(X,MODOS_REAIS(:,1),'b',xnos,modo1,'r')
title('Primeiro modo real')
legend('Modos Assumidos','ANSYS APDL')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])
grid

figure(3)
plot(X,MODOS_REAIS(:,2),'b',xnos,-modo2,'r')
title('Segundo modo real')
legend('Modos Assumidos','ANSYS APDL')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])
grid

figure(4)
plot(X,MODOS_REAIS(:,3),'b',xnos,modo3,'r')
title('Terceiro modo real')
legend('Modos Assumidos','ANSYS APDL')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])
grid

figure(5)
plot(X,MODOS_REAIS(:,4),'b',xnos,-modo5,'r')
title('Quarto modo real')
legend('Modos Assumidos','ANSYS APDL')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])
grid

figure(6)
plot(X,MODOS_REAIS(:,5),'b',xnos,modo6,'r')
title('Quinto modo real')
legend('Modos Assumidos','ANSYS APDL')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])
grid

figure(7)
plot(X,MODOS_REAIS(:,6),'b',xnos,modo6,'r')
title('Sexto modo real')
legend('Modos Assumidos','ANSYS APDL')
xlabel('Posição')
ylabel('Amplitude')
axis([0 sum(L) 1.1*min(min(MODOS_REAIS(:,1:5))) 1.1*max(max(MODOS_REAIS(:,1:5)))])
grid
tempo=toc
%pause
%% Receptância

load PRVAR1.txt
load PRVAR2.txt

F1=PRVAR1(:,1);
F2=PRVAR2(:,1);
H1=PRVAR1(:,2);
H2=PRVAR2(:,2);

% Montagem do vetor de amortecimento

KSI(1:NMODOS+2)=0.002;

% Localização da excitação

xe=0;
auxe=0;

while xe<L1+L2
    auxe=auxe+1;
    xe=X(auxe);
end

% Localização da resposta

xr1=0;
auxr1=0;

while xr1<L1
    auxr1=auxr1+1;
    xr1=X(auxr1);
end

xr2=0;
auxr2=0;

while xr2<L1+L2
    auxr2=auxr2+1;
    xr2=X(auxr2);
end

% Primeiro par excitação e resposta

omega=0:(10*2*pi):3000*2*pi;     %Discretização no domínio da frequência

FORCA=zeros(N+1,1);
FORCA(auxe,1)=1;

for w=1:length(omega)  
    for j=1:NMODOS+2
        D(j,j)=1./(-omega(w)^2+1i*2*KSI(j).*omega(w)*sqrt(u(j))+u(j));
    end
    alpha(:,w)=MODOS_REAIS*D*MODOS_REAIS'*FORCA;
end

figure(8)
plot(omega./(2*pi),20*log10(abs(alpha(auxr1,:))),'b',F1,20*log10(abs(H1)),'r')
legend('Modo Assumido','ANSYS APDL')
title('Comparação entre as receptâncias obtidas - Primeiro ponto')
xlabel('Frequência')
ylabel('Amplitude')
axis([0 F5*1.01 1.1*min(20*log10(abs(H2))) 1.1*max(20*log10(abs(H2)))])
grid

figure(9)
plot(omega./(2*pi),20*log10(abs(alpha(auxr2,:))),'b',F2,20*log10(abs(H2)),'r')
legend('Modo Assumido','ANSYS APDL')
title('Comparação entre as receptâncias obtidas - Segundo ponto')
xlabel('Frequência')
ylabel('Amplitude')
axis([0 F5*1.01 1.1*min(20*log10(abs(H2))) 1.1*max(20*log10(abs(H2)))])
grid