This function fits a power law of the form u = aii*z^n to the given data and returns the constant coefficient aii, the exponent n, and the coef of determination cod. The user has the option to fix the value of n to 1/6 by setting 'fixn' = 1. u is the velocity (normalized by u*), z is the height above bottom (monotonically increasing, normailzed by zo), and h is the normailzed flow depth (by zo).
0001 function [aii,n,cod,upred,zpred,delta] = fitPowerLawFull(u,z,fixn,h) 0002 0003 % This function fits a power law of the form u = aii*z^n to the given data 0004 % and returns the constant coefficient aii, the exponent n, and the coef of determination cod. The user has the 0005 % option to fix the value of n to 1/6 by setting 'fixn' = 1. u is the velocity (normalized by u*), z is 0006 % the height above bottom (monotonically increasing, normailzed by zo), and h is the normailzed flow 0007 % depth (by zo). 0008 0009 % P.R. Jackson, 10-10-10 0010 0011 %Example: 0012 0013 % clear all 0014 % z = 5:50; X = 0.4512*z.^(1/6); 0015 % X = X + (2*rand(size(X))-1).*0.05.*X; 0016 % [aii,n,ssr,Xpred,zpred,delta] = fitPowerLaw(X,z,0,55,1,1,1); 0017 % figure(1); clf; plot(X,z); hold on 0018 % plot(Xpred,zpred,'r-'); hold on 0019 % plot(Xpred+delta,zpred,'r:',Xpred-delta,zpred,'r:'); hold on 0020 0021 if (nargin < 4) 0022 h = max(z); 0023 end 0024 0025 zpred = linspace(0,h,100); 0026 0027 switch fixn 0028 case 0; 0029 coefinit = [1 1/6]; %initial guess at coefficients 0030 pwrfcn = inline('coef(1)*z.^coef(2)', 'coef', 'z'); 0031 [coef,r,J,sig,mse] = nlinfit(z,u,pwrfcn,coefinit); %nonlinear fit 0032 [upred,delta] = nlpredci(pwrfcn,zpred,coef,r,'covar',sig); 0033 aii = coef(1); 0034 n = coef(2); 0035 0036 case 1; 0037 coefinit = [1]; %initial guess at coefficients 0038 pwrfcn = inline('coef(1)*z.^(1/6)', 'coef', 'z'); 0039 [coef,r,J,sig,mse] = nlinfit(z,u,pwrfcn,coefinit); %nonlinear fit 0040 [upred,delta] = nlpredci(pwrfcn,zpred,coef,r,'covar',sig); 0041 aii = coef(1); 0042 n = 1/6; 0043 end 0044 0045 ssr = sum(r.^2); 0046 0047 sstot = sum((u - mean(u)).^2); 0048 0049 cod = 1 - ssr./sstot; 0050