%% ProgFEM012TuboFuenteCerradoConMontajeDiffCentrTiempo

%resultado

clearvars;
close all;
clc;

%% datos iniciales

%tiempo
ti = 0;         %tiempo incial (s)
tf = 1;         %tiempo final  (s)
fs = 2*44100;   %frecuencia de muestreo (Hz)
dt = 1/fs;      %intervalo de tiempo (s)
t = ti:dt:tf;   %vector de tiempo (s)

% dimension tiempo
NTime = length(t); 
Nm = NTime-1; 

%% inicializacion matrices y vectores

Ne = 100;           %numero de elementos
NnD = 0;            %numero de DOF s con condiciones drichelet
N = Ne+1;           %numero de DOFs - funciones de interpolacion
NTDOF = Ne+1-NnD;   %numero de DOFs Totales- 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.0e0;          %longitud tubo m
ax = L/(2*Ne);      %logitus de 1/2 elemento
rad = 0.01;         %radio
A = pi*rad^2;       %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 volumetrica adiabatica

%% matriz de cordenadas

%inicio
Xcord = zeros(N,3);

%coordenadas
deltax = L/Ne;
xcorde = transpose( (0:N-1)*deltax  );

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
nodDr = ones(N,Ndof);   %inicializacion
nodDr(N) = 1;           %cerrado

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    

%% matrices de masa y rigidez

%inicializacion
M = zeros(NTDOF,NTDOF);
K = zeros(NTDOF,NTDOF);

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

    %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;
                        
                         M(NN,MM) = M(NN,MM) + Me( fila, colu );
                         K(NN,MM) = K(NN,MM) + Ke( fila, colu );
                         
                     end
                     
                 end
             end
        end 
    end    
end    
    
%se saca primera fila/columna por fuente
M1 = M(2:NTDOF,2:NTDOF);
K1 = K(2:NTDOF,2:NTDOF);

%% matriz de amortiguamiento proporcional

%coeficientes de proporcionalidad
alfa = 1.0e-5;
beta = 1.0e-3;

%matriz de amortiguamiento porporcional
C = alfa*M + beta*K;

%saca primera fila/columna por fuente
C1 = C(2:NTDOF,2:NTDOF);

%% vector de fuerza por reaccion

%parte vectores de fuerza
F1M = M(2:NTDOF,1);         %parte 1 primera columna que se convierte en fuerza de reaccion inercial
F1K = K(2:NTDOF,1);         %parte 2 primera columna que se convierte en fuerza de reaccion elástica
F1C = C(2:NTDOF,1);         %parte 3 primera columna que se convierte en fuerza de reaccion disipativa

%inicializacion de la presion compleja total
P = zeros(NTDOF-1,NTime);

%condiciones inciales de presion y derivada temporal de la presion
P0    = zeros(NTDOF-1,1);
P0dot = zeros(NTDOF-1,1);

%% fuerza un componente
tf1 = 0.001;               
tf2 = 0.005;
Nf1 = round(tf1/dt);
Nf2 = round(tf2/dt);

F0 = 1e-2;
f(1:NTime) = 0;          %fuerza inicializada con ceros

%pulso rectangular
f(Nf1:Nf2-1) = F0;

%pulso senoide
% f(Nf1:Nf2-1) = F0*sin( 2*pi*t(1:Nf2-Nf1)/( t(Nf2)-t(Nf1) ) );

figure(1)
plot(t,f,'k');
title('fuerza');
xlabel('tiempo t (s)');
ylabel('presion sonora p(t) (Pa)');
legend('f(t)');
axis([t(1) t(NTime) -1.2*max(abs(f)) 1.2*max(abs(f))]);
grid on;
box on;

%% proceso principal

% fuerza inicial t = 0
F = zeros(NTDOF-1,NTime);%N-2
F(1,:) = f;

% aceleracion inicial t = 0
P0dotdot = (M1^-1)*( F(:,1) - C1*P0dot - K1*P0 );

% posicion inicial t = dt
P(:,2) = P0 + P0dot*dt + 1/2*P0dotdot*dt^2;

