%% prog 002 FEM tubo fuente en un extremo cerrado en otro extremos sin montaje

clearvars;
close all;
clc;

%este programa calcula los valores propios (frecuencias naturales) y 
%vectores propios (formas de vibracion) de un tubo con una fuente en un extremo
%y cerrado en el otro extremos de seccion transversal constante A longitud 
%L a una velocidad del sonido constante c
%inicialmente el numero de elementos es 100 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 = 343;        %velocidad del sonido
ro = 1.18;      %densidad del aire
Pin = 1.0;      %presion de la fuente ubicada en x = -L/2 (m)

%% matriz de masa

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

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

%aplicaremos la condicion de contorno para la matriz de masa, es decir
%construiremos una neva matriz eliminando la primera fila y la primera
%columna

M1 = M(2:N,2:N);	%nueva matriz de masa que considera la condicion de 
                    %contorno en el nodo 1

%% matriz de rigidez

%inicializacion
K = zeros(N,N);

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

%aplicaremos la condicion de contorno para la matriz de rigidez, es decir
%construiremos una neva matriz eliminando la primera fila y la primera
%columna

K1 = K(2:N,2:N);	%nueva matriz de rigidez que considera la condicion de 
                    %contorno en el nodo 1


%% 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
%M1 es matriz (N-1 X N-1)
%K1 es matriz (N-1 X N-1)
%phi1 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)

%LAM es la matriz diagonal formada por los valores propios
%LAM = diag(lam1, lam2, ..., lamN-1) 
%LAM = 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('fTEO1 (Hz)');
ylabel('fFEM1 (Hz)');
grid on;
box on;

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

TablaError1 = [fTEO1 fFEM1 errof1];

%% resultados distribucion normalizada a la unidad (1) de presion sonora 

%resultados distribucion normalizada a la unidad (1) de presion sonora 
%para las frecuecias de resonancia

%cordenadas
dx = L/Ne;  %incremento o longitud de cada elemento
xcorde = (0:N-1)*dx - L/2; %nos da las coordenadas desde -L/2 a L/2 

%con el hecho que el primer nodo del tubo posee valor de presion sonora
%prescrito no se debe incluir en los graficos normalizados de la
%distribucion de la presion sonora para cada una de las frecuencias
%naturales, por esa causa eliminaremos la primera coordenada del proceso de
%graficacion

%grafico
figure(3)

subplot(2,2,1)
plot( xcorde(2:N), 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(2:N), 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(2:N), 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(2:N), 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;    
    
%% resolucion del problema [-w^2*M1 + K1]P1(w) = F1(w)
%resolveremos el problema en el dominio de la frecuencia
%desde f = 1(Hz) hasta 1000(Hz)
%debemos recordar que en terminos de condicion de contoreno el valor de la
%presion sonora en el primer nodo Pin = 1 (Pa)

%arnmaremos el vector de fuerzas externas usando parte de la primera
%columna de la matriz de masa y la primera columna de la matriz de rigidez
Fa = M(2:N,1);  %parte de la primera columna de la matriz de masa
Fb = K(2:N,1);  %parte de la primera columna de la matriz de rigidez

%vector de frecuencias
fmin = 1;       %frecuencia inicial
finc = 1;       %incremento de frecuencia
fmax = 1000;    %frecuencia final

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

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

for nf = 1:NF
    
    F1 = w(nf)^2*Pin*Fa - Pin*Fb;
    P1(:,nf) = (-w(nf)^2*M1 + K1)\F1;
    
end;

%esta solucion está incompleta falta incorporar la condicion de contorno
%inicial es decir en el incio del tubo x = L/2 (m) p(L/2) = Pin o dico de
%otra forma P(1) = Pin para el nodo inicial
%la primera fila debe corresponder a P(1,1:NF) = Pin y adosarse a la matriz
%P1 para formar la solucion completa y poder graficarse
Pini(1,1:NF) = Pin;
P = [Pini;P1];

%% resultados

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

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

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

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

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

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

%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