[frqs,modes,x,y,cptim]=runelip(type,a,b,...

function [frqs,modes,x,y,cptim]=runelip(type,a,b,...
                                         nrts,noplot)
% [frqs,modes,x,y,cptim]=runelip(type,a,b,nrts,noplot)
% 
% This program uses separation of variables, elliptical
% coordinates, and Mathieu functions to compute fre-
% quencies and mode shapes for an elliptical membrane
% having zero deflection on the boundary. Modes that
% are symmetric or anti-symmetric about the x-axis are
% allowed. The pertinent boundary value problem is:
% (Dxx+Dyy)u(x,y)+(w^2)u(x,y)=0 for (x/a)^2+(y/b)^2<1
%                      u(x,y)=0 for (x/a)^2+(y/b)^2=1
% If no input data is passed through the argument 
% list, data is read interactively. Output provides
% animated modal plots or contour plots for the 
% various frequencies arranged in increasing order.
%
% type      - use 1 for even modes symmetric about the
%             x-axis. Use 2 for odd modes anti-
%             symmetric about the x-axis.        
% a,b       - the ellipse major and minor semi-
%             diameters along the x and y axes
% nrts      - a vector of the form [nrow,ncol] where
%             nrow is the number of Mathieu functions
%             used and ncols is the number of roots
%             computed for each function. Hence, the
%             number of frequencies computed equals
%             nrow*ncol
% noplot    - give a value for this parameter only
%             when no modal plots are desired. 
%             Otherwise, omit parameter noplot.
%
% frqs      - a vector of natural frequencies
%             arranged in increasing order.
% modes     - a three dimensional array in which
%             modes(:,:,j) defines the modal 
%             deflection surface for frequency
%             frqs(j).
% x,y       - curvilinear coordinate arrays of
%             points in the membrane where modal
%             function values are computed. 
% cptim     - the cpu time in seconds used to 
%             solve for the frequencies and compute
%             the modal surfaces.
%
% User m functions called: 
%             matumodes showmo ce se Mc1 Ms1
%             cev sev Mc1v Ms1v matue evengrid
%             matuwr w2q manyr funcm asymtroe
%             bes funcm

%             by Howard Wilson, 1/2/04
%                [e-mail]
close;
if nargin==0
   disp(' ')
   disp('VIBRATION FREQUENCIES AND MODE SHAPES')
   disp('       OF AN ELLIPTIC MEMBRANE       ')
   disp(' ')          
   v=input(['Input the major and minor ',...
         'semi-diameters > ? '],'s');
   v=eval(['[',v,']']); a=v(1); b=v(2); disp(' ')
   disp('Select the modal symmetry option')
   type=input('1<=>even, 2<=>odd > ? ');
   disp(' ')
   disp('Input the number of Mathieu functions (nrow) and')
   disp('the number of roots for each function (ncol)')
   nrts=input('(typical values are 12,4) > ? ','s');  
   nrts=eval(['[',nrts,']']);
   disp(' ')
   disp('Computing frequencies and modes takes awhile.')
   disp('              PLEASE WAIT.')
end

nxi=20; % This value determines the grid resolution
[frqs,modes,x,y,cptim]=matumodes(type,a,b,nrts,nxi);
indx=ones(length(frqs),1)*type;
if nargin==5, return, end

% Plot a sequence of modal functions
neig=length(frqs);
disp(' '), disp(['Computation time  = ',...
      num2str(sum(cptim)),' seconds.'])
disp(['Number of modes   = ',num2str(neig)]);
disp(['Highest frequency = ',...
      num2str(frqs(end))]), disp(' ')
disp('Press return to see modal plots.')
pause, tpause=.01;
showmo(a,b,type,frqs,x,y,modes,tpause);