showmo(a,b,typ,frqs,x,y,modes,tpause)

function showmo(a,b,typ,frqs,x,y,modes,tpause)
% showmo(a,b,typ,frqs,x,y,modes,tpause)
% ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
% This function makes surface or contour
% plots of the modal surfaces of a general 
% membrane for various frequencies. The 
% interactive part of the data input asks for
% a vector of modal indices (such as 1:5).
% Positive indices give modal animations.
% Negative indices give contour plots. For
% example, -2:3 would yield contour plots of
% modes 2 and 1, as well as animated plots
% of modes 1 thru 3. Also, vectors such as
% [1,-1,2,-2] are allowed.
%
% a,b    - ellipse semi diameters
% typ    - type of boundary condition.
%          1, 2, 3 or 4
% frqs   - vector of sorted frequencies
% x,y    - arrays of points defining the
%          curvilinear coordinate grid
% modes  - array of modal surfaces for 
%          the corresponding frequencies
% tpause - if zeros is input for this value,
%          the computer pauses to wait for
%          the return key to be pressed.
%          The default is tpause=.01 

%          Howard Wilson, 12/04/04

if nargin<8, tpause=0.01; end
xmin=-a; xmax=a; ymin=-b; ymax=b;
zmax=max(abs(modes(:)));
th=linspace(0,2*pi,201)';
xx=a*cos(th); yy=b*sin(th);
nf=25; ft=cos(linspace(0,4*pi,nf));
nfrqs=size(modes,3); ne=num2str(nfrqs);
scalz=0.2/zmax*max(xmax-xmin,ymax-ymin);
range=[xmin,xmax,ymin,ymax,-scalz,scalz];
disp(' ')
disp('   MODAL PLOTS FOR AN ELLIPTIC MEMBRANE')
disp(' ')
disp(['The highest allowable frequency ',...,
      'index is ',ne])
if typ==1, styp='EVEN'; else, styp='ODD'; end  
str1=['FOR ',styp,' MODE ']; str2=',  OMEGA = ';
str3=[',  A = ',num2str(a),',  B = ',...
      num2str(b)];  
while 1
   jlim=[];
   while isempty(jlim), disp(' ')
     disp(['Give a vector of mode ',...
           'indices (such as 1:2:10 or']);
     jlim=input('input 0 to stop) > ? ');
   end
   jj=find(abs(jlim)>nfrqs); 
   if length(jj)>0
     disp(['The largest allowable ',...
           'modal index is ',ne])
     pause(2), jlim(jj)=[];
   end
   if length(jlim)==1 & jlim==0
     disp(' '), disp('All done'), return
   end
   nlim=length(jlim);
   for k=1:nlim
     j=jlim(k); pltype=sign(j); j=abs(j);
     if pltype==0, continue, end
     h=.3*a; u=modes(:,:,j);
     [dumy,k]=max(abs(u(:)));
     u=h/u(k)*u; sj=num2str(j);
     sfrq=num2str(frqs(j),6);
     str=[str1,sj,str2,sfrq,str3];
     if pltype==-1
       % Draw contours  
       %u=modes(:,:,j);
       contour(x,y,u,60);
       axis equal, colormap('default')       
       hold on; plot(xx,yy,'k'); title(str)
       axis off, hold off, shg, pause(4), close
     else % pltype==1
       % Draw animated surface 
       for kk=1:nf
         surf(x,y,ft(kk)*u), axis equal
         axis([xmin,xmax,ymin,ymax,-h,h]);
         xlabel('x axis'), ylabel('y axis')
         zlabel('u(x,y)'), axis off, title(str)
         %colormap([127/255 1 212/255])
         colormap([1 1 0]), shg, pause(.05)
         %colormap('default'), drawnow, shg
       end
       pause(2)      
     end      
     if tpause==0, pause; else, pause(tpause); end
  end
end