% Math 3150-1
%
% Sample code for example 4.3.2
%
% Circular membrane with radius 1, wave velocity 1 and
%  initial shape    u  (r,theta, 0) = f(r) = (1 - r^2)*r*sin(theta)
%                     (in polar coordinates)
%  initial velocity u_t(r,theta, 0) = 0
c = 1; % wave velocity
nT = 30; % number of points for angular discretization
nr = 15; % number of points for radial  discretization
N = 10;
T = (0:nT)*2*pi/nT;
r = linspace(0,1,nr);

% Time steps and final time
nt = 40; Tmax=6;
ts = linspace(0,Tmax,nt);

% whether to take snapshots for several times
take_snapshots=0;
nsnaps = 8;
snapbasename = 'e4_2_2';
snapcount = 0;
if (mod(nt,nsnaps)~=0)
 error('number of snapshots must divide number of time steps');
end;

%%%%% Code starts here

% load Bessel J_n(r) zeros, leaves matrix alpha
% where alpha(n+1,m) is the mth zero of J_n(r)
load bz
% setup grid to plot the function
[rr,TT] = meshgrid(r,T); xx=rr.*cos(TT); yy=rr.*sin(TT);
for it=1:length(ts),
 t = ts(it);
 uu = zeros(size(rr));
 for n=1:N,
  b1n = 16/(alpha(n,2)^3 * besselj(2,alpha(n,2)));
  lambda = alpha(n,2);
  uu = uu  + b1n*cos(alpha(n,2)*t)*sin(TT).*besselj(2,alpha(n,2)*rr);
 end;
 % the rest of the loop handles plotting
 mesh(xx,yy,uu,'edgecolor','k','facecolor','none');
 axis([-1 1 -1 1 -0.4 0.4]);
 title(sprintf('t = %d',t)); xlabel('x'); ylabel('y');
 pause(0.2); % 0.2 is the time we pause
 % take a snapshot if necessary
 if (mod(it-1,nt/nsnaps)==0 & take_snapshots)
  filename = sprintf('%s_%03d.png',snapbasename,snapcount);
  fprintf([ 'saving to file ' filename '\n']);
  print('-dpng','-r50',filename);
  snapcount=snapcount+1;
 end;
end;

