

REstar = [2400];
nR = length(REstar);

Bo=[200 1000 2000 4000 8000];
tbp0 = [1020.2    1020.2 1020.2 1020.2 1020.2 ];
j_index = [-7 -6 -6 -2 3];


Area = nan*ones(5,1);

for i=1:length(Bo), %nR,
    index0=return_basilisk_index(tbp0(i));
    j = j_index(i);
    
    index=index0+j;
    [F,x,y,t]  = read_binary_file(9,REstar,Bo(i),index,'f');
    [P,x,y,t]  = read_binary_file(9,REstar,Bo(i),index,'p');
    %[Z,x,y,t]  = read_binary_file(9,REstar(i),Bo,index0,'z');
    %    ii = find(F>0.2);  F(ii)=1;   ii=find(F<=0.21);  F(ii)=0;

    figure,clf
    imagesc(x,y,F)
    axis xy
    hold on
    axis([min(x)+4.7  max(x)-2.2  -0.1 0.7])
    

    R = read_surface_extract(9,REstar,Bo(i),t);
    [nR,mR] = size(R);
    PC = [R(:,3:4)  0*R(:,3)];
    Rpc = pointCloud( PC );
    Rpc = pcdenoise(Rpc,'Threshold',0.001);
    [labels,nclusters]=pcsegdist(Rpc,0.01);
    ncpoints = zeros(nclusters,1);
    ii = [];
    point_limit = 3;
    if i>=4,
        point_limit = 1;
    end;
     for k=1:nclusters,
         ncpoints(k) = sum( labels==k);
         if ncpoints(k)>point_limit,
             itmp = find(labels==k);
             ii = [ii; itmp];
         end;
     end;
     %     [mxnc,imxnc]=max(ncpoints)
     %ii = find(labels == imxnc);
     nRpc = select(Rpc,ii);
    [Xn,Yn]=sort_surface_extract(nRpc.Location(:,1), nRpc.Location(:,2));     
     %    nRpc = Rpc;
    %    X = [nRpc.Location(:,1); nRpc.Location(end,1); nRpc.Location(1,1)];
    %Y = [nRpc.Location(:,2); -0.3; -0.3];
    nxn = length(Xn);   S=[1:nxn]/nxn';
    indA = [1:nxn]';
    ibad=[];
%     if i==3,
%         ibad = [11930 11931];  % REstar=1200
%     end;
    igood = setdiff(indA,ibad);
    Xn = Xn(igood);  Yn =Yn(igood);
    nxn = length(Xn);   S=[1:nxn]/nxn';    
    scatter(Xn,Yn,4,S,'o')
    colorbar('vert')
    %    plot(Xn,Yn,'m.')

        axis([min(x)+4.7  max(x)-2.2  -0.1 0.7])
    inode_start = [16252 1335 11257  1387 15867];
    inode_end = [16776 1842 11795  1894 16382];    
    %    istart=16860;   iend=17394;  %%REstar=2400
    %istart=11460;  iend=11983;  %% REstar=1800
    %istart = 11842;  iend=12457  %% REstar=1200;
    %istart = 11852;  iend=12470  %% REstar=600;    
    % istart=11903;   iend=12504;  %% REstar=-1800
        plot(Xn(inode_start(i):inode_end(i)),Yn(inode_start(i):inode_end(i)),'m.')
    Area(i) = polyarea(Xn(inode_start(i):inode_end(i)),Yn(inode_start(i):inode_end(i)))
    %    title(sprintf('pcsegdist=0.001, REstar = %d,  index=%d, t=%.3f, j=%d',REstar(i),int32(index),t,j))
    title(sprintf('pcsegdist=0.001, REstar=%d, Bo=%d, t=%.3f, Area=%.3f',REstar,Bo(i),t,Area(i)))
    print('-dpng','-r300',sprintf('Figures/Overturn_Area_RE%d_Bo%d.png',REstar,Bo(i)))
    %    else
    %  title(sprintf('REstar = %d, Bo=%d, index=%d, t=%.3f',REstar,Bo(i),int32(index),t))                
    %end;
end
