0001 function sta = STA_MeanProfileV2(fitProf,units,StrDir)
0002
0003
0004
0005
0006
0007
0008
0009
0010
0011
0012
0013
0014 if nargin < 3
0015 StrDir = [];
0016 end
0017
0018
0019
0020
0021
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
0050 selection = listdlg('ListSize',[300 300],'ListString', files,'Name','Select Data Files');
0051 zFileName = files(selection);
0052
0053
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
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];
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)};
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
0084
0085 j=0;
0086 figure(1); clf
0087
0088
0089 for zi=1:z
0090
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
0099 screenData=1;
0100
0101
0102
0103
0104
0105
0106
0107
0108 ignoreBS = 0;
0109 [A]=tfile(fileName,screenData,ignoreBS);
0110
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
0120 msgbox('Processing Data...','Processing Status','help','replace');
0121
0122
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
0132
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;
0137 sta.vNorth = sta.vNorth./100;
0138 sta.vmag = sta.vmag./100;
0139 sta.vVert = sta.vVert./100;
0140 end
0141
0142
0143 if 0
0144 sta.depth = nanmean(nanmean(A.Nav.depth,2));
0145 sta.vmag = nanmean(A.Wat.vMag,2)./100;
0146 indx = find(~isnan(sta.vmag));
0147 sta.vmag = sta.vmag(indx);
0148 sta.binDepth = nanmean(A.Wat.binDepth,2);
0149 sta.binDepth = sta.binDepth(indx);
0150
0151
0152
0153 sta.vNorth = nanmean(A.Wat.vNorth,2)./100;
0154 sta.vNorth = sta.vNorth(indx);
0155 sta.vEast = nanmean(A.Wat.vEast,2)./100;
0156 sta.vEast = sta.vEast(indx);
0157 end
0158
0159
0160
0161 sta.vDir = 90 - (atan2(sta.vNorth, sta.vEast))*180/pi;
0162 qindx = find(sta.vDir < 0);
0163 if ~isempty(qindx)
0164 sta.vDir(qindx) = sta.vDir(qindx) + 360;
0165 end
0166
0167
0168
0169 if isempty(StrDir)
0170 sta.mvNorth = nanmean(sta.vNorth);
0171 sta.mvEast = nanmean(sta.vEast);
0172
0173 sta.mvDir = 90 - (atan2(sta.mvNorth, sta.mvEast))*180/pi;
0174
0175 if sta.mvDir < 0
0176 sta.mvDir = sta.mvDir + 360;
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
0186 sta.psi = (sta.vDir-sta.Streamwise);
0187
0188
0189 sta.U = cosd(sta.psi).*sta.vmag;
0190 sta.V = sind(sta.psi).*sta.vmag;
0191
0192
0193 sta.z = sta.depth - sta.binDepth;
0194
0195 switch units
0196 case{'metric'}
0197 figure(1); hold on;
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;
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
0222 if fitProf
0223
0224
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
0240 sta.znorm = sta.z./sta.depth;
0241 indxfr = find(sta.znorm >= fitrng(1) & sta.znorm <= fitrng(2) & ~isnan(sta.U));
0242
0243
0244
0245 disp('Log Law Fit')
0246 disp(' u* zo cod ')
0247
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
0254 sta.ustar(zi) = ustar;
0255 sta.zo(zi) = zo;
0256 sta.cod = cod;
0257 sta.ks(zi) = 30*zo;
0258
0259 if 0
0260 sta.cod(zi) = cod;
0261 figure(1); hold on; subplot(ceil(z/3),3,zi)
0262
0263 plot(Xpred,zpred,'r-'); hold on
0264
0265
0266 grid on
0267
0268
0269 disp([ustar zo cod])
0270
0271
0272
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
0283 plot(Xpred.*ustar,zpred.*zo,'r--'); hold on
0284
0285
0286 chenprod = aii*n*2.718281828*0.41;
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
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
0302 plot(Xpred.*ustar,zpred.*zo,'r:'); hold on
0303
0304
0305 chenprod = aii*n*2.718281828*0.41;
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
0318 idx=strfind(fileName,'.');
0319 namecut = fileName(1:idx(1)-5);
0320
0321 clear latlon idx namecut
0322
0323
0324 sta.DAVstws(zi) = VMT_LayerAveMean(sta.z,sta.U);
0325
0326
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
0335 if exist('LastDir.mat') == 2
0336 save('LastDir.mat','stapath','-append')
0337 else
0338 save('LastDir.mat','stapath')
0339 end
0340
0341 end