fitPowerLaw

PURPOSE ^

This function fits a power law of the form X = aii*z^n to the given data

SYNOPSIS ^

function [aii,n,ssr,Xpred,zpred,delta] = fitPowerLaw(X,z,fixn,h,mid,top,bottom)

DESCRIPTION ^

 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.

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

Generated on Thu 21-Aug-2014 10:40:31 by m2html © 2005