fitPowerLawFull

PURPOSE ^

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

SYNOPSIS ^

function [aii,n,cod,upred,zpred,delta] = fitPowerLawFull(u,z,fixn,h)

DESCRIPTION ^

 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).

CROSS-REFERENCE INFORMATION ^

This function calls: This function is called by:

SOURCE CODE ^

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

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