
Bo=4000;
REstarA = [2400 -1800];

for k=1:2,
    REstar = REstarA(k);
for i = 300:300,
    index = i;
    [~,time] = return_basilisk_time(index);
    
% get the surface extract file
[Xs,Ys,Ro] = read_processed_surface_extract(9,REstar,Bo,time);

%R = [X1 X2 Xg Yg U W vort P]

R = Ro(2:end,:);  % remove start and end

P = read_Poffset_file(9,REstar,Bo,time);  %[Xmid Ymid  P]
Pmid = 0.5*( P(1:end-1,:) + P(2:end,:));
Psurf = 0.5*( P(1:end-1,3) + P(2:end,3));

ix1=xindex(Xs,Pmid(1,1));
ix2=xindex(Xs,Pmid(end,1));

% if Xs(ix1)>Pmid(1,1),
%     ix1 = ix1-1;
% end;
% 
% if Xs(ix2)>P(end,1),
%     ix2 = ix2-1;
% end;

Rn = R(ix1:ix2,:);
Rn(:,8) = Psurf;

size(Rn)
size(Psurf)

[~,imax] = max(Rn(:,4)); % location of the highest point
ds = sqrt( (diff(Rn(:,3))).^2 + (diff(Rn(:,4))).^2 );

S = [0; cumsum(ds)];
S = S - S(imax);

Umag = 0.5*sqrt( Rn(:,5).^2 + Rn(:,6).^2);

rho_air= 0.000980392156862745; % really the ratio
Phyd = -rho_air*Rn(:,4);
Pair = Rn(:,8)+Phyd;  % puts the hydrostatic part back in

figure,clf, orient portrait
subplot(2,1,1)
plot(Rn(:,3),Rn(:,4),'*','markersize',2),grid
hold on
%plot(P(:,1),P(:,2),'+') 
plot(Pmid(:,1),Pmid(:,2),'o','markersize',2)
hold off
stitle = sprintf('RE*= %d, dt=%.2f',REstar,time);
title(stitle)
xlabel('$x/h_0$','interpreter','latex')
ylabel('$\eta/h_0$','interpreter','latex')
axis([33 38.5 0 0.8])

subplot(2,1,2)
plot(S,Pair,'b',S,Umag,'r',S,Rn(:,4),'g'),grid
xlabel('$S$','interpreter','latex')
ylabel('dynamic BC terms')
legend('$P_a$','$(1/2)|\vec{u}|^2$','$\eta$','interpreter','latex')
axis([-4 3 -0.05 0.8])

if (k==1),
Rn2 = Rn;
Pair2 = Pair;
S2 = S;
end;

end;
end;

figure,clf, orient tall
subplot(2,1,1)
plot(Rn(:,3),Rn(:,4),'*','markersize',2),grid
hold on
%plot(P(:,1),P(:,2),'+') 
plot(Rn2(:,3),Rn2(:,4),'ro','markersize',2)
hold off
stitle = sprintf('dt=%.2f',time);
title(stitle)
legend('RE*=-1800','RE*=2400','interpreter','latex')
xlabel('$x/h_0$','interpreter','latex')
ylabel('$\eta/h_0$','interpreter','latex')

subplot(2,1,2)
plot(S,Pair,'b',S2,Pair2,'r'),grid
xlabel('$S$','interpreter','latex')
ylabel('$P_\mathrm{air}$','interpreter','latex')
%legend('$P_a$','$(1/2)|\vec{u}|^2$','$\eta$','interpreter','latex')

figure(1),print -dpdf ../Figures/RE2400_DynBC.pdf
figure(2),print -dpdf ../Figures/RE-1800_DynBC.pdf
figure(3),print -dpdf ../Figures/Both_DynBC.pdf
