% solucion de la ecuacion de ondas bi-dimensional
% sobre un rectángulo
% los coeficientes a, b (en forma matricial)
% corresponden a los términos de la expansion en 
% productos de senos sen(n pi x) * sen (m pi y)
% del desplazamiento inicial (a) y de la velocidad inicial (b)
%
% hay 9 casos preprogramados que se eligen poniendo el 
% valor correspondiente a la variable 'caso'

caso = 4;
close all;
L = pi; H = pi;

switch (caso)
    case 1,
        a = 1; b = 0;
    case 2,
        a = 0; b = 1;
    case 3,
        a = [0 1]; b = [0 0];
    case 4,
        a = [0 1]'; b = [0 0]';
    case 5,
        a = [0 1; 0 0]; b = [0 0; 1 0];
    case 6,
        a = [0 0; 0 1]; b = zeros(size(a));
    case 7,
        a = [0 1 ; 0 0]; b = [0 0 ; 1 0];
    case 8,
        a = [1 0.1 ; 0.1 -0.5]; 
        b = zeros (size(a));
    otherwise,
        a = zeros(3,3); 
        b = [1 0.1 0.1 ; 0 0 0 ; 0.1 0 -0.1];
end
        
[N, M] = size(a);

x = [0:0.05:L];
y = [0:0.05:H];

[xx,yy] = meshgrid(x,y);
figure(1)
unos = ones(size(xx));

for t=0:.03:6

  u = zeros(size(xx));
  for n=1:N    
      for m=1:M
          sinx = sin(n*pi/L*xx);
          siny = sin(m*pi/H*yy);
          sinsin = sinx .* siny;
          u = u + ( a(n,m) * cos(sqrt((n*pi)^2+(m*pi)^2)*t) ...
                  + b(n,m) * sin(sqrt((n*pi)^2+(m*pi)^2)*t) )* sinsin;
      end
  end
  surf(x,y,u,unos);
  
  axis([0 L 0 H -1 1])
  if (t==0)
      pause
  else
      pause(0.01)
  end
end