function m = pmm (nu,gamma)
%pmm         This function computes the mach
%            no. that coresponds to a given Prandtl-Meyer
%            angle & gamma.  The form is:
%            [mach no.] = pmm (nu, gamma)
%            (note: nu and gamma should be scalars,
%            and nu should be given in degrees.)


nu=nu*pi/180;              %convert to radians

%set up parameters for bisection method.

mmax = 10;                 %max subsonic mach no.
mmin = 1;
epsilon = .001;            % convergence limit
converg = 0;               % this is a dummy variable

while ~converg

     rtbis = ((mmax-mmin)/2)+mmin;   %set bisection value

     expn = atan((rtbis^2-1)^.5);
     rat=(gamma+1)/(gamma-1);
     g = (rat^.5)*atan((((gamma-1)/(gamma+1))*...
          (rtbis^2-1))^.5)-expn;
     
     %define the convergence function

     func = g-nu; 

     if (func < 0)
          mmin = rtbis;     %compare, define new limits
     else
          mmax = rtbis;
     end

     if abs(func) <= epsilon   %compare epsilon
          converg = 1;
     end
end
m=rtbis;
