

REstar = [2400 1800 1200 600 0 -600 -1200 -1800];
nR = length(REstar);
tbp0_Bo200 =  [1020.12 nan 1020.85 1020.94  1020.96 1020.95 1020.98 1021.23];
tbp0_Bo1000 = [1020.13 1020.51 1020.83 1020.94  1020.96 1020.95 1020.98 1021.23];
tbp0_Bo4000 = [1020.17 1020.60 1020.74 1020.94 1020.97  1021.00  1021.06 1021.22];

BoA=[200 1000 4000];


istart_Bo200 = [11238  nan 11512 11513 9455 9479 11120 11451];
iend_Bo200 = [11618 nan   11920 11926  9866  9872 11526 11856];

% for Bo = 1000
istart_Bo1000 = [11240  11378 11494 11494 9455 7459 11059 11569];
iend_Bo1000 = [11622 11761  11902 11906  9866  7850 11467 11978];

% for Bo=4000;
istart_Bo4000 = [11262  11420  11469 11509  9517  7193  28491 12322 ];
iend_Bo4000 = [11659 11810  11870 11921  9934  7604 28907 12725];

%iBo = 1;

OV_ALL = cell(3,1);

%Bo = BoA(iBo);  % bond number

day = 31;
Area = nan*ones(8,1);

for iBo=3:3,  % for all 3 it would be 1:3

    Bo=BoA(iBo);
    eval( sprintf('tbp0 = tbp0_Bo%d;',Bo) );
    eval( sprintf('istart = istart_Bo%d;',Bo) );
    eval( sprintf('iend = iend_Bo%d;',Bo) );

for i=1:8,
    if isnan(tbp0(i))
        continue;
    end;
    
    index0=return_basilisk_index(tbp0(i));
    
    index=index0;
    [F,x,y,t]  = read_binary_file(day,REstar(i),Bo,index,'f');
    %[P,x,y,t]  = read_binary_file(day,REstar(i),Bo,index,'p');
    %[Z,x,y,t]  = read_binary_file(day,REstar(i),Bo,index0,'z');
    %    ii = find(F>0.2);  F(ii)=1;   ii=find(F<=0.21);  F(ii)=0;

    R = read_surface_extract(day,REstar(i),Bo,t);
    [Xs,Ys,Xraw,Yraw] = read_processed_surface_extract(day,REstar(i),Bo,t);
    ns = size(Xs,1);
    
     if ((i==2)&(Bo==1000)),
         ibad = [11649:11652]';
         is = [1:ns];
         igood = setdiff(is,ibad);
         Xs=Xs(igood);
         Ys=Ys(igood);
     end;
        
    figure(i),clf
    imagesc(x,y,F)
    axis xy
    hold on
    plot(Xs,Ys,'m.') %,Xraw,Yraw,'r.')

    ns = length(Xs);   S=[1:ns]/ns';    
    %    scatter(Xn,Yn,4,S,'o')
    colorbar('vert')
    %    plot(Xn,Yn,'m.')

    xmin = min(x)+4;  xmax=max(x)-1.5;
    zmin=-0.1;  zmax=0.75;
     axis([xmin  xmax  -0.1 0.7])
     plot(Xs(istart(i):iend(i)),Ys(istart(i):iend(i)),'ro')
    %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,Bo=%d; REstar=%d,  t=%.3f,  Area=%.4f',Bo, REstar(i),t,Area(i)))    
    %    print('-dpng','-r300',sprintf('Figures/Overturn_Area_REstar%d_BO1000.png',REstar(i)))

    ii = find( (Xs>xmin)&(Xs<xmax)&(Ys>zmin)&(Ys<zmax) );
    Xov = Xs(istart(i):iend(i));
    Yov = Ys(istart(i):iend(i));    
    %    Xs = Xs(ii);  Ys=Ys(ii);
    if (i==1)
        OV = struct('Xs',Xs,'Ys',Ys,'Xov',Xov,'Yov',Yov);
    else
        OV(i).Xs=Xs;  OV(i).Ys=Ys;  OV(i).Xov =Xov;  OV(i).Yov=Yov;
    end
    %    eval(sprintf('OV%d = OV;',i));
end


OVA = struct('REstar',REstar,'Bo',Bo,'tbp0',tbp0,'istart',istart,'iend',iend,'OV',OV);

OV_ALL{iBo} = OVA;

end


save overturn_analysis/OV_ALL.mat   OV_ALL REstar BoA
