%% prog 007 tubo cerrado en ambos extremos con montaje seccion exponencial

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
A0 = 1e-4;      %area de la seccion transversal
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
beta = 0.025;     %crecimiento exponencial de la seccion transversal

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

%loop principal montaje
for nel = 1:Ne
    
    A(nel) = A0*exp(beta*nel);
    
    %matriz masa elemental
    Me = [ 2, 1 ;
           1, 2];
    Me = (A(nel)*ax)/(3*c^2)*Me;

    
    %matriz rigidez elemental
    Ke = [ 1,-1 ;
          -1, 1];
    Ke = A(nel)/(2*ax)*Ke;

    
    
    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;    


%% problema de valores y vectores propios
%   d^2p
%M*----- + K*p = 0
%   dt^2
%
%asumimos la solucion p(t) = phi*exp(j*w*t) reemplazamos en la solucion
%
% -w^2*M*phi*exp(j*w*t) + K*phi*exp(j*w*t) = 0
%
%podemos simplificar exp(j*w*t) y obtenemos el problema de valores propios
%
% (w^2*M* + K)*phi = 0
%
%para resolver el problema el determinante  de la matriz debe ser nulo
%
% det(w^2*M* + K) = 0 
%
%resolvemos el determinante y obtenemos un polinomio de orden N en w^2
%
% aN*(w^2)^N + aN-1*(w^2)^(N-1) + ... + a2*(w^2)^(2)+ a1*(w^2)^(1) + a0 = 0
%
%al encontrar las raices de este polinomio se determinan  los wn (n = 1,2,...,N)
%que son %las frecuencias de resonancia angulares o frecuencias naturales
%angulares, por supuesto las frecuencias naturales fn = wn/(2*pi)
%
%nos replanteamos los sistemas de ecuaciones
%
%f1 --> w1 --> ((w1)^2*M* + K)*phi1 = 0
%f2 --> w2 --> ((w2)^2*M* + K)*phi2 = 0
%......................................
%fN --> wN --> ((wN)^2*M* + K)*phiN = 0
%
%para poder resover este problema en forma efieciente dentro del contexto
%de matlab o octave se usa el comando [PHI,OMG] = eig(K,M);

%problema de valores propios
[PHI,OMG] = eig(K,M);

%frecuencias naturales o de resonancia FEM
wn = sqrt(diag(abs(OMG)));
fnFEM = wn/(2*pi);

%% frecuencias naturales o de resonancia teoricas TEO
fnTEO = zeros(N,1);
erro = zeros(N,1);
nn = 1:N;

for n = 1:N
    fnTEO(n) = ((n-1)*c)/(2*L);
    erro(n) = 100*abs(fnFEM(n) - fnTEO(n))/fnTEO(n);
end

erro(1) = 0;
TablaErro = [fnFEM,fnTEO,erro];

figure(1)
plot(fnFEM,fnTEO,'b -o');
title('frecuencias naturales FEM vs frecuencias naturales TEO');
xlabel('f_{FEM} (Hz)');
ylabel('f_{TEO} (Hz)');
legend('f_n');
axis([min(fnFEM),max(fnFEM),min(fnTEO),max(fnTEO)]);
grid on;
box on;

figure(2)
plot(nn,erro,'b -o');
title('error frecuencias naturales FEM vs frecuencias naturales TEO');
xlabel('n');
ylabel('error (%)');
legend('error');
axis([min(nn),max(nn),0,max(erro)]);
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, PHI(:,1)/max( abs( PHI(:,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, PHI(:,2)/max( abs( PHI(:,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, PHI(:,3)/max( abs( PHI(:,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, PHI(:,4)/max( abs( PHI(:,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, PHI(:,5)/max( abs( PHI(:,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, PHI(:,6)/max( abs( PHI(:,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, PHI(:,7)/max( abs( PHI(:,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, PHI(:,8)/max( abs( PHI(:,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;

figure(5)
bar(A,'b');
hold on;
bar(-A,'b');
title('figura del area del tubo');
xlabel('x (cm)');
ylabel('A (m^2)');