function Ma = isenMaaas(ar,subsup,k)
% isenMaar  Mach number as fcn of area ratio A/Astar for isentropic
%           compressible flow
%
% Synopsis:  Ma = isenMaaas(ar)
%            Ma = isenMaaas(ar,subsup)
%            Ma = isenMaar(ar,subsup,k)
%
% Input:  ar = A/Astar, area ratio
%         subsup = (optional) flag to select subsonic (subsup<0) or
%                  supersonic (subsup>0) branch of the Ma = f(A/Astar)
%                  curve.  Default: subsup<0, subsonic branch is chosen
%         k = specific heat ratio.  Default:  k = 1.4
%
% Output:  Ma = Mach number

if nargin<2               %  supply defaults for optional input arguments
  subsup = -1;  k = 1.4;
elseif nargin<3
  k = 1.4;
end 

% --- If ar is a vector, the loop below assumes that subsup is a
%     vector of equal length.  That allows each element of ar
%     to be specified as being on a different branch.  If subsup
%     is a scalar, then replicated it to be a vector of the same
%     size as ar.
if length(ar)>1 & length(subsup)==1,  subsup = subsup*ones(size(ar));  end
 
% --- Find A/Astar as function of Ma.  The subsonic or supersonic
%     branch is selected by specifying brackets to fzero.  For the
%     subsonic branch, fzero is told to look for a root in the
%     interval eps <= Ma <= 1.  For the supersonic branch fzero is
%     told to look in the interval 1 <= Ma <= 25.  The upper limit
%     of Ma=25 is arbitrary, but should not pose a problem for nozzle flow
Ma = zeros(size(ar));
for i=1:length(ar)    %  fzero is not vectorized
  if ar(i)==1         % --- Check for special values of ar
    Ma(i) = 1;
  elseif ar(i)>5e6
    warning('Area ratio of zero detected, Ma set to zero');
    Ma(i) = 0;
  else
    if subsup(i)<0
       % pre 5.3  Ma(i) = fzero('aasmaResidual',[5e-6  1],[],[],ar(i),k);  %  subsonic branch
       Ma(i) = fzero('aasmaResidual',[5e-6  1],optimset('Display','none'),ar(i),k);  %  subsonic branch
    else
       % pre 5.e  Ma(i) = fzero('aasmaResidual',[1  25],[],[],ar(i),k);   %  supersonic branch
       Ma(i) = fzero('aasmaResidual',[1  25],optimset('Display','none'),ar(i),k);   %  supersonic branch
    end
  end
end
