function MYO_11_38
% MYO_11_38  Solution to problem 11.38 in Munson, Young and Okiishi

% --- Constants for this analysis
astar = 0.1;                     %  minimum area,  ft^2
k = 1.66;                        %  specific heat ratio for Helium

% --- Geometry of the duct
%       Notice that x = 0 (exactly) is avoided.  This
%       prevents blow up of formulas when Ma = 1 (exactly)
xsub = linspace(-1,1e-4,25)';   %  Position along the duct on subsonic side
xsup = linspace(1e-4,1,25)';    %  Position along the duct on subsonic side
x = [xsub; xsup];               %  Complete x vector;  xsub and xsup are no longer needed
r = (astar + xsub.^2)/pi;       %  Duct radius
aas = 1 + x.^2/astar;           %  Area ratio, A/Astar, as function of position

% --- Compute Mach number as a function of position along the duct
%     The area is given, so Ma is found as a function of A/Astar
%     This requires selecting the appropriate branch of the Ma = f(A/Astar)
%     curve depending on whether the flow is subsonic or supersonic
%     The isenMaar() function finds Ma = f(A/Astar) for isentropic flow.
%
%     Two solutions are found for the flow through the converging-diverging
%     nozzle.  Ma1 is the Mach number as a function of x for subsonic flow
%     in the converging section and supersonic flow in the diverging section.
%     Ma2 is the solution for subsonic flow througout.  The subsonic or
%     supersonic solution is selected by the second input argument to
%     isenMaar.  This argument is a flag:  negative values select the
%     subsonic branch and positive values select the supersonic branch.
%     This is especially convenient for the first solution (subsonic then
%     supersonic) since the x coordinate is zero at the throat, negative
%     upstream of the throat (where the flow is subsonic) and positive
%     downstream of the throat.  Thus the value of x(j) is used as the
%     flag even though the value of x is not a true parameter in
%     Ma = f(A/Astar).
Ma1 = zeros(size(x));    %  first solution:  subsonic, then supersonic
Ma2 = zeros(size(x));    %  second solution:  subsonic throughout
for j=1:length(x)
  Ma1(j) = isenMaaas(aas(j),x(j),k); %  x(j)<0 for subsonic, x(j)>0 for supersonic
  Ma2(j) = isenMaaas(aas(j),-1,k);   %  subsonic flow throughout the nozzle
end

% --- First set of curves:  flow is subsonic for x<0 (before throat)
%     and supersonic after throat.  The Ma number 
tr = isenTT0(Ma1,k);    %  compute T/T0 for given Ma
pr = isenpp0(Ma1,k);    %  p/p0
rr = isenrr0(Ma1,k);    %  rho/rho0

% --- Second set of curves: subsonic flow throughout, Ma2 <= 1 for all x
tr2 = isenTT0(Ma2,k);    %  compute T/T0 for given Ma
pr2 = isenpp0(Ma2,k);    %  p/p0
rr2 = isenrr0(Ma2,k);    %  rho/rho0

% --- plot Ma versus x for two solutions
plot(x,Ma1,'k-',x,Ma2,'k-');   xlabel('x   (ft)');   ylabel('Ma');
title('Munson, Young, and Okiishi, problem 11.38');

% --- plot T/T0, p/p0 and rho/rho0 versus x
f = figure;
plot(x,tr,'k-',x,pr,'r-.',x,rr,'b:',x,tr2,'k-',x,pr2,'r-.',x,rr2,'b:')
legend('T/T_0','p/p_0','\rho/\rho_0',3)
xlabel('x   (ft)');
title('Munson, Young, and Okiishi, problem 11.38');
