% obtain massflux and mean return flow velocity for the LIP-experiment
%
%
%
%
clear all

ij(1,:) = [1 2  3 4  7 8  10  12 15 16 17];
ij(2,:) = [2 13 4 5  7 8  10  13  4  5  6];
id      = [2  3 4 5  7 8  10  3  4  5  6];  % position of pressure sensor 
ib      = [1  2 3 4  7 8  12  12 12 18  18];  % corresponding bottom profile
xc      = [65 102 115 130 145 152 170 102 115 130 138]; % 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];


Tp   = 5; 
nij = length(ij);
for j = 1:nij
   test = '#1b0000'; 
   testb = '%1b003c';
   inst1 = 'DRO00'; 
   inst2 = 'EMF00';
 
   %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
   
   if id(j) < 10
      inst1(5:5)= num2str(id(j))
   else
      inst1(4:5)= num2str(id(j))
   end
   
   [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];
 zz = [0:.1:max(ziemf)];
 ziemf2 = [0 ziemf];
 Um2 = [0 Um];
 ie = find(ziemf2<htr);
 [pp,ss] = polyfit(ziemf2(ie),Um2(ie),3);
 uu = polyval(pp,zz);
 uu = spline(ziemf2(ie),Um2(ie),zz);
 plot(uu,zz,Um2,ziemf2,'o')
 title(num2str(xc(j)))
 pause
 Q1b(j) = trapz(zz,uu)
 U1b(j) = Q1b(j)/htr
 H1b(j) = Hrms
 Z1b(j) = h

end


