%% prog fem010 tubo fuente en un extremos cerrado autovalores con montaje

%este programa calcula los valores propios de un tubo de 
%longitud L con fuente en un extremo cerrado extremo
%area de la seccion transversal constante A
%el numero de elementos es variable 


%se utilizara matrices elementales

clear all;
close all;
clc;

%% inicializacion

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 = 1e-4;       %area seccion transversal tubo cuadrado 0.01 X 0.01 m2
c = 340;        %velocidad del sonido
ro = 1.12;      %densidad aire
B = c^2*ro;     %rigidez de compresibilidad adiabática de volumen

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

% la matriz de masa es normalizada a la velocidad del sonido

Me = [2 1;1 2];
Me = (A*ax)/(3*c^2)*Me;  

%inicializacion
M = zeros(Nt,Nt);

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


%% matriz de rigidez

Ke = [1 -1;-1 1];
Ke = A/(2*ax)*Ke;

%inicializacion
K = zeros(Nt,Nt);

%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;
                         
                         K(NN,MM) = K(NN,MM) + Ke( fila, colu );
                     end;
                     
                 end;
             end;
        end; 
    end;    
end;    
    
%% incorporacion de la fuente en x = -L/2

%tanto para la matriz de masa como la matriz de rigidez extraeremos las
%matrices necesarias para este nuevo sistema

M1 = M(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N
K1 = K(2:N,2:N);%extrae de la matriz principal la submatriz delde la fila/columna 2 a la fila/columna N

%% problema de valores propios

[PHI1,LAM1] = eig(K1,M1);

%frecuencias aproximadas
fFEM1 =  real( sqrt( diag(LAM1) ) )/(2*pi);   %formato compactor

%normalizacion de los vectores propios
mnorm = diag( transpose(PHI1)*M1*PHI1 );

for n1 = 1:Nt-1
    PHI1(:,n1) = PHI1(:,n1)/mnorm(n1);
end;    

%frecuencias teoricas
for n2 = 1:Nt-1
    fTEO(n2,1) = (2*n2-1)*c/(4*L);
end;    

%para visualizar el error
Nerr = Nt-1;

%frecuencias teoricas y FEM
figure(1)
plot(fTEO(1:Nerr),fFEM1(1:Nerr),'-o');
title('Frecuencias Teoricas vs Frecuencias FEM');
xlabel('f_{TEO} (Hz)');
ylabel('f_{FEM} (Hz)');
grid on;
box on;
grid on;

%% calculo del error

%error
errfrec(1:Nt-1,1) = 0;
for n1 = 1:Nt-1
    errfrec(n1) = 100*abs(fTEO(n1) - fFEM1(n1) )/fTEO(n1);
end;


%tabla
TablEerror = [fTEO fFEM1 errfrec];

%eror
figure(2)
plot(errfrec(1:Nerr),'-o');
title('Error Frecuencias Teoricas vs Frecuencias FEM');
xlabel('N de Modo');
ylabel('error [%]');
grid on;
box on;
grid on;


%% modos normales

pini = zeros(1,Nt-1);
PHI1 = [pini;PHI1];

%modos normales
figure(3)
subplot(2,2,1)
plot( xcorde,  PHI1(:,2)/max(PHI1(:,2)) )
title('1er Modo FEM');
grid on;
box on;
grid on;
axis([-L/2 L/2 -1 1]);

subplot(2,2,2)
plot( xcorde,  PHI1(:,3)/max(PHI1(:,3)) )
title('2do Modo FEM');
grid on;
box on;
grid on;
axis([-L/2 L/2 -1 1]);

subplot(2,2,3)
plot( xcorde,  PHI1(:,4)/max(PHI1(:,4)) )
title('3er Modo FEM');
grid on;
box on;
grid on;
axis([-L/2 L/2 -1 1]);

subplot(2,2,4)
plot( xcorde,  PHI1(:,5)/max(PHI1(:,5)) )
title('4to Modo FEM');
grid on;
box on;
grid on;
axis([-L/2 L/2 -1 1]);

%% solución forzada con fuente en x = -L/2    

% resolveremos el problema [-w^2*M1 + K1]P1 = F1
%esto consiste en resolver este sistema de ecuaciones para cada valor de 
%la frecuencia angular w = 2*pi*f, donde f = 1,2,...,1000 (Hz) es la frecuencia
%la presión sonora en x = -L/2 es P(-L/2,t) = PL*exp(j*w*t) con PL = 1 (Pa)

%armaremos los vectores de fueza con las primeras columnas de masa y
%rigidez originales, es deccir antes de colocar las condiciones de contorno
PL = 1;
Pm = M(2:N,1);
Pk = K(2:N,1);

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

%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(Nt-1,NF); %la prelocalización nos ayuda con la velocidad de ejecución

%ciclo principal de calculo
for nf = 1:NF
    
    F1 = w(nf)^2*PL*Pm - PL*Pk;     %vector de fuerza dependiente de la frecuencia
    P1(:,nf) = (-w(nf)^2*M1 + K1)\F1;   %A\b es lo mismo que inv(A)*b pero mas exacta

end;

%despues debo incorporar para todas las frecuencias la condición de
%contorno p(-L/2) = PL
Pini(1,1:NF) = PL;      %fila con valores PL para todas las frecuencias
P = [Pini;P1];          %concatemanos la fila Pini a la matriz P1

%% resultados

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

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

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

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

figure(4)
plot(x,pteo,'b',x,P(:,900),'r -o')
title('comparacion resultados teoricos vs FEM');
xlabel('x (m)');
ylabel('p (Pa)');
legend('pteo','pFEM');
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(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;

figure(6)
surf(Fmesh,Xmesh,20*log10(abs(P)/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('Lp (dB)')
axis([0,f(NF),0,L]);
grid on;
box on;
colormap jet
colorbar;