% This script calculates average velocities for comparison with results
% from hydrodynamical modeling.  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, October 29, 1999.

   itran1=1;
   itran2=run.ntran;
   ntran=itran2-itran1+1;
   sirun.vavfull=zeros(ntran,1); % Average velocity over the full transect.
   sirun.vavcen=zeros(ntran,1); % Average over the center of the transect.
   sirun.vaveast=zeros(ntran,1); % Average velocity for east section.
   sirun.vavwest=zeros(ntran,1); % Average velocity for west section.
   sirun.vavedge=zeros(ntran,1); % Average for east and west sections.
   
   sirun.vavtran=zeros(ntran,1); % Transect number for this good transect.
   sirun.nbin=zeros(ntran,1); % Number of bins used to calculate vavfull.
   sirun.vavtsec=zeros(ntran,1); % Time for this good transect.
   itrgood=0;

% Loop over all transects of the calibration.  Only good transects,
% as identified in the array run.goodtran, will be used in the
% averaging.
   for itran=itran1:itran2
   nens=run.nens(itran);
   if run.goodtran(itran)>=1 % Test scan code for good transects.
      itrgood=itrgood+1;
% In the following, the second component is a bin counter.
      vfull=zeros(2,1); % velocity in center segment.
      veast=zeros(2,1); % velocity in east segment.
      vwest=zeros(2,1); % velocity in west segment.
      vedge=zeros(2,1); % velocity in west and east segments.
      vcen=zeros(2,1); % velocity in center segment.
% dmax is the length of this transect.  dwest is the end of the part
% to use for vwest, similarly for deast.
      dmax=max(run.dist(1:nens,itran));
      dwest=dmax/4;
      deast=dmax*3/4;
      for iens=1:nens
      if run.vflag(iens,itran)>=0 % exclude bad data and empty ensembles.
% Add on the contribution to the sums from a single ensemble.  vpav is
% the average over the water column of the component of the water
% velocity projected onto the channel centerline.  WNmax is the number
% of vertical bins.
	 vfull(1)=vfull(1)+run.vpav(iens,itran)*run.WNmax(iens,itran);
	 vfull(2)=vfull(2)+run.WNmax(iens,itran);
	 if run.dist(iens,itran) < dwest
	    vwest(1)=vwest(1)+run.vpav(iens,itran)*run.WNmax(iens,itran);
	    vwest(2)=vwest(2)+run.WNmax(iens,itran);
	 end 
         if run.dist(iens,itran) > deast 
	    veast(1)=veast(1)+run.vpav(iens,itran)*run.WNmax(iens,itran);
	    veast(2)=veast(2)+run.WNmax(iens,itran);
	 end 
      end
      end
      
      vedge=veast+vwest;
      vcen=vfull-vedge;
      if vfull(1)==0, pause,end
%%%%%% Carry out averages.
      if veast(2)~=0, sirun.vaveast(itrgood)=veast(1)/veast(2);, ...
	 else sirun.vaveast(itrgood)=sirun.vaveast(itrgood-1);, ...
	 ['transect ',int2str(itran),': veast(2)=0'],end
      if vwest(2)~=0, sirun.vavwest(itrgood)=vwest(1)/vwest(2);, ...
	 else sirun.vavwest(itrgood)=sirun.vavwest(itrgood-1);, ...
	 ['transect ',int2str(itran),': vwest(2)=0'],end
      if vfull(2)~=0, sirun.vavfull(itrgood)=vfull(1)/vfull(2);, ...
	 else sirun.vavfull(itrgood)=sirun.vavfull(itrgood-1);, ...
	 ['transect ',int2str(itran),': vfull(2)=0'],end
      if vedge(2)~=0, sirun.vavedge(itrgood)=vedge(1)/vedge(2);, ...
	 else sirun.vavedge(itrgood)=sirun.vavedge(itrgood-1);, ...
	 ['transect ',int2str(itran),': vedge(2)=0'],end
      if vcen(2)~=0, sirun.vavcen(itrgood)=vcen(1)/vcen(2);, ...
	 else sirun.vavcen(itrgood)=sirun.vavcen(itrgood-1);, ...
	 ['transect ',int2str(itran),': vcen(2)=0'],end
      sirun.nbin(itrgood)=vfull(2);
      sirun.vavtran(itrgood)=itran;
      sirun.vavtsec(itrgood)=run.pstsecav(itran);
   end
   end

   ntrgood=itrgood;
   sirun.ntrgood=ntrgood;
% Truncate the arrays to the number of GOOD transects.
   sirun.vavfull=sirun.vavfull(1:ntrgood);
   sirun.vaveast=sirun.vaveast(1:ntrgood);
   sirun.vavwest=sirun.vavwest(1:ntrgood);
   sirun.vavedge=sirun.vavedge(1:ntrgood);
   sirun.vavcen=sirun.vavcen(1:ntrgood);
   sirun.nbin=sirun.nbin(1:ntrgood);
   sirun.vavtran=sirun.vavtran(1:ntrgood);
   sirun.vavtsec=sirun.vavtsec(1:ntrgood);

   plot(sirun.vavcen,sirun.vavedge,'-o');
   axislim=max([max([sirun.vavcen,sirun.vavedge]), ...
      -min([sirun.vavcen,sirun.vavedge])]);
   axis([-axislim,axislim,-axislim,axislim]);
   hold on;
   plot([0,0],[-axislim,axislim],'black-.');
   plot([-axislim,axislim],[0,0],'black-.');
   title(['Edge Effects, ',run.par.infile(length(run.par.infile)-11:...
      length(run.par.infile)-8)],'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;