%% prog 010 tubo cerrado fuente en un extremoimpedancia de radiacion otro con montaje

clearvars;
close all;
clc;

%este programa calcula los valores propios y vectores propios
%los primeros asociados a las frecuencias naturales o de resonancia 
%y los segundos a la distribución de la presion sonora en el tubo cuando
%hay un sonido cuya frecuencia coincide con las frecuencias de resonancia

%% inicializacion

L = 1.0;        %longitud del tubo
a = 0.01;       %lado/radio del tubo
A = pi*a^2;     %area de la seccion transversal del tubo A
ro = 1.18;      %densidad del aire
c = 344.0;      %velocidad del sonido
Ne = 100 ;       %numero de elementos
N = Ne + 1;     %numero de nodos
Ncord = 2;      %numero de cordenadas por elemento
Nnod = 2;       %numero de nodos por elemento
Ndof = 1;       %numero de grados de libertad por nodo
ax = L/(2*Ne);  %longitud de elemento
PL = 1.0;       %fuente de presion en x = -L/2

%% matriz de coordenadas

%inicializacion
Xcord = zeros(N,3);

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

Xcord(:,1) = xcorde;

%% matriz de conectividad

Xconec = zeros(Ne,Ncord);
for n = 1:Ne
    Xconec(n,1) = n;
    Xconec(n,2) = n + 1;
end;

%% matriz de condiciones de Dirichelet

XID = zeros(N,Ndof);    %matriz de condiciones de Dirichelet
nodDr = ones(N,Ndof);

%nodo final y quinto nodo abierto salen cosas raras
% nodDr(N,1) = 0;
% nodDr(5,1) = 0;

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

M = zeros(N,N);         %inicializacion de la matriz de masa global con ceros
K = zeros(N,N);         %inicializacion de la matriz de rigidez global con ceros

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


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

%matriz de amortiguamiento con condicion de contorno de drichelet no nula
C1 = C(2:N,2:N);
 
%% fuente sonora

%en este caso implica eliminar primera fila y primera columna del sistema
%de ecuaciones
%usar parte de la primera fila / columna de cada matriz para armar la
%fuerza reactiva

%eliminamos primera fila primera columna
M1 = M(2:N,2:N);    %nueva matriz de masa
FM1 = M(2:N,1);     %primera columna de la matriz de masa / fuerza reactiva

%eliminamos primera fila primera columna
K1 = K(2:N,2:N);    %nueva matriz de rigidez
FK1 = K(2:N,1);     %primera columna de la matriz de rigidez / fuerza reactiva

%% problema de valores y vectores propios
% 
% %problema de valores propios
% [PHI1,OMG1] = eig(K1,M1);
% 
% %concatenamos una fila de ceros en x = L/2 por la ausencia de fuente
% PHI1 = [zeros(1,N-1);PHI1];
% 
% %frecuencias naturales o de resonancia FEM
% wn1 = sqrt(diag(abs(OMG1)));
% fnFEM1 = wn1/(2*pi);
% 
% %% frecuencias naturales o de resonancia teoricas TEO
% fnTEO1 = zeros(N-1,1);
% erro1 = zeros(N-1,1);
% nn = 1:N-1;
% 
% for n = 1:N-1
%     fnTEO1(n) = ((2*n-1)*c)/(4*L);
%     erro1(n) = 100*abs(fnFEM1(n) - fnTEO1(n))/fnTEO1(n);
% end
% 
% erro(1) = 0;
% TablaErro1 = [fnFEM1,fnTEO1,erro1];
% 
% figure(1)
% plot(fnFEM1,fnTEO1,'b -o');
% title('frecuencias naturales FEM vs frecuencias naturales TEO');
% xlabel('f_{FEM} (Hz)');
% ylabel('f_{TEO} (Hz)');
% legend('f_n');
% axis([min(fnFEM1),max(fnFEM1),min(fnTEO1),max(fnTEO1)]);
% grid on;
% box on;
% 
% figure(2)
% plot(nn,erro1,'b -o');
% title('error frecuencias naturales FEM vs frecuencias naturales TEO');
% xlabel('n');
% ylabel('error (%)');
% legend('error');
% axis([min(nn),max(nn),0,max(erro1)]);
% grid on;
% box on;
% 
%% distribucion de presion sonora

dx = L/Ne;                  %longitud elemento
xcorde = (0:N-1)*dx - L/2;  %ccordenadas en eje x

% figure(3)
% subplot(2,2,1)
% plot(xcorde, PHI1(:,1)/max( abs( PHI1(:,1) ) ),'b');
% title('primer modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_1');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% subplot(2,2,2)
% plot(xcorde, PHI1(:,2)/max( abs( PHI1(:,2) ) ),'b');
% title('segundo modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_2');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% subplot(2,2,3)
% plot(xcorde, PHI1(:,3)/max( abs( PHI1(:,3) ) ),'b');
% title('tercer modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_3');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% subplot(2,2,4)
% plot(xcorde, PHI1(:,4)/max( abs( PHI1(:,4) ) ),'b');
% title('cuarto modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_4');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% 
% figure(4)
% subplot(2,2,1)
% plot(xcorde, PHI1(:,5)/max( abs( PHI1(:,5) ) ),'b');
% title('quinto modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_5');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% subplot(2,2,2)
% plot(xcorde, PHI1(:,6)/max( abs( PHI1(:,6) ) ),'b');
% title('sexto modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_6');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% subplot(2,2,3)
% plot(xcorde, PHI1(:,7)/max( abs( PHI1(:,7) ) ),'b');
% title('septimo modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_7');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 
% subplot(2,2,4)
% plot(xcorde, PHI1(:,8)/max( abs( PHI1(:,8) ) ),'b');
% title('octavo modo FEM normalizado');
% xlabel('x (m)');
% ylabel('p (Pa)');
% legend('PHI_8');
% axis([-L/2,L/2,-1.5,1.5]);
% grid on;
% box on;
% 

%% resolucion del sistema de ecuaciones en el dominio de la frecuencia

%resolveremos el sistema de ecuaciones [-w^2*M1 + K1]*P1 = F1
%las frecuencias van desde f = 1 (Hz) hasta 1000 (Hz)
%debemos ahora recorda que en el primer nodo la presion sonora P(1) = Pl

%vector de frecuencias
fini = 0;
finc = 1;
ffin = 1000;
f = fini:finc:ffin;
w = 2*pi*f;
NF = length(f);

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


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

x = xcorde;

k = 2*pi*f(600)/c;
pcer = PL*cos(k*(L/2 - x))/cos(k*L);
pabr = 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(P1),'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(P1)/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)
%pcer(x,t) = PL*cos[k*(L/2-x)]/cos(k*L)
%pabr(x,t) = PL*sin[k*(L/2-x)]/sin(k*L)

figure(7)
plot(x,abs(pcer),'b',x,abs(pabr),'k',x,abs(P1(:,600)),'r -o');
title(['comparacion solución TEO vs FEM frecuancia f = ',num2str(f(600))] );
xlabel('x (m)');
ylabel('p (Pa)');
legend('p_{TEO_{CER}}(x,f)','p_{TEO_{ABR}}(x,f)','p_{FEM_{DIR}}(x,f)');
axis([-L/2,L/2,-0.1*max(abs(pcer)),2.0*max(abs(pcer))]);
grid on;
box on;