

REstar = [2400 0  -1800];
nR = length(REstar);
BoA=[ 4000 ];
nBo = length(BoA);

ENG = cell(nR);
NN = nan*ones(nR);
MM = nan*ones(nR);
tmax = nan*ones(nR);

for iBo=1:nBo,
    for i=1:nR,
        ENG{i} = read_energy_file(9,REstar(i),BoA(iBo));
        [n,m]=size(ENG{i});
        NN(i) = n;  MM(i)=m;
        tmax(i) = ENG{i}(end,1);
    end;
end;

Nmax = min(NN(:));
tall = ENG{1}(1:Nmax,1);
ind2 = find(tall==1001.5);

TKEw = nan*ones(nR,Nmax);
TKEa = nan*ones(nR,Nmax);
PEw = nan*ones(nR,Nmax);
PEa = nan*ones(nR,Nmax);
maxSw = nan*ones(nR,Nmax);
maxSa = nan*ones(nR,Nmax);
MASSw = nan*ones(nR,Nmax);
MASSa = nan*ones(nR,Nmax);

for i=1:nR,
        TKEw(i,:) = ENG{i}(1:Nmax,2);
        TKEa(i,:) = ENG{i}(1:Nmax,6);
        PEw(i,:) = ENG{i}(1:Nmax,3);
        PEw(i,:) = PEw(i,:)-PEw(i,ind2);
        PEa(i,:) = ENG{i}(1:Nmax,7);
        PEa(i,:) = PEa(i,:)-PEa(i,ind2);
        maxSw(i,:) = ENG{i}(1:Nmax,5);
        maxSa(i,:) = ENG{i}(1:Nmax,9);
        %        MASSw(i,:) = ENG{i}(1:Nmax,10);
        %MASSa(i,:) = ENG{i}(1:Nmax,11);
end;

TKEw(1,2,:)=nan;  TKEa(1,2,:)=nan;  PEw(1,2,:)=nan;  PEa(1,2,:)=nan;
Ew = TKEw+PEw;
Ea = TKEa+PEa;

ichoose = [1 nR];
ichoose = [1:nR];

xlim0 = [1002 1021.2];
figure(1),clf,orient tall
subplot(5,2,1)
plot(tall,squeeze(Ew(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('water E = TKE+dPE')
%legend('2400','1800','1200','600','0','-600','-1200','-1800','location','northwest')

subplot(5,2,2)
plot(tall,squeeze(Ea(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('air E = TKE+dPE')

subplot(5,2,3)
plot(tall,squeeze(TKEw(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('water TKE')
%legend('2400','1800','1200','600','0','-600','-1200','-1800','location','northwest')

subplot(5,2,4)
plot(tall,squeeze(TKEa(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('air TKE')

subplot(5,2,5)
plot(tall,squeeze(PEw(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('water dPE')

subplot(5,2,6)
plot(tall,squeeze(PEa(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('air dPE')

subplot(5,2,7)
plot(tall,squeeze(maxSw(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('water max speed')

subplot(5,2,8)
plot(tall,squeeze(maxSa(ichoose,:)))
xlim(xlim0)
xlabel('time')
ylabel('air max speed')

% subplot(5,2,9)
% plot(tall,squeeze(MASSw(ichoose,:)))
% xlim(xlim0)
% xlabel('time')
% ylabel('water MASS')
% 
% subplot(5,2,10)
% plot(tall,squeeze(MASSa(ichoose,:)))
% xlim(xlim0)
% xlabel('time')
% ylabel('air MASS')


%print -dpdf energy_Bo1000.pdf
clear ENG
save EnergyStats4000.mat
