
REstarA = [2400    -1800];
nR = length(REstarA);

Bo=4000;

    lfile = sprintf('Collected_Surface_Extract_RE%d_Bo%d',REstarA(1),Bo);
    load(lfile)

figure(1),clf

nt = length(time);
Xmax = nan*ones(nt,nR);
Ymax = nan*ones(nt,nR);
Xmax2 = nan*ones(nt,nR);
Ymax2 = nan*ones(nt,nR);
Xmaxslope  = nan*ones(nt,nR);
max_slope  = nan*ones(nt,nR);
Xminslope  = nan*ones(nt,nR);
min_slope  = nan*ones(nt,nR);

Nxbar = nan*ones(nt,nR);
Nxbar_denoise = nan*ones(nt,nR);
istart = 200; % time=0.5
iend = 358; % t=18.3
% index 280  => t=14.00
% index 390  => t=18.30
for k=1:nR,
    REstar = REstarA(k);
    lfile = sprintf('Collected_Surface_Extract_RE%d_Bo%d',REstar,Bo);
    load(lfile)

  for i=istart:iend, %10:390, %280:280, %:nt,
  
      A = SURF{i};
      [ymax,imax]=max(A.Ys);
      Ymax(i,k)=ymax;  Xmax(i,k)=A.Xs(imax);
      ii = find((A.Xs> Xmax(i,k)-6)&(A.Xs<Xmax(i,k)+6));
      
      [fit1,G,O] = fit( A.Xs(ii), A.Ys(ii),'smoothingspline','SmoothingParam',0.99);
      YY=feval(fit1,A.Xs(ii));
      [dYdX,d2YdX2] = differentiate(fit1,A.Xs(ii));

      [Ymax2(i,k),imax2]=max(YY);
      Xmax2(i,k) = A.Xs(ii(imax2));
      
       [mxslope,imaxs] = max(dYdX);
       [minslope,imins] = min(dYdX);      
       max_slope(i,k)= mxslope;
       min_slope(i,k)= minslope;
       Xmaxslope(i,k) = Xbar2(imaxs);
       Xminslope(i,k) = Xbar2(imins);

     
      subplot(3,1,1)
      plot(A.Xs,A.Ys,'r.',A.Xs(ii),YY,'-'),grid
      title(sprintf('RE=%d, Bo=%d, doing t=%.2f',REstar,Bo,time(i)))
      axis([0 50 -0.2 0.8])
      subplot(3,1,2)
      plot(A.Xs(ii),dYdX,'-'),grid
      axis([0 50 -2 0.75])      
      xlabel('x (m)')
      
      Nxbar(i,k) = length(A.Xbar);
      PC1 = pointCloud([ A.Xbar A.slope A.Xbar*0]);
      PC2 = pcdenoise(PC1,'Threshold',0.08,'NumNeighbors',4);
      Nxbar_denoise(i,k) = length(PC2.Location(:,1));
      %min(PC2.Location(:,2))
      pause(0.01)
      Xbar2 = PC2.Location(:,1);
      Slope2 = PC2.Location(:,2);
%       ii = find((Xbar2> Xmax(i,k)-2)&(Xbar2<Xmax(i,k)+2));
%       [mxslope,imaxs] = max(Slope2(ii));
%       [minslope,imins] = min(Slope2(ii));      
%       max_slope(i,k)= mxslope;
%       min_slope(i,k)= minslope;
%       Xmaxslope(i,k) = Xbar2(imaxs);
%       Xminslope(i,k) = Xbar2(imins);
            subplot(3,1,3),
      plot(PC2.Location(:,1),PC2.Location(:,2),'.'),grid on
      axis([0 50 -2.5 0.75])            
      ylabel('slope denoise')
      xlabel('x (m)')

  end;

end;

clear SURF A
save SurfaceExtract4000.mat

