%% prog 002 FEM tubo fuente en extremo 1 cerrado 2 sin montaje

clearvars;
close all;
clc;

%este programa calcula los valores propios (frecuencias naturales) y 
%vectores propios (formas de vibracion) de un tubo cerrado en ambos
%extremos de seccion transversal constante A longitud L a una velocidad del
%sonido constante c
%inicialmente el numero de elementos es 4 pero será modificado a un numero
%variable de elementos

%% inicializacion

Ne = 100;       %numero de elementos
N = Ne + 1;     %numero de nodos / grados de libertad / DOF
L = 1.0;        %longitud del tubo en metros (m)
ax = L/(2*Ne);  %mitad de la longitud de un elemento
A = 0.0001;     %area de la sección transversal (m^2)( 1.0e-4) = 1 X 10^(-4) )
c = 340;        %velocidad del sonido
ro = 1.18;      %densidad del aire

%% matriz de masa

%inicializacion
M = zeros(N,N); %matriz incializada con ceros

%en una primera etapa colocaremos de forma directa todos los valores en la
%matriz 
%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];

%utilizaremos la estructura de la matriz para generalizar resultados para
%tubos de seccion cosntante donde todos los elementos tienen el mismo
%tamaño

%proceso incial
M(1  , 1  )   = 2;
M(N-1, N  )   = 1;
M(N  , N-1)   = 1;
M(N  , N  )   = 2;

%iteracion principal desde la 2da fila/columna a la penultima fila/columna
for n1 = 2:N-1
    M(n1  ,n1    ) = 4;	%diagonal llena de 4
    M(n1-1,n1    ) = 1;	%super diagonal llena de 1
    M(n1  ,n1 - 1) = 1; %sub diagonal llena de 1
end;    

%factores de ajuste
M = (A*ax)/(3*c^2)*M; 

%% matriz de rigidez

%inicializacion
K = zeros(N,N);

%en una primera etapa colocaremos de forma directa todos los valores en la
%matriz 

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

%utilizaremos la estructura de la matriz para generalizar resultados para
%tubos de seccion cosntante donde todos los elementos tienen el mismo
%tamaño

%proceso incial
K(1  , 1  )   =  1;
K(N-1, N  )   = -1;
K(N  , N-1)   = -1;
K(N  , N  )   =  1;

%iteracion principal desde la 2da fila/columna a la penultima fila/columna
for n1 = 2:N-1
    K(n1  ,n1    ) =  2;    %diagonal llena de 2
    K(n1-1,n1    ) = -1;    %super diagonal llena de -1
    K(n1  ,n1 - 1) = -1;    %sub diagonal llena de -1
end;    

%factores de ajuste
K = A/(2*ax)*K; 

%% incorporacion de la fuente en x = -L/2

%tanto para la matriz de masa como la mattriz de rigidez extraeremos las
%matrices necesarias para este nuevo sistema

