% solucion de la ecuacion de ondas bi-dimensional
% en un circulo
% los coeficientes a0, b0 corresponden a la parte 
% radialmente simetrica del desplazamiento y 
% velocidad inicial, respectivamente
%
% los coeficientes a1, a2 (en forma matricial)
% corresponden a los términos de la expansion en 
% cosenos y senos respectivamente del desplazamiento inicial
%
% los coeficientes b1, b2 (en forma matricial)
% corresponden a los términos de la expansion en 
% cosenos y senos respectivamente de la velocidad inicial
%
% hay 8 casos pre-programados, que se eligen poniendo el 
% valor correspondiente a la variable 'caso'

caso = 8 ;
close all
R = 1;
c = 1;
switch (caso)
    case 1,
        a0 = 1; b0 = 0;
        a1 = []; a2 = []; b1 = []; b2 = [];
    case 2,
        a0 = 0; b0 = 1;
        a1 = []; a2 = []; b1 = []; b2 = [];
    case 3,
        a0 = [0 0 .3 .1 .05 .02];
        b0 = zeros(size(a0));
        a1 = []; a2 = []; b1 = []; b2 = [];
    case 4,
        a0 = 0; b0 = 0;
        a1 = 1; a2 = 0;
        b1 = 0; b2 = 1;
    case 5,
        a0 = 0; b0 = 0;
        a1 = 1; a2 = 0;
        b1 = 0; b2 = 0;
    case 6,
        a0 = 0; b0 = 0;
        a1 = 0; a2 = 0;
        b1 = 0; b2 = 1;
    case 7,
        a0 = 0; b0 = 0;
        a1 = [0 1 ; 0 0]; a2 = [0 0 ; 0 0];
        b1 = [0 0 ; 0 0]; b2 = [0 1 ; 0 0];
    otherwise,
        a0 = [0 0 0 0 0 .3 .1 .05 .02 0.01 0.005 0.002];
        b0 = zeros(size(a0));
        a1 = []; a2 = []; b1 = []; b2 = [];
        c = 0.2;
end

figure(1)

M0 = length(a0);
[N, M] = size(a1);

z0 =[ 
2.404825557695773 5.520078110286311 8.653727912911011 11.79153443901428 ...
14.93091770848779 18.07106396791092 21.21163662987926 24.35247153074930 ...
27.49347913204025 30.63460646843198 33.77582021357357 36.91709835366404 ...
40.05842576462824 43.19979171317673 46.34118837166181 49.48260989739782
];
z = [
3.831705970207512 7.015586669815613 10.17346813506272 13.32369193631422 16.47063005087763 19.61585851046824
5.135622301840683 8.417244140399855 11.61984117214906 14.79595178235126 17.95981949498783 21.11699705302185
6.380161895923984 9.761023129981667 13.01520072169843 16.22346616031877 19.40941522643501 22.58272959310444
7.588342434503804 11.06470948850119 14.37253667161759 17.61596604980483 20.82693295696239 24.01901952477111
];

lambda0 = (z0/R).^2;
lambda = (z/R).^2;

r = [0:0.05:R];
th = [-pi:0.05:pi+0.05];

[rr,thth] = meshgrid(r,th);
figure(1)

rcos = rr.*cos(thth);
rsen = rr.*sin(thth);
unos = ones(size(rr));

for t=0:.01:6*pi/lambda0(1)

  u = zeros(size(rr));
  if (M0 > 0)
      for m=1:M0
          
          besselr = besselj(0,sqrt(lambda0(m))*rr);
          
          u = u + ( a0(m) * cos(c*lambda0(m)*t) ...
              + b0(m) * sin(c*lambda0(m)*t) )* besselr;
      end
  end
  for n=1:N
      for m=1:M
          besselr = besselj(n,sqrt(lambda(n,m))*rr);
          cosntheta = cos(n*thth);
          senntheta = sin(n*thth);
          phi1 = besselr.*cosntheta; 
          phi2 = besselr.*senntheta;
          coslambdat = cos(c*lambda(n,m)*t);
          senlambdat = sin(c*lambda(n,m)*t);
          u = u + ( a1(n,m) * coslambdat + b1(n,m) *senlambdat ) * phi1 ...
              + (a2(n,m) * coslambdat + b2(n,m) *senlambdat ) * phi2 ;
      end
  end
  surf(rcos,rsen,u,unos);
  
  axis([-R R -R R -1 1])
  if (t==0)
      pause
  else
      pause(0.01)
  end
end