% obtain massflux and mean return flow velocity for the LIP-experiment
%
%
%
%
inst1 = 'DRO00';
ij(1,:) = [5];
ij(2,:) = [6];
id      = [1 2 3 4 5 6 7 8 9 10];  % position of pressure sensor 
ib      = [1];  % corresponding bottom profile
xc      = [20 65 100 115 130  138  145  152  160  170]; % position of pressure sensor
zidro   = 4.1-[3 3 3.2 3.3 3.55 3.5 3.7 3.6 3.44 3.7];

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

Tp   = 5; 
for j = 1:1
   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 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
   

      % get pressure time series
   for ist=1:10
   inst1 = 'DRO00';
   if (id(ist))<10
    inst1(5:5)= num2str(id(ist))
   else
    inst1(4:5)= num2str(id(ist))
   end
    
   [t,g] = get_series(test,inst1,0);

   
   % 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.6; %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);
   
   h = 4.1 - interp1(x,z,xc(ist))   
   pause
% position of instruments with respect to mean sea level
   zp = zidro(id(ist));
   elev = h - zp;
   k = disper(2*pi*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 of incident waves (high frequency)  

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

 H1a(ist) = Hrms
 H1b(ist) = Hrms_lo

end
end



