

REstar = [2400 1800 1200 600 0 -600 -1200 -1800];
nR = length(REstar);
BoA=[200 1000 4000];


Void_Area = nan*ones(8,3);
Xbp = nan*ones(8,3);
Ybp = nan*ones(8,3);
Xov_back = nan*ones(8,3);
Yov_back = nan*ones(8,3);
Xov_top = nan*ones(8,3);
Yov_top = nan*ones(8,3);
Lw_back = nan*ones(8,3);
Ytop = nan*ones(8,3);
Ytop_lip = nan*ones(8,3);
Area_Lip = nan*ones(8,3);
theta_overturn = nan*ones(8,3);
aspect_overturn = nan*ones(8,3);

theta_fit = nan*ones(8,3);
L_fit = nan*ones(8,3);
W_fit = nan*ones(8,3);

day = 9;

load(sprintf('OV_ALL_%d.mat',day));


tbp0A = nan*ones(8,3);
tbp0A(:,1) = OV_ALL{1}.tbp0';
tbp0A(:,2) = OV_ALL{2}.tbp0';
tbp0A(:,3) = OV_ALL{3}.tbp0';


for iBo = 1:3,
    OVA=OV_ALL{iBo};
    Bo = BoA(iBo);
    figure(iBo),clf,orient tall
   for i=1:nR,
        if isnan(OVA.tbp0(i)),
            continue
        end;
    index0=return_basilisk_index(OVA.tbp0(i));
    
    index=index0;
    [F,x,y,t]  = read_binary_file(9,REstar(i),Bo,index,'f');
    [P,x,y,t]  = read_binary_file(9,REstar(i),Bo,index,'p');

    subplot(4,2,i)
    imagesc(x,y,F)
    axis xy
    hold on
    plot(OVA.OV(i).Xs,OVA.OV(i).Ys,'m.') %,Xraw,Yraw,'r.')
                                         %colorbar('vert')

     xmin = min(x)+0.5;  xmax=max(x)-0.4;
     zmin=-0.1;  zmax=1;
     axis([xmin  xmax  zmin zmax])
     plot(OVA.OV(i).Xov, OVA.OV(i).Yov,'ro')

     % 1) calculate area and Xbp
     Void_Area(i,iBo) = polyarea(OVA.OV(i).Xov, OVA.OV(i).Yov);
     [Xbp0,imax] = max(OVA.OV(i).Xov);
     Xbp(i,iBo)=Xbp0;   Ybp(i,iBo)=OVA.OV(i).Yov(imax);

     
     % 2) estimate (X,Y)ov_back, Ytop and Area_Lip
     [Xback0,iback] = min(OVA.OV(i).Xov);
     Xov_back(i,iBo)=Xback0;   Yov_back(i,iBo)=OVA.OV(i).Yov(iback);
     [dmin, ibacks] = min( (OVA.OV(i).Xs - Xov_back(i,iBo)).^2 +  (OVA.OV(i).Ys - Yov_back(i,iBo)).^2 );

     % 2b top of wave
     iir = find( abs( OVA.OV(i).Xs- Xback0)<4);
     [Ytop0,itopr] = max(OVA.OV(i).Ys(iir));
     Ytop(i,iBo) = Ytop0;
     itop = iir(itopr);
     Xtop(i,iBo) = OVA.OV(i).Xs(itop);

     % 2c : find where back half of wave hits cutoff
     Ycut1 = 0.22;  Ycut2=0.2;
     iw = find( OVA.OV(i).Xs<Xtop(i,iBo));
     iw2 = find( (OVA.OV(i).Ys(iw) < Ycut1)&(OVA.OV(i).Ys(iw) > (Ycut1-0.05)));
     imx = max(iw2);
     Ytail1 = OVA.OV(i).Ys(iw(imx));
     Xtail1 = OVA.OV(i).Xs(iw(imx));                  
     Lw_back(i,iBo) = Xtop(i,iBo)-Xtail1;
     
     %3) area of the lip 
     ii = find( abs(OVA.OV(i).Xs-Xov_back(i,iBo))<0.005);
     if isempty(ii)
         disp('** WHOA badness')
         keyboard
     end;
     
     [Ytop_lip(i,iBo),imax] = max( OVA.OV(i).Ys(ii));     
     ilip_start = ii(imax);
     ilip_end = ibacks;
     Area_Lip(i,iBo) = polyarea(OVA.OV(i).Xs(ilip_start:ilip_end), OVA.OV(i).Ys(ilip_start:ilip_end));

     %4) LH82 function fit for angle and aspect ratio
     [Xn, Zn, params_fit, Xfinal, Zfinal, Rfit] = iterate_LH82_function_fit(OVA.OV(i).Xov,OVA.OV(i).Yov, Xback0);
     theta_fit(i,iBo) = params_fit(3);
     L_fit(i,iBo) = params_fit(1)*6;
     W_fit(i,iBo) = params_fit(2)*2.1773;
     overturn_orientation
     theta_overturn(i,iBo)=theta0;
     aspect_overturn(i,iBo)=aspect0;
     
     plot(OVA.OV(i).Xs(ilip_start:ilip_end), OVA.OV(i).Ys(ilip_start:ilip_end),'k^','markersize',2,'markerfacecolor','k')

     
     plot(Xbp(i,iBo),Ybp(i,iBo),'gd','markersize',6,'markerfacecolor','g')
     plot(Xov_back(i,iBo),Yov_back(i,iBo),'gd','markersize',6,'markerfacecolor','g')
     plot(Xtop(i,iBo),Ytop(i,iBo),'gd','markersize',6,'markerfacecolor','g')     
                  plot(Xov_back(i,iBo),Ytop_lip(i,iBo),'cd','markerfacecolor','c')
     plot(Xtail1,Ytail1,'cd','markerfacecolor','c')                    
     title(sprintf('Bo=%d; REstar=%d,  t=%.2f',Bo,REstar(i),t))    

