# 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