%% prog 001 tubo cerrado en ambos extremos

clear all;
close all;
clc;

%este programa calcula los valores propios y vectores propios
%los primeros asociados a las frecuencias naturales o de resonancia 
%y los segundos a la distribución de la presion sonora en el tubo cuando
%hay un sonido cuya frecuencia coincide con las frecuencias de resonancia

%% inicializacion

L = 1.0;        %longitud del tubo
A = 1e-4;       %area de la seccion transversal
ro = 1.18;      %densidad del aire
c = 344.0;      %velocidad del sonido
Ne = 100;       %numero de elementos
N = Ne + 1;     %numero de nodos
ax = L/(2*Ne);  %longitud de elemento

%% matriz de masa

M = zeros(N,N);     %inicializacion de la matriz de masa con ceros

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

%incializacion
M(1  ,1  ) = 2;
M(N-1,N  ) = 1;
M(N  ,N-1) = 1;
M(N  ,  N) = 2;

%ciclo principal
for n = 2:N-1
    M(n  ,n)   = 4;     %diagonal principal
    M(n  ,n-1) = 1;     %sub diagonal
    M(n-1,  n) = 1;     %super diagonal
end

M = (A*ax)/(3*c^2)*M;

%% matriz de rigidea

K = zeros(N,N);     %inicializacion de la matriz de rigidez con ceros

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

%incializacion
K(1  ,1  ) =  1;
K(N-1,N  ) = -1;
K(N  ,N-1) = -1;
K(N  ,  N) =  1;

%ciclo principal
for n = 2:N-1
    K(n  ,n)   =  2;    %diagonal principal
    K(n  ,n-1) = -1;    %sub diagonal
    K(n-1,  n) = -1;    %super diagonal
end

K = A/(2*ax)*K;


%% problema de valores y vectores propios
%   d^2p
%M*----- + K*p = 0
%   dt^2
%
%asumimos la solucion p(t) = phi*exp(j*w*t) reemplazamos en la solucion
%
% -w^2*M*phi*exp(j*w*t) + K*phi*exp(j*w*t) = 0
%
%podemos simplificar exp(j*w*t) y obtenemos el problema de valores propios
%
% (w^2*M* + K)*phi = 0
%
%para resolver el problema el determinante  de la matriz debe ser nulo
%
% det(w^2*M* + K) = 0 
%
%resolvemos el determinante y obtenemos un polinomio de orden N en w^2
%
% aN*(w^2)^N + aN-1*(w^2)^(N-1) + ... + a2*(w^2)^(2)+ a1*(w^2)^(1) + a0 = 0
%
%al encontrar las raices de este polinomio se determinan  los wn (n = 1,2,...,N)
%que son %las frecuencias de resonancia angulares o frecuencias naturales
%angulares, por supuesto las frecuencias naturales fn = wn/(2*pi)
%
%nos replanteamos los sistemas de ecuaciones
%
%f1 --> w1 --> ((w1)^2*M* + K)*phi1 = 0
%f2 --> w2 --> ((w2)^2*M* + K)*phi2 = 0
%......................................
%fN --> wN --> ((wN)^2*M* + K)*phiN = 0
%
%para poder resover este problema en forma efieciente dentro del contexto
%de matlab o octave se usa el comando [PHI,OMG] = eig(K,M);

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

%frecuencias naturales o de resonancia FEM
wn = sqrt(diag(abs(OMG)));
fnFEM = wn/(2*pi);

%% frecuencias naturales o de resonancia teoricas TEO
fnTEO = zeros(N,1);
erro = zeros(N,1);
nn = 1:N;

for n = 1:N
    fnTEO(n) = ((n-1)*c)/(2*L);
    erro(n) = 100*abs(fnFEM(n) - fnTEO(n))/fnTEO(n);
end

erro(1) = 0;
TablaErro = [fnFEM,fnTEO,erro];

figure(1)
plot(fnFEM,fnTEO,'b -o');
title('frecuencias naturales FEM vs frecuencias naturales TEO');
xlabel('f_{FEM} (Hz)');
ylabel('f_{TEO} (Hz)');
legend('f_n');
axis([min(fnFEM),max(fnFEM),min(fnTEO),max(fnTEO)]);
grid on;
box on;

figure(2)
plot(nn,erro,'b -o');
title('error frecuencias naturales FEM vs frecuencias naturales TEO');
xlabel('n');
ylabel('error (%)');
legend('error');
axis([min(nn),max(nn),0,max(erro)]);
grid on;
box on;

%% distribucion de presion sonora

dx = L/Ne;                  %longitud elemento
xcorde = (0:N-1)*dx - L/2;  %ccordenadas en eje x

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

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

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

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


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

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

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

subplot(2,2,4)
plot(xcorde, PHI(:,8)/max( abs( PHI(:,8) ) ),'b');
title('octavo modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
legend('PHI_8');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;