% Plot number of particles and total biomass in the simulation
%


% close all
clear all
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
inFold = '/home/martin/opendrift/PelagiaMFS09/';
R(1).subFold = 'C2k8_1-1/';
R(1).name = 'C7a';
R(2).subFold = 'C2k8b_1-1/';
R(2).name = 'C7b';
R(3).subFold = 'C2k8c_1-1/';
R(3).name = 'C7c';
fPrefix = 'Pelagia_Med_fwd_';
saveIMG = 1;
IMGfold = 'IMG/total_biomass/';
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

mkdir(IMGfold)


for r=1:length(R)
    inFiles = dir([inFold,R(r).subFold,fPrefix,'*.nc']);
    n = 1;
    
    for nf = 1:length(inFiles)
        inFile = [inFiles(nf).folder,'/',inFiles(nf).name];
        nc = netcdf(inFile);
        nctime = nc{'time'}(:);  % seconds since 1970-01-01 00:00:00
        biomass = nc{'biomass'}(:);  % Actually biomass of each particle
        close(nc)
        if nf==1
            nini = size(biomass,1); % intial number of particles
            mini = sum(biomass(:,1)); % intial total biomass
        end
        biomass(biomass > 1e15) = 0; % remove anomalous particles from the count
        totalMass = sum(biomass)'; % total biomass
        totalN = size(biomass,1) * ones(size(biomass,2),1); % number of particles in the simulation - doesn't change in a year
        dtime = datetime(1970,1,1) + seconds(nctime);
        subTT = timetable(dtime, totalN, totalMass);
        if nf == 1
            R(r).TT = subTT;
        else
            R(r).TT = [R(r).TT; subTT];
        end
    end
    lnames{r} = R(r).name;
    
end


hFig = figure;
set(hFig, 'Position', [10 10 800 400]); % fixed figure dimensions
set(hFig,'Color','w')
clf
hold on
for r=1:length(R)
    plot(R(r).TT.dtime, R(r).TT.totalMass / mini)
end
hold off
legend(lnames)
%     title(subFold,'Interpreter','none')
% title('Total biomass')
% ylim([0 2000*10000*100])

if saveIMG
    % save figure
    fOut = [IMGfold,'Total_biomass_multi'];
    export_fig(fOut, '-png','-pdf');
end

for r=1:length(R)
figure
plot(R(r).TT.dtime, R(r).TT.totalN)
title(['Number of particles - ',R(r).subFold],'Interpreter','none')

% if saveIMG
%     % save figure
%     fOut = [IMGfold,'Num_partricles_',strip(subFold,'/')];
%     export_fig(fOut, '-png','-pdf');
% end
end

