function[Mach_Number]=Newtons_method_Channel(A,Ao,Mo)

%function[Mach_Number]=Newtons_method_Channel(A,Ao,Mo)
%
%This function will attempt to return the mach number at a location along
%the flow channel composed of a small circular arc using Isentropic
%relations and Newton's method 
%
%A   is the Area at the given x value
%Ao  is the Area at the entrance
%Mo  is the Mach Number at the entrance
    
gamma=1.4;                        %Ratio of specific heats
gamma_fac=(gamma-1)/2;            %Coefficient appearing in formula
gamma_exp=(gamma+1)/(gamma-1);    %Exponent appearing in formula
%Initial guess at Mach Number
%M=Mach Number

if Mo < 1 
   M=.001*ones(1,length(A));
else 
M=.001*ones(1,length(A));

%I just figured out what the problem is----Some of these values have two
%roots.  So Newton's method will hop around
num_steps=20

for i=1:num_steps
  f=(A./Ao).^2-(Mo./M).^2.*(((1+gamma_fac*M.^2)./(1+gamma_fac*Mo.^2)).^gamma_exp);
  f_p=-Mo^2/((1+gamma_fac*Mo.^2).^gamma_exp).*...
         (-(2./M.^3).*(1+gamma_fac*M.^2).^gamma_exp+((gamma+1)./M).*(1+gamma_fac*M.^2).^(2./(gamma-1)));
  f_p;
  M=M-f./f_p;   
end

Mach_Number=M;