%% Prog004FEMTuboFuenteAbiertoSinMontaje

clearvars;
close all;
clc;

%este es un programa de elementos finitos que representa un tubo abierto en un extremo 
%con fuete y abierto en el otro extremo
%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 despues determinar la respuesta forzada

%% 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 masa con condicion de contorno de drichelet no nula
M2 = M(2:N-1,2:N-1);
 
 %% 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; 

%matriz de masa con condicion de contorno de drichelet no nula
K2 = K(2:N-1,2:N-1);
 
%% valores y vectores propios
 
%la función eig() determina los valores y vectores propios de la matriz
[PHI2,LAM2] = eig(K2,M2);

%incorporamos la condicion de contorno de tubo abierto con la fuente
%apagada para ello incorporamos una fila llena de ceros e incorporamos otra
%fila de ceros al final representando el tubo abierto
PHI2 = [zeros(1,N-2);PHI2;zeros(1,N-2)];
 
%frecuencias naturales elementos finitos FEM
f_FEM2 = diag( real(sqrt(LAM2)))/(2*pi);
nn = 1:length(f_FEM2);

%frecuencias naturales teóricas TEO
f_TEO2 = zeros(N-2,1);
for n = 1:N-2
     f_TEO2(n) = n*c/(2*L);
end;
 
%error
err2 = 100*abs(f_FEM2 - f_TEO2)./f_TEO2;
 
%tabla
f_FEM2_f_TEO2_err2 = [f_FEM2, f_TEO2, err2];
 
%% figuras frecuencias valores propios

figure(1)
plot(f_TEO2,f_FEM2,'b -o');
title('frecuencias f_{TEO} vs f_{FEM}');
xlabel('f_TEO (Hz)');
ylabel('f_TEO (Hz)');
axis([min(f_TEO2),max(f_TEO2),min(f_FEM2),max(f_FEM2)]);
grid on;
box on;

figure(2)
plot(nn,err2,'b -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, PHI2(:,1)/max( abs(PHI2(:,1) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,2)/max( abs(PHI2(:,2) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,3)/max( abs(PHI2(:,3) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,4)/max( abs(PHI2(:,4) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,5)/max( abs(PHI2(:,5) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,6)/max( abs(PHI2(:,6) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,7)/max( abs(PHI2(:,7) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(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, PHI2(:,8)/max( abs(PHI2(:,8) ) ) ,'b') ;
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM2(8))] );
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2,L/2,-1.5,1.5]);
grid on;
box on;

%% solucion forzada dominio de la frecuencia metodo directo

%resolveremos el sistema de ecuaciones [-w^2*M1 + K1]*P1(w) = F1(w)
%donde w = 2*pi*f y f = 1,2,3,...,1000 (Hz) 
%la pesion sonora en x = -L/2 es PL = exp(-j*w*t) con PL = 1 (Pa)

%componentes del vector de fuerzas
PL = 1;             %amplitud de presion sonora
Pm = M(2:N-1,1);    %primera columna de la matriz de masa
Pk = K(2:N-1,1);    %primera columna de la matriz de rigidez

%frecuencias
fmin = 1;
fmax = 1000;
finc = 1;
f = fmin:finc:fmax;     %vector de frecuencias
w = 2*pi*f;             %vector de frecuencias angulares
NF = length(f);

%prelocalizamos la solucion como matriz
%la primera columna tiene la solución para f = 1 (Hz)
%la segunda columna tiene la solución para f = 2 (Hz)
%...................................................
%la ultima columna tiene la solución para f = 1000 (Hz)
    
P2 = zeros(N-2,NF);

%solucion
for n = 1:NF
    
    %vector de fuerza dependiente de la frecuencia
    F2 = w(n)^2*PL*Pm - PL*Pk;
    
    %usamos A\b que es mejor que inv(A)*b
    P2(:,n) = (-w(n)^2*M2 + K2)\F2;
end;

%incorporamos la solucion PL para x = -L/2
P2 = [PL*ones(1,NF);P2;zeros(1,NF)];

%% comparación resultados teóricos y elementos finitos

k = 2*pi*f(1000)/c;
p = PL*sin(k*(L/2 - x))/sin(k*L);

%grafico espacio - frecuencia - presion sonora
%distribucion del valor absoluto de la presion sonora ontenida usando FEM
%en funcion de la posicion "x" y de la frecuencia "f"

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

figure(5)
surf(Fmesh,Xmesh,abs(P2),'FaceColor','interp',...
    'EdgeColor','none',...
    'FaceLighting','phong')
title('Modulo Presion Sonora - Tubo Cerrado - FEM');
xlabel('f (Hz)');
ylabel('x (m)');
zlabel('|P| (Pa)')
axis([0,f(NF),-L/2,L/2]);
grid on;
box on;
colormap jet
colorbar;

figure(6)
surf(Fmesh,Xmesh,20*log10(abs(P2)/sqrt(2)/2e-5),'FaceColor','interp',...
    'EdgeColor','none',...
    'FaceLighting','phong')
title('Nivel de Presion Sonora - Tubo Cerrado - FEM');
xlabel('f (Hz)');
ylabel('x (m)');
zlabel('|P| (Pa)')
axis([0,f(NF),-L/2,L/2]);
grid on;
box on;
colormap jet
colorbar;

%compararemos con la solucion exacta para la frecuencia de f = 1000 (Hz)
%p(x,t) = PL*sin[k*(x-L/2)]/sin(k*L)

figure(7)
plot(x,p,'b',x,P2(:,1000),'r -o');
title(['comparacion solución TEO vs FEM frecuancia f = ',num2str(f(1000))] );
xlabel('x (m)');
ylabel('p (Pa)');
legend('p_{TEO}(x,f)','p_{FEM_DIR}(x,f)');
axis([-L/2,L/2,-1.5*max(p),1.5*max(p)]);
grid on;
box on;