%% Prog002FEMTuboCerradoCerradoSinMontajeGeneralizado

clearvars;
close all;
clc;

%este es un programa de elementos finitos que representa un tubo cerrado
%en ambos extremos calcularemos las frecuencias naturales (valores propios)
%y la distribucion de la presion sonora (MNV) dentro del tubo a dichas
%frecuencias (vectores propios)
%para este caso comenzaremos con 100 elementos 101 nodos

%% datos inciales
L = 1.0;        %longitud del tubo
Ne = 100;       %numero de elementos
N = Ne + 1;     %numero de nodos grados de libertad DOF
ax = L/(2*Ne);  %semilongitud de cada elemento
a = 0.01;       %lado/radio del tubo
A = pi*a^2;     %area de la seccion transversal del tubo A
c = 344;        %velocidad del sonido

%recordemos que este es un modelo de ondas planas
%esto se cumple cuando la longitud de onda es mucho mayor que las
%dimensiones que forman el area de la seccion transversal

%% matriz de masa

% 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];

%inicializamos
M = zeros(N,N);
M(  1,1) = 2;
M(N-1,N) = 1;
M(N,N-1) = 1;
M(  N,N) = 2;

for n = 2:N-1
    M(  n,  n) = 4;         %diagonal principal
    M(n-1,  n) = 1;         %diagonal superior
    M(  n,n-1) = 1;         %diagonal inferior
end;

 
M = (A*ax)/(3*c^2)*M;
 
 %% matriz de rigidez
 
% 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];

%inicializamos
K = zeros(N,N);
K(  1,1) =  1;
K(N-1,N) = -1;
K(N,N-1) = -1;
K(  N,N) =  1;

for n = 2:N-1
    K(  n,  n) =  2;        %diagonal principal
    K(n-1,  n) = -1;        %diagonal superior
    K(  n,n-1) = -1;        %diagonal inferior
end;

   
K = A/(2*ax)*K; 
 
%% valores y vectores propios
 
%la función eig() determina los valores y vectores propios de la matriz
[PHI,LAM] = eig(K,M);
 
%frecuencias naturales elementos finitos FEM
f_FEM = diag( real(sqrt(LAM)))/(2*pi);
 
%frecuencias naturales teóricas TEO
f_TEO = zeros(N,1);
for n = 1:N
     f_TEO(n) = (n-1)*c/(2*L);
end;
 
%error
err = 100*abs(f_FEM - f_TEO)./f_TEO;
 
%tabla
f_FEM_f_TEO_err = [f_FEM, f_TEO, err];
 
%% figuras frecuencias valores propios

figure(1)
plot(f_TEO,f_FEM,'-o');
title('frecuencias f_{TEO} vs f_{FEM}');
xlabel('f_TEO (Hz)');
ylabel('f_TEO (Hz)');
grid on;
box on;

figure(2)
plot(err,'-o');
title('porcentaje de error frecuencias f_{TEO} vs f_{FEM}');
xlabel('n');
ylabel('error (%)');
grid on;
box on;
 
%% figuras distribución de presion sonora modal

%coordenadas nodales
x = zeros(N,1);
for n = 1:N
    x(n) = -(L/2) + (n-1)*L/Ne;
end;

figure(3)

subplot(2,2,1)
plot(x, PHI(:,1)/max( abs(PHI(:,1) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(1))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;
 
subplot(2,2,2)
plot(x, PHI(:,2)/max( abs(PHI(:,2) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(2))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;

subplot(2,2,3)
plot(x, PHI(:,3)/max( abs(PHI(:,3) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(3))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;

subplot(2,2,4)
plot(x, PHI(:,4)/max( abs(PHI(:,4) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(4))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;


figure(4)

subplot(2,2,1)
plot(x, PHI(:,5)/max( abs(PHI(:,5) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(5))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;
 
subplot(2,2,2)
plot(x, PHI(:,6)/max( abs(PHI(:,6) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(6))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;

subplot(2,2,3)
plot(x, PHI(:,7)/max( abs(PHI(:,7) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(7))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;

subplot(2,2,4)
plot(x, PHI(:,8)/max( abs(PHI(:,8) ) ) ) ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM(8))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;