M1 = M(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N
K1 = K(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N

 
%% problema de valores y vectores propios
%sirve para determinar las frecuencias naturales y la distribucion de la
%presion sonora a dichas frecuencias naturales
%resolvemos la ecuacion
%[lam*M1 + K1]*phi = 0
%M es matriz (N-1 X N-1)
%K es matriz (N-1 X N-1) 
%phi es vector (N-1 X 1)
%lam es escalar lam = -omega^2 (1 X 1)
%fin de que la solucion sea no trivial, es decir distinta de cero
%obtendremos N soluciones para lam (es decir N frecuencias naturales) y N
%soluciones para phi(es decir N distribuciones de presion sonora)

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

%PHI es la matriz (N-1 X N-1)formada por los vectores propios 
%PHI = %[phi1, phi2,..., phiN-1] los phi son vectores columna
%cada phi es (N-1 X 1)

%LAM1 es la matriz diagonal formada por los valores propios
%LAM1 = diag(lam1, lam2, ..., lamN-1) 
%LAM1 = diag(omega1^2, omega2^2, ..., omegaN-1^2)

%frecuencias naturales aproximadas por FEM
%diag estrae la diagonal principal de una matriz
fFEM1 = real( sqrt (diag(LAM1) ) )/(2*pi);

%normalizacion de vectores propios
mnorm = diag(transpose(PHI1)*M1*PHI1); 
for n1 = 1:N-1
    PHI1(:,n1) = PHI1(:,n1)/sqrt( mnorm(n1) );
end;

%% resultados frecuencias

%resultados teoricos
fTEO1 = zeros(N-1,1);
for n1 = 1:N-1
    fTEO1(n1) = (2*n1-1)*c/(4*L);
end;

%error en frecuencias
errof1 = zeros(N-1,1);
for n1 = 1:N-1
    errof1(n1) = 100*abs( (fTEO1(n1) - fFEM1(n1))/fTEO1(n1) );
end;

%grafico de frecuencias
figure(1)
plot(fTEO1,fFEM1,'-o');
title('frecuencias teoricas vs frecuencias fem');
xlabel('fTEO (Hz)');
ylabel('fFEM (Hz)');
grid on;
box on;

%grafico de error
figure(2)
plot(errof1,'-o');
title('error');
xlabel('numero de modo');
ylabel('error (%)');
grid on;
box on;

TablaError1 = [fTEO1 fFEM1 errof1];

%% resultados distribucion de presion sonora
%%%%arreglar desde acá
%cordenadas
dx = L/Ne;  %incremento o longitud de cada elemento
xcorde = (1:N-1)*dx - L/2; %nos da las coordenadas desde -L/2 a L/2 

%grafico
figure(3)

subplot(2,2,1)
plot( xcorde, PHI1(:,1)/max( abs( PHI1(:,1) ) ) );
title('1er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;

subplot(2,2,2)
plot( xcorde, PHI1(:,2)/max( abs( PHI1(:,2) ) ) );
title('2er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;

subplot(2,2,3)
plot( xcorde, PHI1(:,3)/max( abs( PHI1(:,3) ) ) );
title('3er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;    
    
subplot(2,2,4)
plot( xcorde, PHI1(:,4)/max( abs( PHI1(:,4) ) ) );
title('4er modo FEM normalizado');
xlabel('x (m)');
ylabel('p (Pa)');
axis([-L/2 L/2 -1 1]);
grid on;
box on;    
    
%% solución forzada con fuente en x = -L/2    

% resolveremos el problema [-w^2*M1 + K1]P1 = F1
%esto consiste en resolver este sistema de ecuaciones para cada valor de 
%la frecuencia angular w = 2*pi*f, donde f = 1,2,...,1000 (Hz) es la frecuencia
%la presión sonora en x = -L/2 es P(-L/2,t) = PL*exp(j*w*t) con PL = 1 (Pa)

%armaremos los vectores de fueza con las primeras columnas de masa y
%rigidez originales, es deccir antes de colocar las condiciones de contorno
PL = 1;
Pm = M(2:N,1);
Pk = K(2:N,1);

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

%prelocalizacion de los resultados de la presion sonora
%se guardan los resultados en una matriz de N-1 filas que indican la
%posicion en el eje "x" y NF columnas que son para las frecuencias

P1 = zeros(N-1,NF); %la prelocalización nos ayuda con la velocidad de ejecución

%ciclo principal de calculo
for nf = 1:NF
    
    F1 = w(nf)^2*PL*Pm - PL*Pk;     %vector de fuerza dependiente de la frecuencia
    P1(:,nf) = (-w(nf)^2*M1 + K1)\F1;   %A\b es lo mismo que inv(A)*b pero mas exacta

end;

%despues debo incorporar para todas las frecuencias la condición de
%contorno p(-L/2) = PL
Pini(1,1:NF) = PL;      %fila con valores PL para todas las frecuencias
P = [Pini;P1];          %concatemanos la fila Pini a la matriz P1

%% resultados

%compararemos resultados de la distribucion de la presion sonora
%para una frecuencia especifica que no sea de resonancia f = 900 (Hz)
%en la longitud del tubo

%posicion x a partir de x = 0
x = xcorde + L/2;
x = [0,x];
NX = length(x);

%numero de onda
k = 2*pi*f(900)/c;

%resultado teórico
for nx = 1:NX
    pteo(nx) = PL*cos(k*(L - x(nx)))/cos(k*L);
end;

figure(4)
plot(x,pteo,'b',x,P(:,900),'r -o')
title('comparacion resultados teoricos vs FEM');
xlabel('x (m)');
ylabel('p (Pa)');
legend('pteo','pFEM');
grid on;
box on;

%error cuadratico promedio
erro2prom = sum( (pteo - transpose(P(:,500)) ).^2   )/NX;

%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(P),'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),0,L]);
grid on;
box on;
colormap jet
colorbar;

figure(6)
surf(Fmesh,Xmesh,20*log10(abs(P)/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),0,L]);
grid on;
box on;
colormap jet
colorbar;


%tarea hacer lo mismo para el tubo con fuente en el incio y abierto al
%final p(-L/2)= P(1) = Pini y p(L/2) = 0