%% 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
% 
% p(0,t) = 0
% 
% dp/dx(L,t) = 0
% 
% condiciones iniciales
% 
%          { 1           0 <= x <= L/100
% p(x,0) = {
%          { 0        L/100 < x <= L
%          
%           
%              {  0             0 <= x <= 98*L/100
% dp/dt(x,0) = {  100     98*L/100 < x <= 99*L/100     
%              { -100     99*L/100 < x <= L
         
         
%% datos inciales
L = 1;          %longitud de la cuerda
A = 1;          %amplitud
c = 1;          %velocidad de propagacion de ondas
b = 1;
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:5

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

    for n = 1:N
        
        % propiedades modales del sistema
        
        kn = (2*n-1)*pi/(2*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 = wn*sqrt(1 - epsilon^2);	   %frecuencia natuaral amortiguada modal
       
        %coeficientes
        An = (2-2*cos(kn*L/100))/(kn*L);
        Cn = 100*(-4*cos(99*kn*L/100) + 2*cos(98*kn*L/100) + 2*cos(kn*L))/(kn*L);
        Bn = (Cn + An*beta)/wnd;

        %presion sonora actualizada
        pxt(1:Nx) = pxt(1:Nx) + exp(-beta*t)*( An*cos(wnd*t) + Bn*sin(wnd*t) )*sin(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 1]);
    grid on;
    box on;

    pause(0.01);

end