% coeficientes 
cof1 = M1 / (dt^2) + C1 / (2*dt);   %matriz
cof2 = 2 * M1 /(dt^2) - K1;         %matriz 
cof3 = M1 / (dt^2) - C1 / (2*dt);   %matrizmatriz M(3x3)

% loop para los demas puntos
for tn = 2:NTime 
    F(1) = f(tn);
    P(:,tn+1) = (cof1^-1) * ( F(:,tn) + cof2* P(:,tn) - cof3* P(:,tn-1) );      
end

Pini = zeros(1,NTime+1);
%Pfin = zeros(1,NTime+1);

P = [Pini;P];
%% graficos en el dominio del tiempo

p_1 = P(2,1:NTime);            %presion sonora nodo incial + 1
p_N = P(Ne,1:NTime);           %presion sonora nodo final - 1

% superposición soluciones Pn (n = 1:Ne-1)
figure(2)
plot(t,p_1,'b');
title('solucion numerica respecto al tiempo primer nodo');
xlabel('tiempo t (s)');
ylabel('presion sonora p(t) (Pa)');
legend('p_1(t)');
grid on;
box on;

figure(3)
plot(t,p_N,'r');
title('solucion numerica respecto al tiempo ultimo nodo');
xlabel('tiempo t (s)');
ylabel('presion sonora p(t) (Pa)');
legend('p_N(t)');
grid on;
box on;

figure(4)
plot(t,p_1,'b');
hold on; 
plot(t,p_N,'r');
hold on;
plot(t,f,'k');
hold off;
legend('p_1(t)','p_N(t)','f(t)');
title('fuerza y solucion numerica respecto al tiempo primer y ultimo nodo');
xlabel('tiempo t (s)');
ylabel('presion sonora p(t) (Pa)');
grid on;
box on;

%% analisis de frecuencia

%frecuencias
dfrec = fs/NTime;      %intervalo de muestreo - frecuencia
frec = 0:dfrec:fs-dfrec; %frecuencia

%fft
P_1 = fft(p_1)/NTime;
P_N = fft(p_N)/NTime;
F_N = fft(f)/NTime;

%% grafico en el dominio de la frecuencia

figure(5)
semilogy(frec,abs(F_N),'k');
hold off
title('espectro unilateral fuerza');
xlabel('f (Hz)');
ylabel('P(f) (Pa)');
legend('F(f)');
grid on;
box on;
axis([0, 2500, 10e-6 1.1*max(abs(P_1))]);

figure(6)
semilogy(frec,abs(P_1),'b');
title('espectro unilateral presion sonora primer nodo');
xlabel('f (Hz)');
ylabel('P(f) (Pa)');
legend('P_1(f)');
grid on;
box on;
axis([0, 2500, 10e-6 1.1*max(abs(P_1))]);

figure(7)
semilogy(frec,abs(P_N),'r');
title('espectro unilateral presion sonora ultimo nodo');
xlabel('f (Hz)');
ylabel('P(f) (Pa)');
legend('P_N(f)');
grid on;
box on;
axis([0, 2500, 10e-6 1.1*max(abs(P_1))]);

figure(8)
semilogy(frec,abs(P_1),'b');
hold on
semilogy(frec,abs(P_N),'r');
hold on
semilogy(frec,abs(F_N),'k');
hold off
title('espectro unilateral - escala logaritmica');
xlabel('f (Hz)');
ylabel('P(f) (Pa)');
legend('P_1(f)','P_N(f)','F(f)');
grid on;
box on;
axis([0, 2500, 10e-6 1.1*max(abs(P_1))]);


%% animacion temporal

x = xcorde;
for tn = 1:round(NTime/50)
    
    figure(9)
    plot(x,P(:,tn),'b');
    title('solución numerica respecto al tiempo todos los nodos');
    xlabel('espacio x (m)');
    ylabel('presion sonora p(x,t) (Pa)');
    legend('p(x,t)');
    axis([0 L -5.0e-1 5.0e-1]);
    grid on;
    box on;
    pause(0.001);

end