%% prog 003 tubo fuente en un extremo abierto otro extremo

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
Pl = 1.0;       %presion sonora en x = -L/2 es decir primer nodo
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

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

%eliminamos primera fila primera columna
M2 = M(2:N-1,2:N-1);    %nueva matriz de masa
FM2 = M(2:N-1,1);       %primera columna de la matriz de masa / fuerza reactiva

%% matriz de rigidea

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

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

%eliminamos primera fila primera columna
K2 = K(2:N-1,2:N-1);    %nueva matriz de rigidez
FK2 = K(2:N-1,1);       %primera columna de la matriz de rigidez / fuerza reactiva

%% problema de valores y vectores propios

%problema de valores propios
[PHI2,OMG2] = eig(K2,M2);

%concatenamos una fila de ceros en x = -L/2 por la ausencia de fuente
%concatenamos una fila de ceros en x =  L/2 por la tubo abierto
PHI2 = [zeros(1,N-2);PHI2;zeros(1,N-2)];

%frecuencias naturales o de resonancia FEM
wn2 = sqrt(diag(abs(OMG2)));
fnFEM2 = wn2/(2*pi);

%% frecuencias naturales o de resonancia teoricas TEO
fnTEO2 = zeros(N-2,1);
erro2 = zeros(N-2,1);
nn = 1:N-2;

for n = 1:N-2
    fnTEO2(n) = (n*c)/(2*L);
    erro2(n) = 100*abs(fnFEM2(n) - fnTEO2(n))/fnTEO2(n);
end

TablaErro2 = [fnFEM2,fnTEO2,erro2];

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

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

%% distribucion de presion sonora en cada frecuencia de resonancia

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

figure(3)
subplot(2,2,1)
plot(xcorde, PHI2(:,1)/max( abs( PHI2(:,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, PHI2(:,2)/max( abs( PHI2(:,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, PHI2(:,3)/max( abs( PHI2(:,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, PHI2(:,4)/max( abs( PHI2(:,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, PHI2(:,5)/max( abs( PHI2(:,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, PHI2(:,6)/max( abs( PHI2(:,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, PHI2(:,7)/max( abs( PHI2(:,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, PHI2(:,8)/max( abs( PHI2(:,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;

%% resolucion del sistema de ecuaciones en el dominio de la frecuencia

%resolveremos el sistema de ecuaciones [-w^2*M1 + K1]*P1 = F1
%las frecuencias van desde f = 1 (Hz) hasta 1000 (Hz)
%debemos ahora recorda que en el primer nodo la presion sonora P(1) = Pl

%vector de frecuencias
fini = 0;
finc = 1;
ffin = 1000;
f = fini:finc:ffin;
w = 2*pi*f;
NF = length(f);

%prelocalizamos la matriz de solucion donde la columna representa frecuencia
%y la fila posicion 
P2 = zeros(N-2,NF);

for nf = 1:NF
    F2 = w(n)^2*FM2*Pl - FK2*Pl;
    P2(:,nf) = (-w(nf)^2*M2 + K2)\F2;
end

%colocamos una fila de condicion de contorno primer nodo P(1,1:NF) = Pl = 1
%colocamos una fila de condicion de contorno primer nodo P(N,1:NF) = 0
Pini(1,1:NF) = Pl;
Pfin(1,1:NF) = 0;
P2 = [Pini;P2;Pfin];

%% resultados

[Fmesh,Xmesh] = meshgrid(f,x);

figure(5)
surf(Fmesh,Xmesh,abs(P2),'FaceColor','interp','EdgeColor','none','FaceLighting','phong');
title('modulo de la presion sonora - tubo cerrado - FEM');
xlabel('f (Hz)');
ylabel('x (x)');
zlabel('|P(f,x)| (Pa)');
axis([0,f(NF),-L/2,L/2]);
grid on;
box on;
colormap jet
colorbar

%comparamos a una frecuencia fija f = 500 (Hz)
Pf500FEM = P2(:,501);

%presion teorica
f500 = 500;
k500 = 2*pi*f500/c;

Pf500TEO = Pl*sin(k500*(L/2-x))/sin(k500*L);


figure(6)
plot(xcorde,Pf500FEM,'b', xcorde,Pf500TEO,'r -o');
title('comparación presión sonora f = 500 (Hz) FEM vs TEO');
xlabel('x (m)');
ylabel('p(x) (Pa)');
legend('p(x)');
axis([-L/2,L/2,-1.1*max(Pf500TEO),1.1*max(Pf500TEO)]);
grid on;
box on;

%presion sonora en función de la frecuencia en x = L/2
PL2 = P2(Ne,:);
LpL2 = 20*log10(PL2/2e-5/sqrt(2));

figure(7)
plot(f,LpL2,'b');
title('nivel de presion sonora en x = L/2');
xlabel('f (Hz)');
ylabel('p(f) (Pa)');
legend('p(f)');
axis([fini,ffin,0,1.1*max(LpL2)]);
grid on;
box on;