% function R = read_surface_extract(day,REstar,Bo,time)
%
%  reads in the raw surface extract file for a particular simulation and time
%
% format of E is [time KE_w GPE_w  diss_w   max(u)_w KE_a  GPE_a  diss_a  max(|U|)_a  mass_water mass_air

function E = read_energy_file(day,REstar,Bo)

if day==18,  %day = 9;   
    dirstart = 'RE';   %
elseif day==9 || day==31,
    dirstart = 'shoal_RE';
else
    disp('*** Bad day to read_energy_file() ***')
    E=[];
    return
end;

prefix=sprintf('/home/kghanson/basilisk/shoaling/runs/test_2022_11_%02d/%s%d_BO%d',day,dirstart,REstar,Bo);
energy_file = sprintf('%s/EnergyBudget.dat',prefix);
if exist(energy_file)==2,
    E = load( energy_file);
else
    disp('* could not find energy file')
    E = [];
end;

% still need to do some stuff with making energy dissipation have rho on it
return

%windfile1=sprintf('%s/shoal_horiz_profs',prefix1);
%windfile2=sprintf('%s/shoal_horiz_profs',prefix2);
%W1 = load(windfile1)

%TT = load( sprintf('%s/times',prefix));

% for 11_18  BO=1000
%tbreak1 = 1020.72;  tbreak2 = 1020.90;

% for 11_09  BO=1000;
%tbreak1 = 1020.51;  tbreak2 = 1020.98;

%results_file2 = sprintf('%s/results%.3f_0.dat',prefix2,tbreak2)



return
%R2 = load( results_file2);

figure(1),clf
plot(R1(:,3),R1(:,4),'k.',R2(:,3),R2(:,4),'r.')
%scatter(R1(:,3),R1(:,4),[],R1(:,8),'.')
grid
%hold on
%scatter(R2(:,3),R2(:,4),[],R2(:,8),'o')
%colorbar('vert')
axis([39 43  -0.2 0.8])
xlabel('x/h')
ylabel('z/h')
stitle = sprintf('REstar = (%d, %d) (black,red): Bo=%d: moment of impact: t1=%.3f, t2=%.3f',REstar1,REstar2,Bo1,tbreak1-1000,tbreak2-1000);
title(stitle)
%title(sprintf('ind = %d,  file = %s',ind,DD(ind).name))
%pause(0.02)


% DD = dir( sprintf('%s/*.dat',prefix));
% 
% 
% for ind=350:450,
% R = load( sprintf('%s/%s',DD(ind).folder,DD(ind).name));
% 
% scatter(R(:,3),R(:,4),[],R(:,8),'.')
% grid
% colorbar('vert')
% axis([0 60 -0.2 0.8])
% title(sprintf('ind = %d,  file = %s',ind,DD(ind).name))
% pause(0.02)
% end;
