% obtain massflux and mean return flow velocity for the LIP-experiment
%
%
%
%
inst1 = 'DRO00';
inst2 = 'EMF00';
ij(1,:) = [1 2 3 4 5 6 7 8 9 10 11];
ij(2,:) = [2 3 4 5 6 7 8 9 11 12 13];
id      = [2 3 4 5 6 7 8 9  9  6  3];  % position of pressure sensor 
ib      = [1 2 3 4 5 6 7 7  7 12 12];  % corresponding bottom profile
xc      = [65 100 115 130  138  145  152  160  156  141  102]; % position of carriage
ziemf   = [0.1 0.2 0.4 0.7 1.1]; % height above the bed
zidro   = 4.1-[3 3 3.2 3.3 3.55 3.5 3.7 3.6 3.44 3.4];

test = '#1a0000';
testb = '%1a003c';

Tp   = 5; 
nij = length(ij);
for j = 2:nij
%for j = 2:2
   if ij(1,j) < 10
      test(5:5) = num2str(ij(1,j))
   else
      test(4:5) = num2str(ij(1,j))
   end
   if ij(2,j) < 10
      test(7:7) = num2str(ij(2,j))
   else
      test(6:7) = num2str(ij(2,j))
   end
   
   % get velocity time series
   
   for ist=1:5
      inst2(4:5) = num2str(10+2*ist-1)
      [t,g] = get_series(test,inst2,0);
      Um(ist) = mean(g);
   end
   
   % get pressure time series
   
   inst1(5:5)= num2str(id(j))
   [t,g] = get_series(test,inst1,0);
   
   % get bottom profile
   
   if ib(j) < 10
      testb(5:5) = num2str(ib(j))
   else
      testb(4:5) = num2str(ib(j))
   end

   [x,z] = get_profile(testb,0);
   
   % local depth
   
   h = 4.1 - interp1(x,z,xc(j))   
   
   % obtain local wave height 
   dt = t(2)-t(1)
   m1 = 2048;
   m2 = m1/2;
%   df = 1/(m1*dt);

% pressure spectrum
   ng = length(g);
   ng1 = floor(ng/10);
   ng2 = ng-ng1;
   P = detrend(g(ng1:ng2))/9.81;
   [pf, f] = crosgk(P,P,m1,1,dt,1);
   pp = pf(:,1);
   nf = length(pf);
   df = f(2)-f(1);
   
   f_high = 0.5; %cut off frequency
   ii = find(f>f_high);
   ni = length(ii); 
   f_low = 1/Tp/2; % high pass filter
   jj = find(f<f_low);
   nj = length(jj);
   

% position of instruments with respect to mean sea level
   zp = zidro(id(j));
   elev = h - zp;
   k = disper(f,h);
   
   ht = ones(nf,1)*h;
   hz = ones(nf,1)*elev;
   trans = (cosh(k.*ht)./(cosh(k.*hz))).^2;
   trans(ii,:) = zeros(size(trans(ii,:)));
   trans(jj,:) = zeros(size(trans(jj,:)));

% surface elevation spectrum
   pt = pp.*trans;
   
% wave height   

  Hrms = 4*sqrt(sum(abs(pt)).*df)/sqrt(2)
  
% estimate the mass-flux
% fit a parabola through the wet points

 htr = h-Hrms/2;
 zz = [0:.1:htr];
 ie = find(ziemf<htr);
 [pp,ss] = polyfit(ziemf(ie),Um(ie),2);
 uu = polyval(pp,zz);
 plot(uu,zz,Um,ziemf,'o')
 pause(1)
 Q1a(j) = trapz(zz,uu)
 U1a(j) = Q1a(j)/htr
 H1a(j) = Hrms
 Z1a(j) = h

end


