# STA_MeanProfileV2

## PURPOSE

Compute the mean profile from stationary measurements at a point.

## SYNOPSIS

function sta = STA_MeanProfileV2(fitProf,units,StrDir)

## DESCRIPTION

     Compute the mean profile from stationary measurements at a point.
     Currently assumes incomming units are metric--can convert using the 'units' input.

## CROSS-REFERENCE INFORMATION

This function calls:

- ComputeNormalizedProfile This function computes the normalized velocity profile from basic velocity
- fitLogLawV2 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_MeanProfileV2(fitProf,units,StrDir)
    0002 % Compute the mean profile from stationary measurements at a point.
    0003 % Currently assumes incomming units are metric--can convert using the 'units' input.
    0004
    0005 % fitProf = 1 for log and power law fitting of the profile
    0006 % units = 'metric or 'english' (for output plots)
    0007
    0008 %V2 modifies the original code to use the streamwise component of velocity
    0009 %rather than the velocity magnitude. It aslo allows the user to fit only a range of the flow profile.
    0010 %P.R. Jackson, 2-10-11
    0011
    0012 % P.R. Jackson 10.12.10
    0013
    0014 if nargin < 3
    0015     StrDir = [];
    0016 end
    0017
    0018 %% Read and Process the Data
    0019
    0020 % Determine Files to Process
    0021 % Prompt user for directory containing files
    0022 defaultpath = 'C:\';
    0023 stapath = [];
    0024 if exist('VMT\LastDir.mat') == 2
    0025     load('VMT\LastDir.mat');
    0026     if exist(stapath) == 7
    0027         stapath = uigetdir(stapath,'Select the Directory Containing ASCII Output Data Files (*.txt)');
    0028     else
    0029         stapath = uigetdir(defaultpath,'Select the Directory Containing ASCII Output Data Files (*.txt)');
    0030     end
    0031 else
    0032     stapath = uigetdir(defaultpath,'Select the Directory Containing ASCII Output Data Files (*.txt)');
    0033 end
    0034 zPathName = stapath;
    0035 Files = dir(zPathName);
    0036 allFiles = {Files.name};
    0037 filefind = strfind(allFiles,'ASC.TXT')';
    0038 filesidx=nan(size(filefind,1),1);
    0039 for i=1:size(filefind,1)
    0040     filesidx(i,1)=size(filefind{i},1);
    0041 end
    0042 filesidx=find(filesidx>0);
    0043 files=allFiles(filesidx);
    0044
    0045 if isempty(files)
    0046     errordlg(['No ASC.TXT files found in ' stapath '.  Ensure data files are in the form "*_ASC.TXT" (Standard WRII naming convention).']);
    0047 end
    0048
    0049 % Allow user to select which files are to be processed
    0050 selection = listdlg('ListSize',[300 300],'ListString', files,'Name','Select Data Files');
    0051 zFileName = files(selection);
    0052
    0053 % Determine number of files to be processed
    0054 if  isa(zFileName,'cell')
    0055     z=size(zFileName,2);
    0056     zFileName=sort(zFileName);
    0057 else
    0058     z=1;
    0059     zFileName={zFileName}
    0060 end
    0061 msgbox('Loading Data...Please Be Patient','Processing Status','help','replace');
    0062
    0063 % Set the fitting region (new 2-10-11, PRJ)
    0064 if fitProf
    0065     choice = questdlg('Do you want to fit the full profile?','Fitting Setup',...
    0066         'Yes','No','Yes');
    0067         switch choice
    0068             case 'Yes'
    0069                 fitrng = [0 1]; %z/H (normalized) so z/H = 0 at bed and 1 at surface
    0070             case 'No'
    0071                 disp('Enter Normalized (z/H) Fit Range')
    0072                 prompt = {'Lower Bound','Upper Bound'};
    0073                 dlg_title = 'Fitting Setup';
    0074                 num_lines = 1;
    0075                 def = {num2str(0),num2str(0.2)};  %Wall region only by default (Nezu and Nakagawa 1993)
    0076                 answer = inputdlg(prompt,dlg_title,num_lines,def);
    0077                 [ftl, status1] = str2num(answer{1});
    0078                 [ftu, status2] = str2num(answer{2});
    0079                 fitrng = [ftl ftu];
    0080         end
    0081 end
    0082
    0083 % Read in Selected Files
    0084 % Initialize row counter
    0085 j=0;
    0086 figure(1); clf
    0087
    0088 % Begin master loop
    0089 for zi=1:z
    0090     % Open txt data file
    0091     if  isa(zFileName,'cell')
    0092         fileName=strcat(zPathName,'\',zFileName(zi));
    0093         fileName=char(fileName);
    0094     else
    0095         fileName=strcat(zPathName,zFileName);
    0096     end
    0097
    0098     % screenData 0 leaves bad data as -32768, 1 converts to NaN
    0099     screenData=1;
    0100
    0101     % tfile reads the data from an RDI ASCII output file and puts the
    0102     % data in a Matlab data structure with major groups of:
    0103     % Sup - supporing data
    0104     % Wat - water data
    0105     % Nav - navigation data including GPS.
    0106     % Sensor - Sensor data
    0107     % Q - discharge related data
    0108     ignoreBS = 0;
    0109     [A]=tfile(fileName,screenData,ignoreBS);
    0110     %Extract only Lat lon data
    0111     latlon(:,1)=A.Nav.lat_deg(:,:);
    0112     latlon(:,2)=A.Nav.long_deg(:,:);
    0113
    0114     switch A.Sup.units(1,:)
    0115         case{'ft'}
    0116             error('Units must be metric to start')
    0117     end
    0118
    0119     % Write the data to a file
    0120     msgbox('Processing Data...','Processing Status','help','replace');
    0121
    0122     %Rescreen data to remove missing data (30000 value)
    0123     indx1 = find(abs(latlon(:,1)) > 90);
    0124     indx2 = find(abs(latlon(:,2)) > 180);
    0125     latlon(indx1,1)=NaN;
    0126     latlon(indx2,2)=NaN;
    0127
    0128     sta.lat = nanmean(latlon(:,1));
    0129     sta.lon = nanmean(latlon(:,2));
    0130
    0131     %Compute the median profile using a normalized profile (added 3-20-12,
    0132     %PRJ)
    0133     if 1
    0134         [sta.binDepth,zm,sta.vEast,sta.vNorth,sta.vVert,sta.vmag,sta.obsav,sta.depth] = ComputeNormalizedProfile(nanmean(A.Nav.depth,2)',A.Wat.binDepth,A.Wat.vEast,A.Wat.vNorth,A.Wat.vVert);
    0135         sta.binDepth = sta.binDepth';
    0136         sta.vEast = sta.vEast./100; %convert to m/s from cm/s
    0137         sta.vNorth = sta.vNorth./100; %convert to m/s from cm/s
    0138         sta.vmag = sta.vmag./100; %convert to m/s from cm/s
    0139         sta.vVert = sta.vVert./100; %convert to m/s from cm/s
    0140     end
    0141
    0142     %OLD METHOD--Compute the mean profile (magnitude)
    0143     if 0
    0144         sta.depth = nanmean(nanmean(A.Nav.depth,2)); %Mean Depth
    0145         sta.vmag  = nanmean(A.Wat.vMag,2)./100;  %Mean v magnitude in m/s
    0146         indx = find(~isnan(sta.vmag));
    0147         sta.vmag = sta.vmag(indx); % Use only data with no nans
    0148         sta.binDepth = nanmean(A.Wat.binDepth,2);
    0149         sta.binDepth = sta.binDepth(indx);
    0150
    0151         %Compute the streamwise component (need Vnorth, Veast to get average
    0152         %Vdir
    0153         sta.vNorth  = nanmean(A.Wat.vNorth,2)./100;  %Mean v north in m/s
    0154         sta.vNorth  = sta.vNorth(indx); % Use only data with no nans
    0155         sta.vEast   = nanmean(A.Wat.vEast,2)./100;  %Mean v east in m/s
    0156         sta.vEast   = sta.vEast(indx); % Use only data with no nans
    0157     end
    0158
    0159     % Compute the average direction from mean ve and vn (so i don't have to
    0160     % average directions on a 0-360 scale)
    0161     sta.vDir = 90 - (atan2(sta.vNorth, sta.vEast))*180/pi; %Compute the atan from the velocity componentes, convert to radians, and rotate to north axis
    0162     qindx = find(sta.vDir < 0);
    0163     if ~isempty(qindx)
    0164         sta.vDir(qindx) = sta.vDir(qindx) + 360;  %Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis
    0165     end
    0166
    0167     % Define the streamwise direction as the mean flow direction unless
    0168     % provided by user
    0169     if isempty(StrDir)
    0170         sta.mvNorth = nanmean(sta.vNorth);  %Mean vnorth (single value)
    0171         sta.mvEast  = nanmean(sta.vEast);   %Mean veast (single value)
    0172         %Compute the overall mean flow direction
    0173         sta.mvDir = 90 - (atan2(sta.mvNorth, sta.mvEast))*180/pi; %Compute the atan from the velocity components, convert to radians, and rotate to north axis
    0174
    0175         if sta.mvDir < 0
    0176             sta.mvDir = sta.mvDir + 360;  %Must add 360 deg to Quadrant 4 values as they are negative angles from the +y axis
    0177         end
    0178         sta.Streamwise = sta.mvDir;
    0179         disp(['Using Streamwise Direction of ' num2str(sta.Streamwise) ' degrees from true north'])
    0180     else
    0181         sta.Streamwise = StrDir;
    0182         disp(['Using Streamwise Direction of ' num2str(sta.Streamwise) ' degrees from true north'])
    0183     end
    0184
    0185     % Determine the deviation of a vector from streamwise velocity
    0186     sta.psi = (sta.vDir-sta.Streamwise);
    0187
    0188     % Determine streamwise (U) and transverse (V)
    0189     sta.U = cosd(sta.psi).*sta.vmag;
    0190     sta.V = sind(sta.psi).*sta.vmag;
    0191
    0192     % Orient with origin at the bed
    0193     sta.z = sta.depth - sta.binDepth;
    0194
    0195     switch units
    0196         case{'metric'}
    0197             figure(1); hold on; %subplot(ceil(z/3),3,zi)
    0198             plot(sta.U,sta.z,'ko','MarkerFaceColor','k'); hold on
    0199             xlabel('Velocity (m/s)')
    0200             ylabel('Height above bottom (m)')
    0201             ylim([0 ceil(sta.depth)])
    0202             xlim([0 max(sta.U)+0.1])
    0203             plot([0 max(sta.U)+0.1],[sta.depth sta.depth],'k--')
    0204             text((max(sta.U)+0.1)/2,sta.depth,'Water Surface')
    0205         case{'english'}
    0206             sta.depth = sta.depth*3.281;
    0207             sta.U = sta.U*3.281;
    0208             sta.binDepth = sta.binDepth*3.281;
    0209             sta.z = sta.z*3.281;
    0210             figure(1); hold on; %subplot(ceil(z/3),3,zi)
    0211             plot(sta.U,sta.z,'ko','MarkerFaceColor','k'); hold on
    0212             xlabel(upper('Velocity, in feet per second'))
    0213             ylabel(upper('Height above bottom, in feet'))
    0214             ylim([0 ceil(sta.depth)])
    0215             xlim([0 max(sta.U)+0.1])
    0216             plot([0 max(sta.U)+0.1],[sta.depth sta.depth],'k--')
    0217             text((max(sta.U)+0.1)/2,sta.depth,'Water Surface')
    0218     end
    0219
    0220
    0221     %Fit the profile
    0222     if fitProf
    0223
    0224         % Determine if any manual editing is necessary
    0225         choice = questdlg('Do you want to manually edit the data?','Data Editor',...
    0226         'Yes','No','No');
    0227         switch choice
    0228             case 'Yes'
    0229                 manedit = 1;
    0230             case 'No'
    0231                 manedit = 0;
    0232         end
    0233
    0234         if manedit
    0235             disp('Select profile data points to remove using mouse')
    0236             [sta.U,sel_indx] = DataEditor(sta.U);
    0237         end
    0238
    0239         % Apply the fitting range
    0240         sta.znorm = sta.z./sta.depth;
    0241         indxfr = find(sta.znorm >= fitrng(1) & sta.znorm <= fitrng(2) & ~isnan(sta.U));
    0242
    0243         % First, fit a log law to get u* and zo
    0244
    0245         disp('Log Law Fit')
    0246         disp('   u*       zo        cod   ')
    0247         %figure(1); clf; plot(sta.U(indxfr),sta.z(indxfr),'ko-'); pause
    0248         [ustar,zo,cod] = fitLogLawV2(sta.U(indxfr),sta.z(indxfr),sta.depth);
    0249         disp([num2str(ustar,3) '   ' num2str(zo,3) '      ' num2str(cod,3)])
    0250         zeval = linspace(0,sta.depth);
    0251         ueval = ustar./0.41.*log(zeval./zo);
    0252         figure(1); hold on; plot(ueval,zeval,'r-')
    0253         %[ustar,zo,ks,cod,Xpred,zpred,delta] = fitLogLaw(sta.U(indxfr)',sta.z(indxfr)');
    0254         sta.ustar(zi) = ustar;
    0255         sta.zo(zi) = zo;
    0256         sta.cod = cod;
    0257         sta.ks(zi) = 30*zo;
    0258         %sta.ks(zi) = ks;
    0259         if 0
    0260             sta.cod(zi) = cod;
    0261             figure(1); hold on; subplot(ceil(z/3),3,zi)
    0262             %plot(X,z,'ko','MarkerFaceColor','k'); hold on
    0263              plot(Xpred,zpred,'r-'); hold on
    0264             % plot(Xpred+delta',zpred,'r:',Xpred-delta',zpred,'r:'); hold on
    0265              %plot([0 1],[sta.depth sta.depth],'k--'); hold on  %Mean flow depth for period
    0266              grid on
    0267              %ylabel('Mean Height Above Bottom (ft)')
    0268              %xlabel('Mean Velocity (ft/s)')
    0269              disp([ustar zo cod])
    0270
    0271
    0272             % Now, fit with power law (fixed at 1/6)
    0273
    0274             disp('Fixed Exponent Power Law Fit')
    0275             disp('     a      n      cod     ChenProd')
    0276
    0277             [aii,n,cod,Xpred,zpred,delta] = fitPowerLawFull(sta.U'./ustar,sta.z'./zo,1,sta.depth./zo);
    0278
    0279
    0280
    0281              figure(1); hold on; subplot(ceil(z/3),3,zi)
    0282              % plot(X,z,'bo','MarkerFaceColor','b'); hold on
    0283              plot(Xpred.*ustar,zpred.*zo,'r--'); hold on
    0284              %plot(Xpred+delta,zpred,'r:',Xpred-delta,zpred,'r:'); hold on
    0285              %plot([0 1],[25.00 25.00],'k--'); hold on
    0286              chenprod = aii*n*2.718281828*0.41; %product of a*n*e*kappa(Von Karman cst) == 1 for exact fit with log law
    0287               disp([aii n cod chenprod])
    0288               sta.aii_p1(zi) = aii;
    0289               sta.n_p1(zi) = n;
    0290               sta.cod_p1(zi) = cod;
    0291               sta.cp_p1(zi) = chenprod;
    0292
    0293              % Now, fit with power law (variable exponent)
    0294
    0295             disp('Variable Exponent Power Law Fit')
    0296             disp('     a      n      cod     ChenProd')
    0297
    0298             [aii,n,cod,Xpred,zpred,delta] = fitPowerLawFull(sta.U'./ustar,sta.z'./zo,0,sta.depth./zo);
    0299
    0300              figure(1); hold on; subplot(ceil(z/3),3,zi)
    0301              % plot(X,z,'bo','MarkerFaceColor','b'); hold on
    0302              plot(Xpred.*ustar,zpred.*zo,'r:'); hold on
    0303              %plot(Xpred+delta,zpred,'r:',Xpred-delta,zpred,'r:'); hold on
    0304              %plot([0 1],[25.00 25.00],'k--'); hold on
    0305              chenprod = aii*n*2.718281828*0.41; %product of a*n*e*kappa(Von Karman cst) == 1 for exact fit with log law
    0306              disp([aii n cod chenprod])
    0307               sta.aii_p2(zi) = aii;
    0308               sta.n_p2(zi) = n;
    0309               sta.cod_p2(zi) = cod;
    0310               sta.cp_p2(zi) = chenprod;
    0311         end
    0312
    0313     end
    0314
    0315     clear A
    0316
    0317     % Determine the file name
    0318     idx=strfind(fileName,'.');
    0319     namecut = fileName(1:idx(1)-5);
    0320
    0321     clear latlon idx namecut
    0322
    0323     % Compute the depth average velocity (streamwise)
    0324     sta.DAVstws(zi) = VMT_LayerAveMean(sta.z,sta.U);
    0325
    0326     % Compute the depth average vertical velocity
    0327     sta.DAVvert(zi) = VMT_LayerAveMean(sta.z,sta.vVert);
    0328
    0329 end
    0330
    0331 msgbox('Processing Complete','Processing Status','help','replace');
    0332
    0333
    0334 % Save the paths
    0335 if exist('LastDir.mat') == 2
    0336     save('LastDir.mat','stapath','-append')
    0337 else
    0338     save('LastDir.mat','stapath')
    0339 end
    0340
    0341 end