%% prog fem014 tubo fuente en un extremo cerrado autovalores con montaje EXPONENCIAL

%este programa calcula los valores propios de un tubo de 
%longitud L con fuente en un extremo cerrado extremo
%area de la seccion transversal exponencial
%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
rad = sqrt(A)/pi;

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

%inicializacion
M = zeros(Nt,Nt);

%loop principal montaje
for nel = 1:Ne
    
    %area transversal exponencial
    A = pi*(rad*exp(0.01*nel))^2;
    % la matriz de masa elemental
    Me = [2 1;1 2];
    Me = (A*ax)/(3*c^2)*Me;  
    
    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

%inicializacion
K = zeros(Nt,Nt);

%loop principal montaje
for nel = 1:Ne

    %area transversal exponencial
    A = pi*(rad*exp(0.01*nel))^2;
    % la matriz de rigidez elemental
    Ke = [1 -1;-1 1];
    Ke = A/(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;
                         
                         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;    

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

%para visualizar el error
Nerr = N;

%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)');
grid on;
box on;
grid on;

%% calculo del error

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


%tabla
TablEerror = [fTEO fFEM 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

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

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

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

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