%      if (i==7)&(iBo==3),
%          keyboard
%      end;
   end
         print('-dpng','-r300',sprintf('../Figures/Overturn_Area_BO%d_AllREstar.png',Bo))
end;


%%%%

load ../mat/wind_1109.mat

UC7 = Umean(:,:,7)/C
UC3 = Umean(:,:,3)/C
UC7 = UC7';
UC3 = UC3';

figure(10),clf,orient tall

subplot(3,2,1)
plot(REstar,Xbp(:,3),'*'),grid
%xlabel('Uw/C at y=3')
%xlabel('RE^*')
ylabel('Xbp')

% subplot(3,2,2)
% plot(REstar,tbp0A,'*'),grid
% %xlabel('Uw/C at y=3')
% %xlabel('RE^*')
% ylabel('time bp')

subplot(3,2,2)
plot(REstar,Ytop(:,3),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('y Top of the wave')

subplot(3,2,3)
plot(REstar,Void_Area(:,3)./(Ytop(:,3).^2),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('Void Area/H^2')

subplot(3,2,4)
plot(REstar,W_fit(:,3)./L_fit(:,3),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('Aspect Ratio W/L')

subplot(3,2,5)
plot(REstar,Area_Lip(:,3)./(Ytop(:,3).^2),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('Lip Area/H^2')
xlabel('RE^*')

subplot(3,2,6)
plot(REstar,theta_fit(:,3),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('angle of overturn (deg)')
xlabel('RE^*')

print('-dpdf','../Figures/overturn_geometry_REstar.pdf')

figure(4),clf
subplot(3,2,1)
plot(REstar,L_fit,'*')
ylabel('L')

subplot(3,2,2)
plot(REstar,W_fit,'*')
ylabel('W')

UC=UC3;

figure(11),clf,orient tall

subplot(3,2,1)
plot(UC3,Xbp,'*'),grid
%xlabel('Uw/C at y=3')
%xlabel('RE^*')
ylabel('Xbp')

% subplot(3,2,2)
% plot(UC3,tbp0A,'*'),grid
% %xlabel('Uw/C at y=3')
% %xlabel('RE^*')
% ylabel('time bp')

subplot(3,2,2)
plot(UC3,Ytop,'*'),grid
%xlabel('Uw/C at y=3')
ylabel('y Top of the wave')

subplot(3,2,3)
plot(UC3,Void_Area./(Ytop.^2),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('Void Area/H^2')

subplot(3,2,4)
plot(UC3,W_fit./L_fit,'*'),grid
%xlabel('Uw/C at y=3')
ylabel('Aspect Ratio W/L')

subplot(3,2,5)
plot(UC3,Area_Lip./(Ytop.^2),'*'),grid
%xlabel('Uw/C at y=3')
ylabel('Lip Area/H^2')
xlabel('U/C @ y=3')

subplot(3,2,6)
plot(UC3,theta_fit,'*'),grid
%xlabel('Uw/C at y=3')
ylabel('angle of overturn (deg)')
xlabel('U/C @ y=3')

print('-dpdf','../Figures/overturn_geometry_UC3.pdf')

save OVERTURN_ANALYSIS.mat
