%% membrana

clear all; 
close all; 
clc;

%% parametros inciales

%dimensiones membrana
Lx = 1.00;
Ly = 2.00;
dx = Lx/200;
dy = Ly/400;
A = 100.0;

%membrana
[x,y] = meshgrid(0:dx:Lx,0:dy:Ly);
Ndimx = length(0:dx:Lx);
Ndimy = length(0:dy:Ly);

%velocidad del sonido y amortiguamiento membrana
c = 1.0;
b = 5e-1;
dnxny = c^2*b/2;

%% calculo principal

Nx = 200;
Ny = 200;

tini = 0.00;
dt   = 0.01;
tend = 10.00;  

for t = tini:dt:tend
    
    p(1:Ndimy,1:Ndimx) = 0;
    
    for nx = 1:Nx
        
        for ny = 1:Ny
            
            %numeros de onda 
            knx = nx*pi/Lx;
            kny = ny*pi/Ly;
            
            %frecuencias naturales - frecuencias naturales amortiguadas
            wnxny = c*sqrt( knx^2 + kny^2 );
            wdnxny = sqrt( wnxny^2 - dnxny^2 );
            
            %cooefcientes Anxny
            Anxnx = 0;
            
            %cooefcientes Bnxny
            Bnxnx = ( 4*A )/( nx*ny*wdnxny*pi^2 )*( sin(25*nx*pi/50) - sin(24*nx*pi/50) )*( sin(25*ny*pi/50) - sin(24*ny*pi/50) );
            
            %desplazamiento membrana
            pt = exp( -dnxny*t )*( Anxnx*cos( wdnxny*t ) + Bnxnx*sin( wdnxny*t ) );
            px = sin( knx*x );
            py = sin( kny*y );
            
            p = p + pt*px.*py;
            
        end
        
    end
    
    figure(1)
    surf(x,y,p,'FaceColor','interp',...
                'EdgeColor','none',...
                'FaceLighting','phong');
    title('membrana')
    xlabel('x (m)');
    ylabel('y (m)');
    zlabel('p(x,y,t) (N/m^2)');
    legend(['t = ',num2str(t),' (s)']);
    axis([0,Lx,0,Ly,-10.0e0,10.0e0]);
    daspect([1,1,1])
    view(0,90);
    colormap jet;
    colorbar;
    grid on;
    box on;
            
    pause(dt);
            
    
end