[w,f,x,y,cptim]=matumodes(ftype,a,b,nrts,nxi)

function [w,f,x,y,cptim]=matumodes(ftype,a,b,nrts,nxi)
% [wrts,I,J,f,x,y]=matumodes(ftype,a,b,nrts,nxi)

% This function evaluates frequencies and mode shapes
% for an elliptic membrane and plots the mode shapes
% ftype - 1 for modes symmetric about the x-axis
%         2 for modes anti-symmetric about the x-axis
% a,b   - major and minor semi-diameters or the ellipse 
% nrts  - vector [nrow,ncol] giving the number of
%         rows and columns in the wrts matrix
% nxi   - the number of evenly spaced points into which
%         the semi-diameter b is divided to form a
%         grid using elliptic coordinates
% wrts  - matrix of dimension nrow by ncol containing
%         frequency values found by computing roots of
%         radial mathieu functions Mc1 and Ms1
% f     - a three dimensional array for which f(:,:,k)
%         represents the k'th modal function
% x,y   - Cartesian coordinate arrays corresponding to
%         the modal coordinate matrix f
% cptim - the computation time in seconds  

%         HBW 12/01/04

cptim=clock;
if nargin<5, nxi=15; end
if nargin<4, nrts=[12,4,48]; end
nrow=nrts(1); ncol=nrts(2);
if length(nrts)==3, nmax=nrts(3);
else, nmax=nrow*ncol; end 
if nargin==0, ftype=1; a=1; b=.5; end
step=.25; ntrms=50; h=sqrt(a^2-b^2);
[xxi,eeta,x,y]=evengrid(a,b,nxi);
neta=length(eeta); nxi=length(xxi);
xxi=xxi(:)'; eeta=eeta(:);
N=inline('num2str(n)','n');
wrts=matuwr(ftype,a,b,step,nrow,ncol,ntrms);

% Sort frequencies in increasing order
[w,g]=sort(wrts(:)); q=(w*h/2).^2; 
I=(1:nrow)'*ones(1,ncol); I=I(g); 
f=zeros([size(x),nmax]);
if ftype==1, n=I-1; else, n=I; end
nf=25; ct=a/3*cos(linspace(0,4*pi,nf));
str1=['A = ',num2str(a),',  B = ',...
      num2str(b),',  TYPE = ',num2str(ftype)];  
  
% Compute modal surfaces      
for k=1:nmax
  wk=w(k); qk=q(k); ii=I(k); nk=n(k); 
  if ftype==1
    fk=ce(eeta,nk,qk)*Mc1v(xxi,nk,qk)';
  else
    fk=se(eeta,nk,qk)*Ms1v(xxi,nk,qk)';
  end
  [dum,kbig]=max(abs(fk(:))); fk=fk/fk(kbig);
  f(:,:,k)=fk; 
end
cptim=etime(clock,cptim);