% Math 3150-1
% 
% Sample code for example 3.7.1: 
%
% Square membrane with side 1, wave velocity c=1/pi and
%     initial shape     u  (x,y,0) = f(x,y) = x*(1-x)*y*(1-y)
%     initial velocity  u_t(x,y,0) = g(x,y) = 0

% number of terms in the expansion
N=20; M=20;

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

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

%%%% Code starts here

% setup a grid to plot the function
Nx=20; Ny=20; % number of points in x and y directions
x=linspace(0,1,Nx); y=linspace(0,1,Ny); 
[xx,yy]=meshgrid(x,y);

for it=1:length(ts),
 t = ts(it);
 ss = zeros(size(xx));
 for n=1:N,
  for m=1:M,
    % Sine Series coefficients for initial shape
    Bmn = 16*((-1)^n-1)*((-1)^m-1)/(pi*n)^3/(pi*m)^3;
    % (essentially) Sine Series coefficients for initial velocity
    Bmnstar = 0;
    lambdamn = sqrt(m^2+n^2);
    ss = ss + sin(m*pi*xx).*sin(n*pi*yy) * ...
             (Bmn*cos(lambdamn*t) + Bmnstar*sin(lambdamn*t));
   end;
  end;
  % the rest of the loop handles plotting
  mesh(xx,yy,ss,'edgecolor','k','facecolor','none');
  axis([0 1 0 1 -1/16 1/16]); 
  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;
