%% prog 011 FEM Tubo Fuente Abierto Con Montaje en el tiempo
% usar con coudado no esta refinado

clearvars;
close all;
clc;

%este es un programa de elementos finitos que representa un tubo abierto en un extremo 
%con fuete y cerrado en el otro extremo
%en ambos extremos calcularemos las frecuencias naturales (valores propios)
%y la distribucion de la presion sonora (MNV) dentro del tubo a dichas
%frecuencias (vectores propios)
%para despues determinar la respuesta forzada

%% datos inciales
Ne = 100;       %numero de elementos
NnD = 1;        %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 = 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
P0 = 2e-5;

%recordemos que este es un modelo de ondas planas
%esto se cumple cuando la longitud de onda es mucho mayor que las
%dimensiones que forman el area de la seccion transversal

%se utilizara matrices elementales

%% tiempo

fs = 88.2e3;                %frecuencia de muestreo
dt = 1/fs;                  %intervalo de muestreo
tini = 0;                   %tiempo incial 
tfin = 1;                   %tiempo final
t = tini:dt:tfin;           %vector de tiempo
Ntiempo = length(t);        %numero de muestras total de tiempo


%% 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 (el último nodo esta abierto)
nodDr = ones(N,Ndof);
nodDr(N,Ndof) = 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

%inicializacion matriz de masa global
M = zeros(Nt,Nt);

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

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

%matriz de 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 proporcional

alfa = 1e-4;
beta = 1e-4;

C = alfa*M + beta*K;
 
%% condiciones inciales

p0 =  zeros(Nt,1);     %presion inicial
dp0 = zeros(Nt,1);     %primera variación de presion inicial

%% fuerza externa

%componentes de incializacion
f = zeros(Nt,Ntiempo);

tf1 = 0.001;            %tiempo1
tf2 = 0.010;            %tiempo2
nf1 = round(tf1/dt);    %muestra1
nf2 = round(tf2/dt);    %muestra2

f(1,nf1:nf2) = 2e-4;

%% desplazamiento

%inicializamos el desplazamiento
P = zeros(Nt,Ntiempo);

%paso 1 determinar la segunda variación de presion inicial

f0 = f(:,1);
d2p0 = M\(f0 - p0 - dp0);

%paso 2 determinar la presion p1 ubicado en la segunda columna de P
P(:,2) = p0 + dp0*dt + (1/2)*d2p0*dt^2;

%paso 3 determinamos el resultado mediante recurrencia
P(:,1) = p0;


C1 = M/(dt^2) + C/(2*dt);
C2 = K - (2*M)/(dt^2);
C3 = M/(dt^2) - C/(2*dt);

for n = 2:Ntiempo-1
    P(:,n+1) = C1\( f(:,n) - C2*P(:,n) - C3*P(:,n-1) );
end;


sound(P(90,:),fs);