function MYO9_13
% MYO9_13  Solution to problem 9.13 in Munson, Young, and Okiishi

% --- Given data
%       y = position  (mm)
%       h = manometer reading, (mm H20)
y = [0   2.1   4.3   6.4  10.7  15.0  19.3  23.6  26.8  29.3  32.7];
h = [0  10.6  21.1  25.6  32.5  36.9  39.4  40.5  41.0  41.0  41.0];

gammaw = 999;     %  specific weight of water,  N/m^3
rhoa = 1.23;      %  density of air

% --- Convert manometer readings to velocities
%         0.5*rhoa*u^2 = gammaw*h => u = sqrt(2*gammaw*h/rhoa)
u = sqrt(2*gammaw*(h/1000)/rhoa);   %  h/1000 is manometer reading in meters
U = max(u);                         %  Free stream velocity
plot(u,y,'o-');   xlabel('u  (m/s)');   ylabel('y (mm)');

% --- Print a nice table of data
fprintf('\n   y (mm)      h (mm)    u (m/s)       u/U\n');
for i=1:length(y)
  fprintf('%7.1f     %7.1f   %9.2f   %9.3f\n',y(i),h(i),u(i),u(i)/U);
end

% --- Find del99 by interpolation in the velocity data
iU = min(find(u==U));         %  iU is index of first element with u(y) = U
del99 = interp1(u(1:iU)/U,y(1:iU),0.99); %  Find y value where u/U = 0.99

% --- use numerical integration to compute delStar and theta
delStar = trapz(y,1-u/U);       %  Integral of (1-u/U) w.r.t. y
theta = trapz(y,(u/U).*(1-u/U));   %  Integral of (u/U)*(1-u/U) w.r.t. y
fprintf('\ndelta99 = %6.2f (mm)\ndelta*  = %6.2f (mm)\ntheta   = %6.2f (mm)\n',...
        del99,delStar,theta);

% --- Add horizontal lines at y = delta, delStar, and theta
hold on
plot([0 U],del99*[1 1],'r--');             %  horizontal line at y = del99
theText = sprintf('\\delta_{99} = %4.2f',del99);
text(0.1*U,1.05*del99,theText,'FontSize',16);    %  label the line

% --- Add horizontal line at y =
uds = interp1(y,u,delStar);              %  interpolate to find u(y)(delStar)
plot([0 U],[delStar delStar],'r--');    %  line at y = delStar
theText = sprintf('\\delta^* = %4.2f',delStar);
text(0.1*U,1.2*delStar,theText,'FontSize',16);   %  and label the line

ut = interp1(y,u,theta);                     %  interpolate to find u(theta)
plot([0 U],[theta theta],'r--');             %  line at y = theta
theText = sprintf('\\theta = %4.2f',theta);
text(0.1*U,1.35*theta,theText,'FontSize',16);        %  and label the line
hold off
xlabel('u   (m/s)');   ylabel('y   (mm)');
