% This function carries out calculations from raw data in the
% structure ens.  Here only quantities are calculated which 
% cannot be calculated later from information stored in the run structure.
function ens = caladens(ens)
   thetap= ens.thetap; % Angle from magnetic north to channel axis.
   theta0=thetap*pi/180.;
   stheta0=sin(theta0);
   ctheta0=cos(theta0);
   ens.vflag=0; % velocity error flag; 0 --> no error.
   ens.vw=zeros(5,ens.WN); % water velocity with respect to the bottom.
   % The fifth velocity component is for something (velocity magnitude,
   % velocity projected on a particular direction) for later use in
   % the analysis.
   for ibin=1:ens.WN
      ens.vw(1:3,ibin)=ens.vtx(1:3,ibin)-ens.BTV(1:3);
      ens.vw(4,ibin)=sqrt(ens.vtx(4,ibin).^2+ens.BTV(4).^2)* ...
	 sign(ens.vtx(4,ibin));
      ens.vw(5,ibin)=ens.vw(1,ibin)*stheta0+ ...
	 ens.vw(2,ibin)*ctheta0; % Projected horizontal velocity.
   end
% Calculate the depth array.
   ens.depth=([1:ens.WN]-1)*ens.WS + ens.dis1 + ens.dtrans;
% Determine how far down to keep the data.
   daver = mean(ens.BTR);
   ens.dbot = daver+ens.dtrans; % Use this for the distance to 
				       % the bottom.
   ens.dlimit = ens.dbot*cos(20.*pi/180.)-ens.WS;
% Find the largest bin number to use.
   ibin=1;
   while ens.depth(ibin) < ens.dlimit, ibin=ibin+1;, end
   ens.WNmax=ibin-1;
   ens.vw(1:5,ens.WNmax+1:ens.WN)=0; % Zero velocities beyond bin WNmax.
   ens=fixbadv(ens); % Deal with bad velocity measurements.  Patch up
		     % single bins, reject bad bottom velocities.

   % Calculate average velocity error.
   if ens.WNmax>0, ens.verav=mean(ens.vw(4,1:ens.WNmax));, end;
   if ens.WNmax>0, ens.verrms=norm(ens.vw(4,1:ens.WNmax))/ ...
      sqrt(ens.WNmax);, end

% Calculate normalized velocities; scale so that vE has an average value
% of +/- 1.
%   maxbin=max(1,ens.WNmax); 
%   ens.vnorm(1:5,1:maxbin)=ens.vw(1:5,1:maxbin)/abs(mean( ...
%      ens.vw(2,1:maxbin)));
%   ens.dnorm(1:maxbin)=ens.depth(1:maxbin)/ens.dbot;
% Calculate a velocity function following the y^(1/6) law, with average
% value equal to 1.
%   ens.d1over6=(0:100)/100.;
%   ens.v1over6=7./6.*(1.-ens.d1over6).^(0.16666666);

% Calculate the ensemble time in sec (at end of ensemble).
   tsec1=((ens.RTCdy*24.+ens.RTChr)*60.+ens.RTCmn)*60.+ ...
      ens.RTCsc + ens.RTChd/100.;
   if ens.iens == 1 % special case.
      ens.tsec=tsec1;
      ens.xens=[0,0,0]';
   end
   ens.etsec=tsec1-ens.tsec; % Ensemble elapsed time.
   ens.tsec=tsec1; %NOW we can update the time.
   ens.delx=ens.BTV(1:3)*ens.etsec; % Displacement during ensemble.
   ens.xens=ens.xens+ens.delx; % Accumulate displacement.


% ens=fixbadv(ens) does something with velocities flagged as bad.
% Individual bad bins are populated by interpolation or (constant)
% extrapolation, even if there is only one good bin in the ensemble.
% If all bins are bad or if the bottom velocity is bad, the measurements
% are unchanged but the ensemble is flagged as bad.
% ens.vflag=0   good bin
%           +n	n bad bins fixed
%          =-m  bad ensemble; -1=bad bottom velocity measurement
%			      -2=all bins bad
%                             -4=no usable bins (WNmax=0)
function ens=fixbadv(ens);
   if ens.BTV(4)==-32.768
      ['Bad bottom velocity, trans & ens = ',int2str(ens.itran), ...
	 ', ',int2str(ens.iens)]
      ens.vflag=-1;
   elseif ens.WNmax<=0
      ens.vflag=-4;
   else
      nbin=ens.WNmax;
      ig1=0; % First good bin.
      ig2=999; % Second good bin.
      nbad=0; % Bad bin counter.
      bad=0; % Flag for bad bin.
      for ibin=1:nbin
	 % if abs(ens.vtx(4,ibin))>.2,ens.vtx(4,ibin), end
	 if ens.vtx(4,ibin)==-32.768
	    nbad=nbad+1;
	    if bad==0
	       % start of bad section
   	       bad=1; % flag for "shit coming down"
	       if ibin>1, ig1=ibin-1;, end
	    end
	 elseif bad==1
	    % end of bad section
	    bad=0;
	    ig2=ibin;
	    ens=vinterp(ig1,ig2,ens); % interpolate (or extrapolate)
	    ig2=999; % Set up for next bad section.
	 end
      end
      if bad==1
	 if ig1>0
	    ens=vinterp(ig1,ig2,ens);
	    ens.vflag=nbad;
	 else
	    ens.vflag=-2; % all bins are bad
	 end
      else
	 ens.vflag=nbad;
      end
   end

% ens=vinterp(ig1,ig2,ens) fills in values for bad velocities.
% if ig1 and ig2 are valid bin numbers, linear interpolation is used.
% if either a final bin or an initial bin is missing,
function ens=vinterp(ig1,ig2,ens);
%   ['Bad section: itran, iens = ',int2str(ens.itran), ...
%      ', ',int2str(ens.iens),', ig1, ig2 = ',int2str(ig1), ...
%      ', ',int2str(ig2)]
   if (ig1>0)&(ig2<999) % interpolation
      for ibin=ig1+1:ig2-1
	 a=(ig2-ibin)/(ig2-ig1);
	 b=(ibin-ig1)/(ig2-ig1);
	 ens.vw(:,ibin)=a*ens.vw(:,ig1)+b*ens.vw(:,ig2);
      end
   elseif ig1>0  % extrapolation down to bottom
      for ibin=ig1+1:ens.WNmax
	 ens.vw(:,ibin)=ens.vw(:,ig1);
      end
   elseif ig2<999 % extrapolation up to surface
      for ibin=1:ig2-1
	 ens.vw(:,ibin)=ens.vw(:,ig2);
      end
   else 'ERROR:  bad call to function vinterp.', end