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

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