%% prog 009 tubo cerrado fuente en un extremo cerrado otro con montaje

clearvars;
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
Ncord = 2;      %numero de cordenadas por elemento
Nnod = 2;       %numero de nodos por elemento
Ndof = 1;       %numero de grados de libertad por nodo
ax = L/(2*Ne);  %longitud de elemento
PL = 1.0;       %fuente de presion en x = -L/2

%% matriz de coordenadas

%inicializacion
Xcord = zeros(N,3);

%coordenadas en x
deltaX = L/Ne;
xcorde = transpose( (0:N-1)*deltaX - L/2 );

Xcord(:,1) = xcorde;

%% matriz de conectividad

Xconec = zeros(Ne,Ncord);
for n = 1:Ne
    Xconec(n,1) = n;
    Xconec(n,2) = n + 1;
end;

%% matriz de condiciones de Dirichelet

XID = zeros(N,Ndof);    %matriz de condiciones de Dirichelet
nodDr = ones(N,Ndof);

%nodo final y quinto nodo abierto salen cosas raras
% nodDr(N,1) = 0;
% nodDr(5,1) = 0;

cont = 1;               %contador

for n = 1:N
    
    if nodDr(n) == 0
        XID(n,1) = 0;
    else
        XID(n,1) = cont;
        cont = cont + 1;
    end;
    
end;



%% matriz de masa y rigidez

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

%matriz masa elemental
Me = [ 2, 1 ;
       1, 2];
Me = (A*ax)/(3*c^2)*Me;


%matriz rigidez elemental
Ke = [ 1,-1 ;
      -1, 1];
Ke = A/(2*ax)*Ke;

%loop principal montaje
for nel = 1:Ne
    
    for noi = 1:Nnod
        for ngi = 1:Ndof
             for noj = 1:Nnod
                 for ngj = 1:Ndof
                    
                     NN = XID( Xconec(nel,noi) , ngi );
                     MM = XID( Xconec(nel,noj) , ngj );
                     
                     if (NN ~= 0) && (MM ~= 0) 
                         
                         fila = Ndof*(Nnod-1)*(noi - 1) + ngi;
                         colu = Ndof*(Nnod-1)*(noj - 1) + ngj;
                         
                         M(NN,MM) = M(NN,MM) + Me( fila, colu );
                         K(NN,MM) = K(NN,MM) + Ke( fila, colu );
                     
                     end;
                     
                 end;
             end;
        end; 
    end;    
end;    

%% fuente sonora

%en este caso implica eliminar primera fila y primera columna del sistema
%de ecuaciones
%usar parte de la primera fila / columna de cada matriz para armar la
%fuerza reactiva

%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

%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

erro(1) = 0;
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

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

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

x = xcorde;

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