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: 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

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