This function fits a power law of the form X = aii*z^n to the given data and returns the constant coefficient aii, the exponent n, and the sum of the squared residuals ssr. The user has the option to fix the value of n to 1/6 by setting 'fixn' = 1. X is the velocity cross product, z is the height above bottom (monotonically increasing), and h is the flow depth; mid, top and bottom are TF (1,0) values that tell the code to fit the middle of the profile, top of the profile, and bottom of the profile, respectively.
0001 function [aii,n,ssr,Xpred,zpred,delta] = fitPowerLaw(X,z,fixn,h,mid,top,bottom) 0002 0003 % This function fits a power law of the form X = aii*z^n to the given data 0004 % and returns the constant coefficient aii, the exponent n, and the sum of the squared residuals ssr. The user has the 0005 % option to fix the value of n to 1/6 by setting 'fixn' = 1. X is the velocity cross product, z is 0006 % the height above bottom (monotonically increasing), and h is the flow 0007 % depth; mid, top and bottom are TF (1,0) values that tell the code to fit the 0008 % middle of the profile, top of the profile, and bottom of the profile, respectively. 0009 0010 % P.R. Jackson, 5-1-08 0011 0012 %Example: 0013 0014 % clear all 0015 % z = 5:50; X = 0.4512*z.^(1/6); 0016 % X = X + (2*rand(size(X))-1).*0.05.*X; 0017 % [aii,n,ssr,Xpred,zpred,delta] = fitPowerLaw(X,z,0,55,1,1,1); 0018 % figure(1); clf; plot(X,z); hold on 0019 % plot(Xpred,zpred,'r-'); hold on 0020 % plot(Xpred+delta,zpred,'r:',Xpred-delta,zpred,'r:'); hold on 0021 0022 zpred = []; 0023 0024 switch mid 0025 case 0; zpred = []; 0026 case 1; zpred = z; 0027 end 0028 0029 switch top 0030 case 0; zpred = zpred; 0031 case 1; zpred = horzcat(zpred,linspace(z(end),h,10))'; 0032 end 0033 0034 switch bottom 0035 case 0; zpred = zpred; 0036 case 1; if top; zpred = horzcat(zpred',linspace(0,z(1),100))'; 0037 else zpred = horzcat(zpred,linspace(0,z(1),100))'; 0038 end 0039 end 0040 0041 zpred = sort(zpred); 0042 0043 switch fixn 0044 case 0; 0045 coefinit = [1 1/6]; %initial guess at coefficients 0046 pwrfcn = inline('coef(1)*z.^coef(2)', 'coef', 'z'); 0047 [coef,r,J,sig,mse] = nlinfit(z,X,pwrfcn,coefinit); %nonlinear fit 0048 [Xpred,delta] = nlpredci(pwrfcn,zpred,coef,r,'covar',sig); 0049 aii = coef(1); 0050 n = coef(2); 0051 0052 case 1; 0053 coefinit = [1]; %initial guess at coefficients 0054 pwrfcn = inline('coef(1)*z.^(1/6)', 'coef', 'z'); 0055 [coef,r,J,sig,mse] = nlinfit(z,X,pwrfcn,coefinit); %nonlinear fit 0056 [Xpred,delta] = nlpredci(pwrfcn,zpred,coef,r,'covar',sig); 0057 aii = coef(1); 0058 n = 1/6; 0059 end 0060 0061 ssr = sum(r.^2); 0062 0063 0064