%% prog 001 FEM tubo cerrado en ambos extremos sin montaje

clearvars;
close all;
clc;

%este programa calcula los valores propios (frecuencias naturales) y 
%vectores propios (formas de vibracion) de un tubo cerrado en ambos
%extremos de seccion transversal constante A longitud L a una velocidad del
%sonido constante c
%inicialmente el numero de elementos es 4 pero será modificado a un numero
%variable de elementos

%% inicializacion

Ne = 100;      %numero de elementos
N = Ne + 1;     %numero de nodos / grados de libertad / DOF
L = 1.0;        %longitud del tubo en metros (m)
ax = L/(2*Ne);  %mitad de la longitud de un elemento
A = 0.0001;     %area de la sección transversal (m^2)( 1.0e-4) = 1 X 10^(-4) )
c = 340;        %velocidad del sonido
ro = 1.18;      %densidad del aire

%% matriz de masa

%inicializacion
M = zeros(N,N); %matriz incializada con ceros

%en una primera etapa colocaremos de forma directa todos los valores en la
%matriz 
%M = [2  1  0  0  0;
%     1  4  1  0  0;
%     0  1  4  1  0;
%     0  0  1  4  1;
%     0  0  0  1  2];

%utilizaremos la estructura de la matriz para generalizar resultados para
%tubos de seccion cosntante donde todos los elementos tienen el mismo
%tamaño

%proceso incial
M(1  , 1  )   = 2;
M(N-1, N  )   = 1;
M(N  , N-1)   = 1;
M(N  , N  )   = 2;

%iteracion principal desde la 2da fila/columna a la penultima fila/columna
for n1 = 2:N-1
    M(n1  ,n1    ) = 4;	%diagonal llena de 4
    M(n1-1,n1    ) = 1;	%super diagonal llena de 1
    M(n1  ,n1 - 1) = 1; %sub diagonal llena de 1
end;    

%factores de ajuste
M = (A*ax)/(3*c^2)*M; 

%% matriz de rigidez

%inicializacion
K = zeros(N,N);

%en una primera etapa colocaremos de forma directa todos los valores en la
%matriz 

%K = [ 1 -1  0  0  0;
%     -1  2 -1  0  0;
%      0 -1  2 -1  0;
%      0  0 -1  2 -1;
%      0  0  0 -1  1];

%utilizaremos la estructura de la matriz para generalizar resultados para
%tubos de seccion cosntante donde todos los elementos tienen el mismo
%tamaño

%proceso incial
K(1  , 1  )   =  1;
K(N-1, N  )   = -1;
K(N  , N-1)   = -1;
K(N  , N  )   =  1;

%iteracion principal desde la 2da fila/columna a la penultima fila/columna
for n1 = 2:N-1
    K(n1  ,n1    ) =  2;    %diagonal llena de 2
    K(n1-1,n1    ) = -1;    %super diagonal llena de -1
    K(n1  ,n1 - 1) = -1;    %sub diagonal llena de -1
end;    

%factores de ajuste
K = A/(2*ax)*K; 
 
%% problema de valores y vectores propios
%sirve para determinar las frecuencias naturales y la distribucion de la
%presion sonora a dichas frecuencias naturales
%resolvemos la ecuacion
%[lam*M + K]*phi = 0
%M es matriz (N X N)
%K es matriz (N X N) 
%phi es vector (N X 1)
%lam es escalar lam = -omega^2 (1 X 1)
%fin de que la solucion sea no trivial, es decir distinta de cero
%obtendremos N soluciones para lam (es decir N frecuencias naturales) y N
%soluciones para phi(es decir N distribuciones de presion sonora)

%eig soluciona el problema de valores propios
[PHI, LAM] = eig(K,M);

%PHI es la matriz (N X N)formada por los vectores propios 
%PHI = %[phi1, phi2,..., phiN] los phi son vectores columna
%cada phi es (N X 1)

%LAM es la matriz diagonal formada por los valores propios
%LAM = diag(lam1, lam2, ..., lamN) 
%LAM = diag(omega1^2, omega2^2, ..., omegaN^2)

%frecuencias naturales aproximadas por FEM
%diag estrae la diagonal principal de una matriz
fFEM = real( sqrt (diag(LAM) ) )/(2*pi);

%normalizacion de vectores propios
mnorm = diag(transpose(PHI)*M*PHI); 
for n1 = 1:N
    PHI(:,n1) = PHI(:,n1)/sqrt( mnorm(n1) );
end;

%% resultados frecuencias

%resultados teoricos
fTEO = zeros(N,1);
for n1 = 1:N
    fTEO(n1) = (n1 - 1)*c/(2*L);
end;

%error en frecuencias
errof = zeros(N,1);
for n1 = 2:N
    errof(n1) = 100*abs( (fTEO(n1) - fFEM(n1))/fTEO(n1) );
end;

%grafico de frecuencias
figure(1)
plot(fTEO,fFEM,'-o');
title('frecuencias teoricas vs frecuencias fem');
xlabel('fTEO (Hz)');
ylabel('fFEM (Hz)');
grid on;
box on;

%grafico de error
figure(2)
plot(errof,'-o');
title('error');
xlabel('numero de modo');
ylabel('error (%)');
grid on;
box on;

TablaError = [fTEO fFEM errof];

%% resultados distribucion de presion sonora

%cordenadas
dx = L/Ne;  %incremento o longitud de cada elemento
xcorde = (0:N-1)*dx - L/2; %nos da las coordenadas desde -L/2 a L/2 

%grafico
figure(3)

subplot(2,2,1)
plot( xcorde, PHI(:,2)/max( abs( PHI(:,2) ) ) );
title('1er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;

subplot(2,2,2)
plot( xcorde, PHI(:,3)/max( abs( PHI(:,3) ) ) );
title('2er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;

subplot(2,2,3)
plot( xcorde, PHI(:,4)/max( abs( PHI(:,4) ) ) );
title('3er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;    
    
subplot(2,2,4)
plot( xcorde, PHI(:,5)/max( abs( PHI(:,5) ) ) );
title('4er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;