# STA_MeanProfile

## PURPOSE

Compute the mean priole from stationary measurements at a point.

## SYNOPSIS

function sta = STA_MeanProfile(fitProf,units)

## DESCRIPTION

     Compute the mean priole from stationary measurements at a point.
     Currently assumes units are metric.

## CROSS-REFERENCE INFORMATION

This function calls:

- fitLogLaw This function fits a log law of the form u/u* = 1/kappa*ln(z/zo) to the given data
- fitPowerLawFull This function fits a power law of the form u = aii*z^n to the given data
- tfile tfile reads the data from an RDI ASCII output file and puts the

This function is called by:

## SOURCE CODE

    0001 function sta = STA_MeanProfile(fitProf,units)
    0002 % Compute the mean priole from stationary measurements at a point.
    0003 % Currently assumes units are metric.
    0004
    0005 % P.R. Jackson 10.12.10
    0006
    0007 %% Read and Process the Data
    0008
    0009 % Determine Files to Process
    0010 % Prompt user for directory containing files
    0011 zPathName = uigetdir('','Select the Directory Containing ASCII Output Data Files (*.txt)');
    0012 Files = dir(zPathName);
    0013 allFiles = {Files.name};
    0014 filefind=strfind(allFiles,'ASC.TXT')';
    0015 filesidx=nan(size(filefind,1),1);
    0016 for i=1:size(filefind,1)
    0017     filesidx(i,1)=size(filefind{i},1);
    0018 end
    0019 filesidx=find(filesidx>0);
    0020 files=allFiles(filesidx);
    0021
    0022 % Allow user to select which files are to be processed
    0023 selection = listdlg('ListSize',[300 300],'ListString', files,'Name','Select Data Files');
    0024 zFileName = files(selection);
    0025
    0026 % Determine number of files to be processed
    0027 if  isa(zFileName,'cell')
    0028     z=size(zFileName,2);
    0029     zFileName=sort(zFileName);
    0030 else
    0031     z=1;
    0032     zFileName={zFileName}
    0033 end
    0034 msgbox('Loading Data...Please Be Patient','Processing Status','help','replace');
    0035 % Read in Selected Files
    0036 % Initialize row counter
    0037 j=0;
    0038 figure(1); clf
    0039
    0040 % Begin master loop
    0041 for zi=1:z
    0042     % Open txt data file
    0043     if  isa(zFileName,'cell')
    0044         fileName=strcat(zPathName,'\',zFileName(zi));
    0045         fileName=char(fileName);
    0046     else
    0047         fileName=strcat(zPathName,zFileName);
    0048     end
    0049
    0050     % screenData 0 leaves bad data as -32768, 1 converts to NaN
    0051     screenData=1;
    0052
    0053     % tfile reads the data from an RDI ASCII output file and puts the
    0054     % data in a Matlab data structure with major groups of:
    0055     % Sup - supporing data
    0056     % Wat - water data
    0057     % Nav - navigation data including GPS.
    0058     % Sensor - Sensor data
    0059     % Q - discharge related data
    0060     ignoreBS = 0;
    0061     [A]=tfile(fileName,screenData,ignoreBS);
    0062     %Extract only Lat lon data
    0063     latlon(:,1)=A.Nav.lat_deg(:,:);
    0064     latlon(:,2)=A.Nav.long_deg(:,:);
    0065
    0066     switch A.Sup.units(1,:)
    0067         case{'ft'}
    0068             error('Units must be metric to start')
    0069     end
    0070
    0071     % Write the data to a file
    0072     msgbox('Processing Data...','Processing Status','help','replace');
    0073
    0074     %Rescreen data to remove missing data (30000 value)
    0075     indx1 = find(abs(latlon(:,1)) > 90);
    0076     indx2 = find(abs(latlon(:,2)) > 180);
    0077     latlon(indx1,1)=NaN;
    0078     latlon(indx2,2)=NaN;
    0079
    0080     sta.lat = nanmean(latlon(:,1));
    0081     sta.lon = nanmean(latlon(:,2));
    0082
    0083     %Compute the mean profile
    0084     sta.depth = nanmean(nanmean(A.Nav.depth,2)); %Mean Depth
    0085     sta.vmag  = nanmean(A.Wat.vMag,2)./100;
    0086     indx = find(~isnan(sta.vmag));
    0087     sta.vmag = sta.vmag(indx);
    0088     sta.binDepth = nanmean(A.Wat.binDepth,2);
    0089     sta.binDepth = sta.binDepth(indx);
    0090
    0091     % Orient with origin at the bed
    0092     sta.z = sta.depth - sta.binDepth;
    0093
    0094     switch units
    0095         case{'metric'}
    0096             figure(1); hold on; subplot(ceil(z/3),3,zi)
    0097             plot(sta.vmag,sta.z,'ko','MarkerFaceColor','k')
    0098             xlabel('Velocity (m/s)')
    0099             ylabel('Height above bottom (m)')
    0100             ylim([0 ceil(sta.depth)])
    0101         case{'english'}
    0102             sta.depth = sta.depth*3.281;
    0103             sta.vmag = sta.vmag*3.281;
    0104             sta.binDepth = sta.binDepth*3.281;
    0105             sta.z = sta.z*3.281;
    0106             figure(1); hold on; subplot(ceil(z/3),3,zi)
    0107             plot(sta.vmag,sta.z,'ko','MarkerFaceColor','k')
    0108             xlabel(upper('Velocity, in feet per second'))
    0109             ylabel(upper('Height above bottom, in feet'))
    0110             ylim([0 ceil(sta.depth)])
    0111     end
    0112
    0113
    0114     %Fit the profile
    0115     if fitProf
    0116         % First, fit a log law to get u* and zo
    0117
    0118         disp('Log Law Fit')
    0119         disp('     u*      zo      cod   ')
    0120         [ustar,zo,ks,cod,Xpred,zpred,delta] = fitLogLaw(sta.vmag',sta.z');
    0121         sta.ustar(zi) = ustar;
    0122         sta.zo(zi) = zo;
    0123         sta.ks(zi) = ks;
    0124         sta.cod(zi) = cod;
    0125         figure(1); hold on; subplot(ceil(z/3),3,zi)
    0126         %plot(X,z,'ko','MarkerFaceColor','k'); hold on
    0127          plot(Xpred,zpred,'r-'); hold on
    0128         % plot(Xpred+delta',zpred,'r:',Xpred-delta',zpred,'r:'); hold on
    0129          %plot([0 1],[sta.depth sta.depth],'k--'); hold on  %Mean flow depth for period
    0130          grid on
    0131          %ylabel('Mean Height Above Bottom (ft)')
    0132          %xlabel('Mean Velocity (ft/s)')
    0133          disp([ustar zo cod])
    0134
    0135
    0136         % Now, fit with power law (fixed at 1/6)
    0137
    0138         disp('Fixed Exponent Power Law Fit')
    0139         disp('     a      n      cod     ChenProd')
    0140
    0141         [aii,n,cod,Xpred,zpred,delta] = fitPowerLawFull(sta.vmag'./ustar,sta.z'./zo,1,sta.depth./zo);
    0142
    0143
    0144
    0145          figure(1); hold on; subplot(ceil(z/3),3,zi)
    0146          % plot(X,z,'bo','MarkerFaceColor','b'); hold on
    0147          plot(Xpred.*ustar,zpred.*zo,'r--'); hold on
    0148          %plot(Xpred+delta,zpred,'r:',Xpred-delta,zpred,'r:'); hold on
    0149          %plot([0 1],[25.00 25.00],'k--'); hold on
    0150          chenprod = aii*n*2.718281828*0.41; %product of a*n*e*kappa(Von Karman cst) == 1 for exact fit with log law
    0151           disp([aii n cod chenprod])
    0152           sta.aii_p1(zi) = aii;
    0153           sta.n_p1(zi) = n;
    0154           sta.cod_p1(zi) = cod;
    0155           sta.cp_p1(zi) = chenprod;
    0156
    0157          % Now, fit with power law (variable exponent)
    0158
    0159         disp('Variable Exponent Power Law Fit')
    0160         disp('     a      n      cod     ChenProd')
    0161
    0162         [aii,n,cod,Xpred,zpred,delta] = fitPowerLawFull(sta.vmag'./ustar,sta.z'./zo,0,sta.depth./zo);
    0163
    0164          figure(1); hold on; subplot(ceil(z/3),3,zi)
    0165          % plot(X,z,'bo','MarkerFaceColor','b'); hold on
    0166          plot(Xpred.*ustar,zpred.*zo,'r:'); hold on
    0167          %plot(Xpred+delta,zpred,'r:',Xpred-delta,zpred,'r:'); hold on
    0168          %plot([0 1],[25.00 25.00],'k--'); hold on
    0169          chenprod = aii*n*2.718281828*0.41; %product of a*n*e*kappa(Von Karman cst) == 1 for exact fit with log law
    0170          disp([aii n cod chenprod])
    0171           sta.aii_p2(zi) = aii;
    0172           sta.n_p2(zi) = n;
    0173           sta.cod_p2(zi) = cod;
    0174           sta.cp_p2(zi) = chenprod;
    0175
    0176     end
    0177
    0178     clear A
    0179
    0180     % Determine the file name
    0181     idx=strfind(fileName,'.');
    0182     namecut = fileName(1:idx(1)-5);
    0183
    0184     clear latlon idx namecut
    0185
    0186     % Compute the depth average velocity
    0187     sta.dam(zi) = -1*VMT_LayerAveMean(sta.z,sta.vmag);
    0188 end
    0189
    0190 msgbox('Processing Complete','Processing Status','help','replace');
    0191
    0192
    0193 end