% This illustrates normal modes of vibration for a simple model of a
% "square drum," i.e., a vibrating membrane on the unit square in the
% plane.  The code can be edited just below to change the length and width 
% of the drum, i.e., to make it rectangular. The user is queried for the 
% mode(s) of interest. 
%
% The initial-boundary value problem is 
%
%    u_tt = u_xx + u_yy in D = [0,1]x[0,1], 
%    u = 1, u_t = 0 in D when t = 0, 
%    u = 0 on the boundary of D for t >= 0. 
%
% The discretization is by centered differences on a uniform mx(times)my
% grid. The following important parameters can be changed by editing:
%   Lx = domain width (x-direction)
%   Ly = domain height (y-direction)
%   mx = number of gridpoints in the x-direction
%   my = number of gridpoints in the y-direction
%
% The "First, try this" line below sets Lx = Ly = 1, mx = my = 32.
% Uncommenting the "Now, try this" line sets Lx = 3, Ly = 1, mx = 96, 
% my = 32. The user can change these to try other values as desired;
% however, making the grid much finer may undesirably increase the time
% required to compute the eigenvalues and eigenvectors. 

% First, try this.
Lx = 1; Ly = 1;

% Now, try this: 
%Lx = 3; Ly = .5;

my = 32;
mx = Lx*my;

% Discretize the laplacian. 
L = -laplacian(mx,Lx,my,Ly);

% Compute the eigenvectors and eigenvalues. Note that L is in sparse
% storage mode and must be converted to full storage mode for the eig
% routine.  
tic;[V,D] = eig(full(L)); 
fprintf('\n Time for computing eigenvalues/eigenvectors = %g\n',toc);

mu = sqrt(diag(D));

hx = Lx/(mx+1);
hy = Ly/(my+1);
x = 0:hx:Lx;
y = 0:hy:Ly;

mode = 1;

while mode ~=0 
    mode = input('Type desired mode or zero to quit:   ');

    if mode == 0, break; end 
    
    vbar = reshape(V(:,mode),mx,my);
    z = zeros(mx,1);
    vbar = [z,vbar,z];
    z = zeros(my+2,1);
    vbar = [z,vbar',z];

    vmax = max(max(abs(vbar)));

    figure(1); clf; 
    surf(x,y,vbar); 
    axis([0 Lx 0 Ly -1*vmax 1*vmax]); 
    title(['Normal mode ',num2str(mode), ' (stationary).']);

    t = 0:.01:pi;

    figure(2); clf;
 
    for j = 1:length(t),
        surf(x,y,cos(mu(mode)*t(j))*vbar); 
%        axis([0 1 0 1 -1*vmax 1*vmax]); 
        axis([0 Lx 0 Ly -1*vmax 1*vmax]); 
        title(['Normal mode ',num2str(mode), ' (in motion), t = ',num2str(t(j))]); 
        drawnow; pause(.01);
    end
end




