%% prog 004 FEM tubo fuente en un extremo cerrado cin impedancia otro extremo sin montaje

clearvars;
close all;
clc;

%este programa calcula la distribución de la presion sonora de untubo 
%con una fuente en un extremo y cerrado en el otro extremo con una 
%impedancia arbitraria.
%el tubo es 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 = 500;       % 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
rd = 0.01;      %radio del tubo (m)
A = pi*rd^2;    %area de la sección transversal  (m^2)
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

%% matriz de amortiguamiento

%inicializacion
C = zeros(N,N);

%impedancia
zl = 2000 + j*500;
yl = 1/zl;

C(N,N) = ro*yl*A;

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

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


                    
    
%% resolucion del problema [-w^2*M1 +j*w*C + 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
Fc = C(2:N,1);  %parte de la primera columna de la matriz de amortiguamiento

%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 - j*w(nf)*Pin*Fc - Pin*Fb;
    P1(:,nf) = (-w(nf)^2*M1 + j*w(nf)*C1 + 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

%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 

%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(1)
plot(x,abs(pteo),'b',x,abs(P(:,500)),'r -o')
title('comparacion resultados teoricos vs FEM con impedancia zl');
xlabel('x (m)');
ylabel('|p| (Pa)');
legend('pteo','pFEM zl');
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(2)
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