% This script calculates average velocities to use in index curves. % It is assumed that a run structure % has been loaded by MatLab. It is assumed that the column averages % have already been carried out and appended to the run structure. % RWB, June 26, 1998. L1adcp=run.par.L1adcp; % Distance from west transect buoy to shore (m). L2adcp=run.par.L2adcp; % Distance from east transect buoy to shore (m). L1uvm=run.par.L1uvm; % Distance from west pier to shore (m). L2uvm=run.par.L2uvm; % Distance from east pier to shore (m). d0uvm=run.par.d0uvm; % Depth of UVM path below sea level. thchan=run.par.thchan; % Bearing of channel centerline, in degrees. thdev=run.par.devmag; % Magnetic deviation, in degrees. % Interpolate for logger values of stage and UVM velocity. run.hcs=interp1(run.cs20.pstseccs,run.cs20.hcs,run.pstsecav); run.vcs=interp1(run.cs20.pstseccs,run.cs20.vcs,run.pstsecav); % % Load in scan codes for transects. % infile=run.infile; % Name of last file in run. % scanfile=[infile(1:length(infile)-8),'scan']; % load (scanfile) itran1=1; itran2=run.ntran; ntran=itran2-itran1+1; run.vavcen=zeros(ntran,1); % Average velocities over the transect. run.nbin=zeros(ntran,1); % Number of bins used to calculate vavcen. run.vavfull=zeros(ntran,1); % Average velocity including missing ends. run.vaveast=zeros(ntran,1); % Average velocity for missing east end. run.vavwest=zeros(ntran,1); % Average velocity for missing est end. run.vavedge=zeros(ntran,1); % vaveast + vavwest run.vavuvm=zeros(ntran,1); % Average velocity along UVM path. run.vavtran=zeros(ntran,1); % Transect number for this good transect. run.vavtsec=zeros(ntran,1); % Time for this good transect. itrgood=0; for itran=itran1:itran2 if run.goodtran(itran)>=1 % Test scan code for good transects. itrgood=itrgood+1; rlandau=sqrt(run.madegood(itran)*run.par.lambda); % Landau radius. % In the following, the second component is a bin counter. vcen=zeros(2,1); % velocity in center segment. veast=zeros(2,1); % velocity on east edge. vwest=zeros(2,1); % velocity on west edge. vuvm=zeros(2,1); % velocity on UVM path. % The following eight are accumulators for linear fits % to the velocity dists. near the two edges. s1v=0; s1d=0; s1d2=0; s1vd=0; n1=0; s2v=0; s2d=0; s2d2=0; s2vd=0; n2=0; % Determine west end (#1) and east end (#2) of path. if run.dist(1,itran)==0 i1ens=1; i2ens=run.nens(itran); else i1ens=run.nens(itran); i2ens=1; end % Calculate start and end distances for the UVM path. throtrad=-(thchan-thdev)*pi/180.; % Angle of rotation from y=mag north to % y = along channel axis. mgperp=(run.xens(1,i2ens,itran)-run.xens(1,i1ens,itran))* ... cos(throtrad) + (run.xens(2,i2ens,itran)-run.xens(2,i1ens,itran))* ... sin(throtrad); Lstart=(L1uvm-L1adcp)/mgperp*run.madegood(itran); Lend=(1.-(L2uvm-L2adcp)/mgperp)*run.madegood(itran); % Depth of UVM path below free water surface. duvm=d0uvm+run.hcs(itran); for iens=1:run.nens(itran) if run.vflag(iens,itran)>=0 % exclude bad data and empty % ensembles. f1=1.; f2=1.; % Correction factors for missing top and bottom layers. % A y^1.6 power law for v is assumed. y1=run.dbot(iens,itran)-run.depth(run.WNmax(iens,itran),itran)- ... run.ens.WS/2.; % height of lower edge of measured layer. y2=run.dbot(iens,itran)-run.depth(1,itran)+ ... run.ens.WS/2.; % height of upper edge of measured layer. f1=run.dbot(iens,itran)^1.16667/(y2^1.16667-y1^1.16667); f2=run.dbot(iens,itran)/(y2-y1); if (y1<0)|(y2<0) [y1,y2] pause end if run.dist(iens,itran)(run.madegood(itran)-L2adcp) s2v=s2v+run.vpav(iens,itran)*f1; s2d=s2d+run.dist(iens,itran); s2d2=s2d2+run.dist(iens,itran)^2; s2vd=s2vd+run.vpav(iens,itran)*f1*run.dist(iens,itran); n2=n2+1; % [run.vpav(iens,itran),run.WNmax(iens,itran), ... % run.dist(iens,itran)] end vcen(1)=vcen(1)+run.vpav(iens,itran)*f1*run.WNmax(iens,itran)*f2; vcen(2)=vcen(2)+run.WNmax(iens,itran)*f2; %Sums for UVM line-average velocity. if (run.dist(iens,itran)>Lstart)&(run.dist(iens,itran)=0) rsound=2.*rlandau*min(run.dist(iens,itran)-Lstart, ... Lend-run.dist(iens,itran))/(Lend-Lstart); rsound = max(rsound,ens.WS); % At least +/- one bin. ibin1=fix((duvm-rsound-run.depth(1,itran))/run.ens.WS+1.5); ibin2=fix((duvm+rsound-run.depth(1,itran))/run.ens.WS+1.5); ibin1=max(1,min(run.WNmax(iens,itran),ibin1)); ibin2=max(1,min(run.WNmax(iens,itran),ibin2)); dl=run.dist(iens,itran)-run.dist(max(1,iens-1),itran); vuvm(1)=vuvm(1)+mean(run.vp(ibin1:ibin2,iens,itran))*dl; vuvm(2)=vuvm(2)+dl; end end end %%%%%% End of transect; evaluate various averages. % Linear fits for missing end sections. s1d=s1d-L1adcp; s1d2=s1d2+L1adcp^2; n1=n1+1; a1=(s1vd*s1d-s1v*s1d2)/(s1d^2-n1*s1d2); b1=(s1v*s1d-n1*s1vd)/(s1d^2-n1*s1d2); v1av=a1-b1/3*L1adcp; nbin1=run.WNmax(i1ens,itran)/2.* ... run.nens(itran)/run.madegood(itran)*L1adcp; vwest(1)=v1av*nbin1; vwest(2)=nbin1; % 's1d, s1d2, s1v, s1vd, a1, b1, v1av, nbin1 = ' % [s1d,s1d2,s1v,s1vd,a1,b1,v1av,nbin1] s2d=s2d+L2adcp+run.dist(i2ens,itran); s2d2=s2d2+(L2adcp+run.dist(i2ens,itran))^2; n2=n2+1; a2=(s2vd*s2d-s2v*s2d2)/(s2d^2-n2*s2d2); b2=(s2v*s2d-n2*s2vd)/(s2d^2-n2*s2d2); dlast=run.dist(i2ens,itran); v2av=(-b2/3*((dlast+L2adcp)^3-dlast^3)+ ... (b2*(dlast+L2adcp)-a2)/2*((dlast+L2adcp)^2-dlast^2) + ... a2*(dlast+L2adcp)*L2adcp)/ ... ((dlast+L2adcp)*L2adcp-((dlast+L2adcp)^2-dlast^2)/2); nbin2=run.WNmax(i2ens,itran)/2.* ... run.nens(itran)/run.madegood(itran)*L2adcp; veast(1)=v2av*nbin2; veast(2)=nbin2; % 's2d, s2d2, s2v, s2vd, a2, b2, v2av, nbin2 = ' % [s2d,s2d2,s2v,s2vd,a2,b2,v2av,nbin2] % pause vedge=veast+vwest; vfull=vcen+vedge; if vcen(1)==0, pause,end %%%%%% Carry out averages. if vfull(2)~=0, run.vavfull(itrgood)=vfull(1)/vfull(2);, ... else run.vavfull(itrgood)=run.vavfull(itrgood-1);, ... ['transect ',int2str(itran),': vfull(2)=0'],end if veast(2)~=0, run.vaveast(itrgood)=veast(1)/veast(2);, ... else run.vaveast(itrgood)=run.vaveast(itrgood-1);, ... ['transect ',int2str(itran),': veast(2)=0'],end if vwest(2)~=0, run.vavwest(itrgood)=vwest(1)/vwest(2);, ... else run.vavwest(itrgood)=run.vavwest(itrgood-1);, ... ['transect ',int2str(itran),': vwest(2)=0'],end if vedge(2)~=0, run.vavedge(itrgood)=vedge(1)/vedge(2);, ... else run.vavedge(itrgood)=run.vavedge(itrgood-1);, ... ['transect ',int2str(itran),': vedge(2)=0'],end if vcen(2)~=0, run.vavcen(itrgood)=vcen(1)/vcen(2);, ... else run.vavcen(itrgood)=run.vavcen(itrgood-1);, ... ['transect ',int2str(itran),': vcen(2)=0'],end run.nbin(itrgood)=vcen(2); if vuvm(2)~=0, run.vavuvm(itrgood)=vuvm(1)/vuvm(2);, ... else run.vavuvm(itrgood)=run.vavuvm(itrgood-1);, ... ['transect ',int2str(itran),': vuvm(2)=0'],end run.vavtran(itrgood)=itran; run.vavtsec(itrgood)=run.pstsecav(itran); end end % Truncate the arrays to the number of GOOD transects. ntrgood=itrgood; run.vavcen=run.vavcen(1:ntrgood); run.nbin=run.nbin(1:ntrgood); run.vavfull=run.vavfull(1:ntrgood); run.vaveast=run.vaveast(1:ntrgood); run.vavwest=run.vavwest(1:ntrgood); run.vavedge=run.vavedge(1:ntrgood); run.vavuvm=run.vavuvm(1:ntrgood); run.vavtran=run.vavtran(1:ntrgood); run.vavtsec=run.vavtsec(1:ntrgood); % Interpolate for logger values of stage and UVM velocity. run.vavhcs=interp1(run.cs20.pstseccs,run.cs20.hcs,run.vavtsec); run.vavvcs=interp1(run.cs20.pstseccs,run.cs20.vcs,run.vavtsec); plot(run.vavcen,run.vavedge,'-o'); axislim=max([max([run.vavcen,run.vavedge]), ... -min([run.vavcen,run.vavedge])]); axis([-axislim,axislim,-axislim,axislim]); hold on; plot([0,0],[-axislim,axislim],'black-.'); plot([-axislim,axislim],[0,0],'black-.'); title(['Edge Effects, Tms4, With y^1^/^6 Top-And-Bottom Correction'], ... 'fontsize',14); xlabel('velocity averaged over measured channel width (m/s)','fontsize',12); ylabel('estimated velocity for missing edges (m/s)','fontsize',12); hold off;