[wrts,funvals]=matuwr(ftype,a,b,...

function [wrts,funvals]=matuwr(ftype,a,b,...
                         step,nrow,ncol,ntrms)
% [wrts,funvals]=matuwr(ftype,a,b,step,...
%                             nrow,ncol,ntrms)

% This function computes an array of w values 
% which are approximate frequencies related to
% the values of q giving zeros of the modified
% Mathieu functions Mc1(z,n,q), Ms1(z,n,q),
% Mc1p(z,n,q), or Ms1p(z,n,q).
% ftype   -  1 ,  2 ,  3  ,  4  for 
%           Mc1, Ms1, Mc1p, Ms1p, respectively
% a,b     - ellipse semi-diameters
% nrow    - number of rows in the root matrix
% ncol    - number of columns in the root matrix
% ntrms   - number of terms used in the series
%           expansions for Mc and Ms 
% wrts    - a rectangular array of frequencies. 
%           For function Mc1, successive rows
%           correspond to orders 0,1,...,nrows-1.
%           wrts(i,j) gives the j'th frequency 
%           relating to the i'th function.
%           For function Ms1, successive rows 
%           correspond to orders 1,...,nrow.
% funvals - values of the function at the
%           corresponding roots estimates

%           HBW 12/01/4
if nargin<7, ntrms=50; end
if nargin==0
  ftype=1; a=cosh(2); b=sinh(2); 
  nrow=7; ncol=8; ntrms=50;   
end
if ftype<3, c=.25; else, c=.1; end
wmin=c*asymtroe(ftype,a,b,0:nrow-1,1); 
wmax=2*asymtroe(ftype,a,b,nrow-1,ncol);
tol=1e-12; h2=sqrt(a^2-b^2)/2;
z=atanh(b/a); wrts=zeros(nrow,ncol); 
for j=1:nrow
  wdj=[wmin(j),step,wmax]; 
  switch ftype
    case 1      
      rts=manyr(@funcm,wdj,ncol,...
          z,j-1,h2,1,ntrms);
          %  disp([j,rts(:)']), pause(3)
      lenj=length(rts); wrts(j,1:lenj)=rts;
      qrts=w2q(rts,a,b);
      funvals(j,1:lenj)=squeeze(Mc1v(z,j-1,qrts))';
    case 2
      rts=manyr(@funcm,wdj,ncol,...
          z,j,h2,2,ntrms);
      %  disp([j,rts(:)']), pause(3)
      lenj=length(rts); wrts(j,1:lenj)=rts;
      qrts=w2q(rts,a,b);
      funvals(j,1:lenj)=squeeze(Ms1v(z,j,qrts))';
    case 3
      rts=manyr(@funcm,wdj,ncol,...
          z,j-1,h2,3,ntrms);
      %  disp([j,rts(:)']), pause(3)
      lenj=length(rts); wrts(j,1:lenj)=rts;
      qrts=w2q(rts,a,b);
      funvals(j,1:lenj)=squeeze(Mc1pv(z,j-1,qrts))';
    case 4
      rts=manyr(@funcm,wdj,ncol,...
          z,j,h2,4,ntrms);
      %  disp([j,rts(:)']), pause(3)
      lenj=length(rts); wrts(j,1:lenj)=rts;
      qrts=w2q(rts,a,b);
      funvals(j,1:lenj)=squeeze(Ms1pv(z,j,qrts))';
  end
end