%% prog 002 tubo fuente en un extremo cerrado 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
M1 = M(2:N,2:N);    %nueva matriz de masa
FM1 = M(2:N,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
K1 = K(2:N,2:N);    %nueva matriz de rigidez
FK1 = K(2:N,1);     %primera columna de la matriz de rigidez / fuerza reactiva

%% problema de valores y vectores propios

%problema de valores propios
[PHI1,OMG1] = eig(K1,M1);

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

%frecuencias naturales o de resonancia FEM
wn1 = sqrt(diag(abs(OMG1)));
fnFEM1 = wn1/(2*pi);

%% frecuencias naturales o de resonancia teoricas TEO
fnTEO1 = zeros(N-1,1);
erro1 = zeros(N-1,1);
nn = 1:N-1;

for n = 1:N-1
    fnTEO1(n) = ((2*n-1)*c)/(4*L);
    erro1(n) = 100*abs(fnFEM1(n) - fnTEO1(n))/fnTEO1(n);
end

TablaErro1 = [fnFEM1,fnTEO1,erro1];

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

figure(2)
plot(nn,erro1,'b -o');
title('error frecuencias naturales FEM vs frecuencias naturales TEO');
xlabel('n');
ylabel('error (%)');
legend('error');
axis([min(nn),max(nn),0,max(erro1)]);
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, PHI1(:,1)/max( abs( PHI1(:,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, PHI1(:,2)/max( abs( PHI1(:,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, PHI1(:,3)/max( abs( PHI1(:,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, PHI1(:,4)/max( abs( PHI1(:,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, PHI1(:,5)/max( abs( PHI1(:,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, PHI1(:,6)/max( abs( PHI1(:,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, PHI1(:,7)/max( abs( PHI1(:,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, PHI1(:,8)/max( abs( PHI1(:,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 
P1 = zeros(N-1,NF);

for nf = 1:NF
    F1 = w(n)^2*FM1*Pl - FK1*Pl;
    P1(:,nf) = (-w(nf)^2*M1 + K1)\F1;
end

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

%% resultados

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

figure(5)
surf(Fmesh,Xmesh,abs(P1),'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 = P1(:,501);

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

Pf500TEO = Pl*cos(k500*(x-L/2))/cos(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.5,1.5]);
grid on;
box on;

%presion sonora en función de la frecuencia en x = L/2
PL2 = P1(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;