function out=wave2disk(FF,ij1,ij2)
% WAVE2DISK(FF) or: WAVE2(FF,[i1,j1],[i2,j2]) or WAVE2() or WAVE2(n).
% Animations of solution of 2dim wave equation on a
% circular (disk) domain, determined  by matrix FF 
% of Fourier-Bessel 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 
%     (j oscillations in theta ==> j-th Bessel of 2ns kind,
%      i oscillations radially -- i-th root of BesselJ(j,.)
%
%     Initially work with FF=zeros(4,4); FF(2,1)=1; etc.
% 
% If called with a scalar, a default example will be displayed. 
%
% This function uses a table of the first 10 zeros of 
% each of the Bessel functions J(0, .) to J(10, . ).
%
% All rights reserved: Matthias Kawski, October 23, 2000. 
% Contact: kawski@asu.edu, http://math.la.asu.edu/~kawski
%
% Updates: None yet.
% 
% Part 0: Table look-up -- j-th zero of (i-1)-th Bessel 
%         of 2nd kind  
%
Jjroots=...
[ 2.404825558, 5.520078110, 8.653727913, 11.79153444, ...
 14.93091771, 18.07106397, 21.21163663, 24.35247153, ...
 27.49347913, 30.63460647, 33.77582021;...
 3.831705970, 7.015586670, 10.17346814, 13.32369194, ...
 16.47063005, 19.61585851, 22.76008438, 25.90367209, ...
 29.04682854, 32.18967991, 35.33230755;...
 5.135622302, 8.417244140, 11.61984117, 14.79595178,... 
 17.95981949, 21.11699705, 24.27011231, 27.42057355,...
 30.56920450, 33.71651951, 36.86285651;...
 6.380161896, 9.761023130, 13.01520072, 16.22346616,...
 19.40941523, 22.58272959, 25.74816670, 28.90835078,...
 32.06485241, 35.21867074, 38.37047243;...
 7.588342435, 11.06470949, 14.37253667, 17.61596605,...
 20.82693296, 24.01901952, 27.19908777, 30.37100767,...
 33.53713771, 36.69900113, 39.85762730;...
 8.771483816, 12.33860420, 15.70017408, 18.98013388,...
 22.21779990, 25.43034115, 28.62661831, 31.81171672,...
 34.98878129, 38.15986856, 41.32638325;...
 9.936109524, 13.58929017, 17.00381967, 20.32078921,...
 23.58608444, 26.82015198, 30.03372239, 33.23304176,...
 36.42201967, 39.60323942, 42.77848161;...
 11.08637002, 14.82126873, 18.28758283, 21.64154102,...
 24.93492789, 28.19118846, 31.42279419, 34.63708935,...
 37.83871738, 41.03077369, 44.21540851;...
 12.22509226, 16.03777419, 19.55453643, 22.94517313,...
 26.26681464, 29.54565967, 32.79580004, 36.02561506,...
 39.24044800, 42.44388774, 45.63844418;...
 13.35430048, 17.24122038, 20.80704779, 24.23388526,...
 27.58374896, 30.88537897, 34.15437792, 37.40009998,...
 40.62855372, 43.84380142, 47.04870074;...
 14.47550069, 18.43346367, 22.04698536, 25.50945055,...
 28.88737506, 32.21185620, 35.49990921, 38.76180702,...
 42.00419024, 45.23157410, 48.44715139;...
 15.58984788, 19.61596690, 23.27585373, 26.77332255,...
 30.17906118, 33.52636408, 36.83357134, 40.11182327,...
 43.36836095, 46.60813268, 49.83465351];
% 
% Part 1: Interpreting the input parameters.  
%
nargs=nargin; 
if nargin < 1
   % set default    
   N=3;    
   FF=zeros(N,N);    
   FF(1,3)=1;    
   FF(2,1)=1;    
   i1=1;    
   j1=3;    
   i2=2;    
   j2=1;    
   nargs=3;
   repeats = 5; 
elseif max(size(FF))==1       
   FF=zeros(2,2); 
   FF(1,2)=1;  
   FF(2,1)=1;       
   N=2;    
   i1=1;    
   j1=2;   
   i2=2;    
   j2=1;    
   nargs=3; 
elseif nargin == 1
   % set default    
   N=min(size(FF))  
   nargs=1;
   repeats = 5; 
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.
%
numR=24; 
numT=36; 
rr=[0:1/numR:1]; 
th=[0:2*pi/numT:2*pi]; 
[rrr,tth]=meshgrid(rr,th);
xxx=rrr.*cos(tth);
yyy=rrr.*sin(tth);

for i=[1:N]    
   for j=[1:N]     
      if FF(i,j) ~= 0         
         u(i,j,:,:)=FF(i,j)*besselj(j,Jjroots(j+1,i)*rrr).*sin(j*tth);
         uumax(i,j)=max(max(abs(squeeze(u(i,j,:,:)))))
      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(uumax)); 
if umax == 0     
   umax = 1; 
end 
axis([-1 1 -1 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([-1 1 -1 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([-1 1 -1 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)*squeeze(u(i,j,:,:))*...
               sin(Jjroots(i+1,j)*t*Tscale); 
         end       
      end    
   end    
   figure(p0)    
   surf(xxx,yyy,ut);    
   axis([-1 1 -1 1 -umax umax]);    
   shading interp
   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([-1 1 -1 1 -umax umax]); 
      shading interp    
   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([-1 1 -1 1 -umax umax]);    
      shading interp    
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