function [Usig,twoA,thetas,period] = zeroups(vx,vy,zp,F)
% ZEROUPS - Zero up-crossing puv wave statistics
% [Usig,twoA,thetas,period] = zeroups(vx,vy,zp,F)
%
% This routine takes two axis velocity and a pressure height
% m/s m/s and m(H2O) and does counting-type wave statistics
% Doug Wilson, OHSU

dfifty = 0.00017;
vabs = sqrt(vx.*vx + vy.*vy);
vabsmean = mean(vabs);
vabs2mean = mean(vabs.*vabs);

vabsstd = std(vabs);

% CUMSUM of the velocity is a way to convert velocity to 
% excursion distance s = vdt.  

xpos  = cumsum(vx)./F;
ypos  = cumsum(vy)./F;
excur = sqrt(xpos.^2 + ypos.^2);

% ---------------------------------------------------------------------- %
%  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *  *
%  To get a distribution of wave heights, I find the ZERO UP CROSSINGS!
zpstd  = std(zp);		 % RMS wave height.
zpmean = mean(zp);		 % mean water depth above sensor
iof = find(zp<zpmean);  % iof is the locations of all z < mean.

len_iof = length(iof);

zptest = zp;		 		 % create a fool around vector.
zpl    = length(zptest);
zptest(iof) = ones(len_iof,1).*1000;		 % make all zp<mean = 1000.
dzp_dt = zptest(2:zpl) - zptest(1:zpl-1); 		 % differentiate by 
		 		 		 % subtracting i+1th from ith
cr = find(dzp_dt<-900);		 		 % find all really negative slopes
% (since realzp-1000 is -999.?)  this differentiates negative slopes
% due to noise in the data from negative slopes due to zero-up crossings
		 
lcros = length(cr);		 		 % number of up crossings
%		 stf     = zeros(200,1);
%        stf(cr) = ones(lcros,1);
%		 plot([1:200],zp(1:200),'*',[1:200],stf(1:200)-1.2,'-')

% 		 Now calc. all wave heights since we have indices for all the up
% 		 crossings.  Find the max and min of zp within crossing(q) and (q+1)
%		 the p is just for multiple sections of the file (which was done
% 		 because the mean zp tends to wander, and crossings is somewhat 
% 		 sensitive to zpmean.
p = 1;
for q=1:lcros-1		 
  pera(p) = cr(q+1) - cr(q);		 		 
  hite(p) = max(zp(cr(q):cr(q+1)))- min(zp(cr(q):cr(q+1)));   
  vtmax(p) = max(vabs(cr(q):cr(q+1)));
  vtmin(p) = min(vabs(cr(q):cr(q+1)));
  tvec     = detrend(excur(cr(q):cr(q+1)));
  exmax(p) = max(tvec);
  exmin(p) = min(tvec);
  twoa(p)  = exmax(p) -exmin(p);
 p = p+1;
end

% and find H1/3 by sorting hites and averaging the highest one third
ll       = length(hite); 
hsort    = sort(hite);
vmaxsort = sort(vtmax);
vminsort = sort(-vtmin);
twoasort = sort(twoa);
% persort  = sort(pera);
twothird = ll - fix((ll)/3)+1;
Honethird = mean(hsort(twothird:ll));
vmaxthird = mean(vmaxsort(twothird:ll));
vminthird = mean(vminsort(twothird:ll));
twoathird = mean(twoasort(twothird:ll));
% perathird = mean(persort(twothird:ll));
%DeltaV = vmaxthird - vminthird;
Usig   = vmaxthird;
twoA   = twoathird;
period = mean(pera)./F;
%B(2) = Usig(k);

%twoA = Usig*peraray/pi;
fw = exp(5.213*(2.5*dfifty*2./twoA).^0.194 - 5.977);
thetas = fw*Usig.^2/(32.34*dfifty);