load y.asc
load U.asc
load Umean.asc

U = Umean;

N = length(y);
N2 = floor(N/2);

nu = 1/4000;
ustar = sqrt(nu*abs((U(1)-U(2))/(y(1)-y(2))))
a = y(N);
b = y(1);
c = ustar/nu;
yplus = c*(y-a);
Uplus = U/ustar;

yplus = yplus(N2:N);
Uplus = Uplus(N2:N);
%U = (U + U(Ny:-1:1))/2;

% wall law
ywall = 1e-1: 0.5 : 19;

% core law
ycore = 2:250;

semilogx(yplus,Uplus,'k', ywall,ywall, 'b--', ycore, 2.5*log(ycore)+5,'b--');

axis([1e-1 1e3 0 25])

text(150,19.5,'2.5 ln y_+ + 5')
text(17,19.5, 'y_+');
text(16,8.5, 'U_+');
text(.3,20, strcat('u^* = ', num2str(ustar)));

dir = regexprep(pwd, '_', '-');
title(dir);

xlabel('y_+')
ylabel('U_+')

print -depsc2 walllaw.eps