%% prog fem011 tubo abierto abierto autovalores con montaje 

%este programa calcula los valores propios de un tubo de 
%longitud L abierto en ambos extremos
%area de la seccion transversal variable radioA(nel) = a*exp( 2.0*xcorde(nel) );
%el numero de elementos es variable 


%se utilizara matrices elementales

clear all;
close all;
clc;

%% inicializacion

Ne = 100;       %numero de elementos
NnD = 2;        %numero de DOF s con condiciones drichelet (tubo abierto ambos extremos)
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

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

nodDr(1) = 0;   %p(-L/2,t) = 0
nodDr(N) = 0;   %p( L/2,t) = 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 - rigidez elemental y global
    
%matriz de masa elemental
Me = [2 1;1 2];
Me = (A*ax)/(3*c^2)*Me;  
    
%inicializacion matriz de masa global
M = zeros(Nt,Nt);

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

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

%inicializacion de radios
radioA = zeros(1,Ne);

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


%% problema de valores propios

[PHI,LAM] = eig(K,M);

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

%normalizacion de los vectores propios
mnorm = diag( transpose(PHI)*M*PHI );

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

%incorporar a la matriz modal la condiciones de contorno %p(-L/2,t) = p(L/2,t) = 0
%esto significa incorporar una fila de zeros al principio y otra al final
PHI = [zeros(1,Nt);PHI;zeros(1,Nt)];

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

%para visualizar el error
Nerr = Nt;

%frecuencias teoricas y FEM
figure(1)
plot(fTEO(1:Nerr),fFEM(1:Nerr),'-o');
title('frecuencias teoricas vs frecuencias FEM');
xlabel('f_{TEO} (Hz)');
ylabel('f_{FEM} (Hz)');
axis([min(fTEO),max(fTEO),min(fFEM),1.1*max(fFEM)]);
grid on;
box on;
grid on;

%% calculo del error

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


%tabla
TablEerror = [fTEO fFEM errfrec];

nn = 1:length(errfrec);

%error
figure(2)
plot(nn,errfrec(1:Nerr),'-o');
title('error frecuencias teoricas vs frecuencias FEM');
xlabel('n de Modo');
ylabel('error [%]');
axis([1,length(errfrec),min(errfrec),1.1*max(errfrec)]);
grid on;
box on;
grid on;


%% modos normales

%figuras

figure(3)
subplot(2,2,1)
plot( xcorde,  PHI(:,1)/max( abs(PHI(:,1)) ) )
title('1er Modo FEM');
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(1))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);

subplot(2,2,2)
plot( xcorde,  PHI(:,2)/max( abs(PHI(:,2)) ) )
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(2))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);

subplot(2,2,3)
plot( xcorde,  PHI(:,3)/max( abs(PHI(:,3)) ) )
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(3))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);

subplot(2,2,4)
plot( xcorde,  PHI(:,4)/max( abs(PHI(:,4)) ) )
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(4))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);



figure(4)
subplot(2,2,1)
plot( xcorde,  PHI(:,5)/max( abs(PHI(:,5)) ) )
title('1er Modo FEM');
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(5))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);

subplot(2,2,2)
plot( xcorde,  PHI(:,6)/max( abs(PHI(:,6)) ) )
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(6))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);

subplot(2,2,3)
plot( xcorde,  PHI(:,7)/max( abs(PHI(:,7)) ) )
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(7))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);

subplot(2,2,4)
plot( xcorde,  PHI(:,8)/max( abs(PHI(:,8)) ) )
title(['distribucion de presion sonora normalizada para f_{FEM} = ',num2str(fFEM(8))] );
grid on;
box on;
grid on;
axis([-L/2 L/2 -1.5 1.5]);