%% programa cuerda/presion 02

clear all;
close all;
clc;

% ecuacion de onda
%
% d^2p    1  d^2p      dp
% ---- = ---*---- + b*----
% dx^2   c^2 dt^2      dt
% 
% condiciones de contorno
% 
% dp/dx(0,t) = 0
% 
% dp/dx(L,t) = 0
% 
% condiciones iniciales
% 
%          { A           0 <= x <= L/20
% p(x,0) = {
%          { 0        L/20 < x <= L
%          
%           
        
%% datos inciales
L = 1;          %longitud de la cuerda
A = 1;          %amplitud
c = 1;          %velocidad de propagacion de ondas
b = 0.5;
d = b/2;
x = 0:L/1000:L;  %vector de distancia

Nx = length(x); %longitud de vector x

%% ciclo principal de animacion

N = 500;          %numero de terminos de la sumatoria
kn  = zeros(N,1);
wn  = zeros(N,1);
wnd = zeros(N,1);
An  = zeros(N,1);
Bn  = zeros(N,1);


pxt(1:Nx) = 0;          %incializacion de presion

for t = 0:0.01:10

    pxt(1:Nx) = 0;      %incializacion de presion

    for n = 1:N
        
        % propiedades modales del sistema
        
        kn = n*pi/L;                %numero de onda modal
        wn = kn*c;                  %frecuencia natural modal
        beta = c^2*b/2;             %factor de amortiguamiento
        epsilon = beta/wn;          %amortiguamiento modal
        wnd = sqrt(wn^2 - beta^2);  %frecuencia natuaral amortiguada modal
       
        %coeficientes
        An = ((2*A)/(n*pi))*sin(n*pi/20);
        Cn = 0;
        Bn = (An*beta)/wnd;

        %presion sonora actualizada
        pxt(1:Nx) = pxt(1:Nx) + exp(-beta*t)*( An*cos(wnd*t) + Bn*sin(wnd*t) )*cos(kn*x);
    
    end;



    figure(1)
    plot(x,pxt)
    title('Onda de Presion Sonora');
    xlabel('x (m)');
    ylabel('p(x,t) (Pa)');
    legend(['t = ',num2str(t),' (s)']);
    axis([0 L -1.5 1.5]);
    grid on;
    box on;

    pause(0.01);

end