% script to test log fit of velocity profiles

% 'real' values of u* (cm/s) and zo (cm):
ustr = 1.5; zo = 1;

% elevations (cm)
z = [50;40;30;20;10];
%z = [50;45;40;35;30;25;20;15;10];
vk = 0.41;

% real profile
u = ustr/vk .* log( z./zo );

% standard deviation of u measurements
sd = 0.05.*u;
% ...used to generate random errors in measured profile
randn('state',0)
rn = randn(length(z),1).*sd;
% measured profile
ur = u+rn;

figure(1)
clf
semilogy([u;0],[z;zo],'-b') 
hold on
semilogy(ur,z,'+k')
axis([0 50 .1 100])

% least-squared fit, log(z) is independent var., ur is dependent var.
[a,b,r2,sa,sb]=lsfit(log(z),ur,0)
% back out estimtes of u* and zo
ustr_hat = vk*b
lnzo_hat = -a/b
% 'measured' profile
uhat = (ustr_hat/vk)*(log(z)-lnzo_hat);
semilogy([uhat;0],[z;exp(lnzo_hat)],'--r')

% error bars on u* and ln(zo)
dustr = vk*sb
dlnzo = abs(lnzo_hat)*sqrt( sa^2/abs(a) + sb^2/abs(b) )

ebar = [exp(lnzo_hat+dlnzo),exp(lnzo_hat-dlnzo)];
h=semilogy([0 0],ebar,'-r')
set(h,'linewidth',4)