function out=wave2(FF,ij1,ij2)

% WAVE2(FF) or: WAVE2(FF,[i1,j1],[i2,j2]) or WAVE2() or WAVE2(n).

% Animations of solution of 2dim wave equation determined
% by matrix FF of Fourier coefficients FF. If additional 
% parameters [i1,j1],[i2,j2] are provided then in separate
% windows the corresponding eigen-modes are displayed.
%

% FF is typically a sparse small matrix, e.g. size 4 x 4
%     whose (i,j) entry is the amplitude corresponding to
%     the (i,j)-th eigen mode / pair of eigenvalues.
%
% If called with a scalar, a default example will be displayed.
%

% The more interesting cases are nonperiodic vibrations
% corresponding to sums of eigen modes such that
% sqrt((i1^2+j1^2)/(i2^2+j2^2)) is not rational.
%

% All rights reserved: Matthias Kawski, October 13, 2000. 

% Contact: kawski@asu.edu, http://math.la.asu.edu/~kawski

%

% Updates: None yet.


%
% Part 1: Interpreting the input parameters.  

%

nargs=nargin;
if nargin < 1

   % set default
   N=3;
   FF=zeros(N,N);
   FF(1,2)=5;
   FF(3,1)=3;
   i1=1;
   j1=2;
   i2=3;
   j2=1;
   nargs=3;

   repeats = 5;
elseif max(size(FF))==1
      FF=zeros(7,7);
      FF(1,1)=1;
      FF(3,3)=1/9;
      FF(5,5)=1/25;
      FF(7,7)=1/49;
      N=7;
      i1=1;
      j1=1;
      i2=3;
      j2=3;
      nargs=3;
else
   N=min(size(FF))
   i1=ij1(1);

   j1=ij1(2);

   i2=ij2(1);

   j2=ij2(2);

end

%

% Part 2: Basic settings. (x,t)-grid used for plotting etc.
%         Evaluate each eigen mode on this grid.

%

numX=36;
numY=36;
xx=[0:1/numX:1];
yy=[0:1/numY:1];
[xxx,yyy]=meshgrid(xx,yy);

for i=[1:N]
   for j=[1:N]
      if FF(i,j) ~= 0
         u(i,j,:,:)=sin(pi*i*xxx).*sin(pi*j*yyy);
      end
   end
end
%

% Position/size of 1st plot window: left bottom width height

%

scrsz = get(0,'ScreenSize');

p0=figure('Position',[35 82+scrsz(4)*2/5 scrsz(3)/3-5 scrsz(4)*2/5]);

set(p0,'DoubleBuffer','on')

umax=sum(sum(abs(FF)));
if umax == 0 
   umax = 1;
end
axis([0 1 0 1 -umax umax])
%
h=hot;
myhot=h(10:44,:);
colormap(myhot)

shading interp

%%
if nargs > 1
   p1=figure('Position',[35 35 scrsz(3)/3-5 scrsz(4)*2/5]);

   set(p1,'DoubleBuffer','on')

   axis([0 1 0 1 -umax umax])
   colormap(myhot)

   shading interp
   end

if nargs > 2
   p2=figure('Position',[40+scrsz(3)/3 35 scrsz(3)/3-5 scrsz(4)*2/5]);

   set(p2,'DoubleBuffer','on')

   axis([0 1 0 1 -umax umax])
   colormap(myhot)

   shading interp
   end

%
% Part 3: Calculate and display the solution with time
%
Tscale=0.03;
for t=[1:1000]
   ut=zeros(size(xxx));
   for i=[1:N]
      for j=[1:N]
         if FF(i,j) ~= 0
            ut=ut+FF(i,j)*sin(pi*sqrt(i^2+j^2)*t*Tscale)*squeeze(u(i,j,:,:));
         end
      end
   end
   figure(p0)
   surf(xxx,yyy,ut);
   axis([0 1 0 1 -umax umax]);
   if nargs > 1
      ut1=FF(i1,j1)*sin(pi*sqrt(i1^2+j1^2)*t*Tscale)*squeeze(u(i1,j1,:,:));
      figure(p1)
      surf(xxx,yyy,ut1);
      axis([0 1 0 1 -umax umax]);
   end
   if nargs > 2
      ut2=FF(i2,j2)*sin(pi*sqrt(i2^2+j2^2)*t*Tscale)*squeeze(u(i2,j2,:,:));
      figure(p2)
      surf(xxx,yyy,ut2);
      axis([0 1 0 1 -umax umax]);
   end   
   pause(0.05)
end
if nargs > 2
   text(0,1,umax*4/5,'Strike any key to kill all figures and terminate script')

end
if nargs > 1
   figure(p1)

   text(0,1,umax*4/5,'Strike any key to kill all figures and terminate script')

end;
figure(p0)

text(0,1,umax*3/5,'Strike any key to kill all figures and terminate script')

pause

delete(p0)

if nargs > 1
   delete(p1)

end
if nargs > 2
   delete(p2)

end
