%% Prog010FEMTuboFuenteCerradoImpRadConMontaje

clearvars;
close all;
clc;

%este es un programa de elementos finitos que representa un tubo abierto en un extremo
%con fuete y cerrado 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
Ne = 100;       %numero de elementos
NnD = 0;        %numero de DOF s con condiciones drichelet
N = Ne+1;       %numero de DOFs - funciones de interpolacion
Nt = Ne+1-NnD;  %numero de DOFs - funciones de interpolacion con drichelet
Nnod = 2;       %numero de nodos elemento
Ncord = 2;      %numero de coordenadas por elemento
Ndof = 1;       %numero de dof por nodo
L = 1.0;        %longitud tubo m
ax = L/(2*Ne);  %logitud de 1/2 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
ro = 1.18;      %densidad del aire
B = c^2*ro;     %rigidez de compresibilidad adiabática de volumen
P0 = 2e-5;

%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

%se utilizara matrices elementales

%% matriz de cordenadas

%inicio
Xcord = zeros(N,3);

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

Xcord(:,1) = xcorde;

%% matriz de conectividad

%inicio
Xconec = zeros(Ne,Ncord);

for n = 1:N-1
    Xconec(n,1) = n;
    Xconec(n,2) = n+1;
end

%% matriz de ID

%inicio
XID = zeros(N,Ndof);

%nodos con condicion de drichelet (ninguno en este caso)
nodDr = ones(N,Ndof);

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 - rigidez elemental y global

%inicializacion matriz de masa global
M = zeros(Nt,Nt);

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

%inicializacion matriz de rigidez global
K = zeros(Nt,Nt);

%matriz de 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

%% matriz de amortiguamiento

C = zeros(N,N);
C(N,N) = 1;

%% correccion de fuente

M1 = M(2:Nt,2:Nt);
K1 = K(2:Nt,2:Nt);
C1 = C(2:Nt,2:Nt);

% %% valores y vectores propios
%
% %la función eig() determina los valores y vectores propios de la matriz
% [PHI1,LAM1] = eig(K1,M1);
%
% %incorporamos la condicion de contorno de tubo abierto con la fuente
% %apagada para ello incorporamos una fila llena de ceros
% PHI1 = [zeros(1,N-1);PHI1];
%
% %frecuencias naturales elementos finitos FEM
% f_FEM1 = diag( real(sqrt(LAM1)))/(2*pi);
% nn = 1:length(f_FEM1);
%
% %frecuencias naturales teóricas TEO
% f_TEO1 = zeros(N-1,1);
% for n = 1:N-1
%      f_TEO1(n) = (2*n-1)*c/(4*L);
% end;
%
% %error
% err1 = 100*abs(f_FEM1 - f_TEO1)./f_TEO1;
%
% %tabla
% f_FEM1_f_TEO1_err1 = [f_FEM1, f_TEO1, err1];
%
% %% figuras frecuencias valores propios
%
% figure(1)
% plot(f_TEO1,f_FEM1,'b -o');
% title('frecuencias f_{TEO} vs f_{FEM}');
% xlabel('f_TEO (Hz)');
% ylabel('f_TEO (Hz)');
% axis([min(f_TEO1),max(f_TEO1),min(f_FEM1),max(f_FEM1)]);
% grid on;
% box on;
%
% figure(2)
% plot(nn,err1,'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, PHI1(:,1)/max( abs(PHI1(:,1) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,2)/max( abs(PHI1(:,2) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,3)/max( abs(PHI1(:,3) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,4)/max( abs(PHI1(:,4) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,5)/max( abs(PHI1(:,5) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,6)/max( abs(PHI1(:,6) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,7)/max( abs(PHI1(:,7) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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, PHI1(:,8)/max( abs(PHI1(:,8) ) ) ,'b') ;
% title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(f_FEM1(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 +j*w*C1 + 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 = 1e-3;          %amplitud de presion sonora
Pm = M(2:N,1);      %primera columna de la matriz de masa
Pc = C(2:N,1);      %primera columna de la matriz de masa
Pk = K(2:N,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)

P1 = zeros(N-1,NF);

%solucion
for n = 1:NF
    %impedancia de radiacion
    %ZL = ro*c*( 1-2*besselj(1,2*k*a) + 2*j*struve(1,2*k*a) )
    k = w(n)/c;
    ka2 = 2*k*a;
    RL = 1 - 2*besselj(1,ka2)/ka2;
    XL = 2*struve(1,ka2)/ka2;
    ZL = ro*c*( RL + j*XL);

    %impedancia arbitraria
    %ZL = ro*c*( 0.1 + j*0.01 );

    C1(N-1,N-1) = ro*A/ZL;

    %vector de fuerza dependiente de la frecuencia
    F1 = w(n)^2*PL*Pm - j*w(n)*Pc - PL*Pk;

    %usamos A\b que es mejor que inv(A)*b
    P1(:,n) = (-w(n)^2*M1 +j*w(n)*C1 + K1)\F1;
end

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

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

k = 2*pi*f(500)/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(1)
surf(Fmesh,Xmesh,abs(P1),'FaceColor','interp',...
    'EdgeColor','none',...
    'FaceLighting','phong')
title('Modulo Presion Sonora - Tubo Impedancia Radiacion - 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(2)
surf(Fmesh,Xmesh,20*log10((abs(P1) + P0)/sqrt(2)/2e-5),'FaceColor','interp',...
    'EdgeColor','none',...
    'FaceLighting','phong')
title('Nivel de Presion Sonora - Tubo Impedancia Radiacion - 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(3)
plot(x,abs(p),'b',x,abs(P1(:,500)),'r -o');
title(['comparacion solucion TEO vs FEM frecuencia f = ',num2str(f(500))] );
xlabel('x (m)');
ylabel('p (Pa)');
legend('p_{TEO}(x,f)','p_{FEM_DIR}(x,f)');
axis([-L/2,L/2,0.0,1.5*max(p)]);
grid on;
box on;

figure(4)
plot(f,20*log10(P1(1,:)/sqrt(2)/P0),'b',f,20*log10(P1(Nt,:)/sqrt(2)/P0),'r');
title('niveles de presion sonora al incio y al final del tubo');
xlabel('f (Hz)');
ylabel('L_p (dB)');
legend('L_{p1}(f)','L_{p2}(f)');
grid on;
box on;

figure(5)
plot(f,20*log10( P1(1,:)./P1(Nt,:)),'b');
title('reduccion de ruido NR');
xlabel('f (Hz)');
ylabel('NR(f) (dB)');
legend('NR(f)');
grid on;
box on;

figure(6)
plot(f,20*log10( P1(Nt,:)./P1(1,:)),'b');
title('respuesta de frecuencia');
xlabel('f (Hz)');
ylabel('H(f) (dB)');
legend('H (f)');
grid on;
box on;