% Calculate number of individuals, using saved env. values and IBM
% equations.Growth equation - Michaelis–Menten kinetics aka Monod eq. 
% The permanent population areas is a polygon (shape file).
%

clear all
% close all

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% Particle trajectories file
inFile = '/home/martin/opendrift/PelagiaMFS09/baserun/Pelagia_baserun_4y_10k_2000.nc';
%
% Permanent population polygon
shpFile = '/home/martin/opendrift/PelagiaMFS09/mygeodata/Pelagia_polygon_1-polygon.shp';
%
HShallow = 200; % depth limit in m for additional mortality due to shallow water
%
% Parameter step and range
MPR.k1 = 0.001:0.002:0.2; % Max specific growth
MPR.k2 = 0.1:0.05:7; % substrate affinity / half saturation constant
MPR.k3 = -0.07:0.002:0.15; % natural moratlity & predation;  0.12
%
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

nc = netcdf(inFile);
nctime = nc{'time'}(:);  % seconds
lons = nc{'lon'}(:);
lats = nc{'lat'}(:);
pp = nc{'p_production'}(:);
number_ind = nc{'number_ind'}(:);
nearshore_mask = nc{'nearshore_mask'}(:);  % used for additional mortality due to beaching
status = nc{'status'}(:);
h = nc{'sea_floor_depth_below_sea_level'}(:);
close(nc)

dt = (nctime(2)-nctime(1)) / (3600*24) ; % output time step seconds to days

%
% Keep just the particles which are still in the starting area at the end of
% the run
%
%  Read the shape of the polygon
S = shaperead(shpFile); 
inArea = inpolygon(lons(:,end),lats(:,end),S.X,S.Y);
PP = pp(inArea,:);
H = h(inArea,:);
NS = nearshore_mask(inArea,:);
status = status(inArea,:);
Nopendrift = number_ind(inArea,end);

figure
hold on
plot(S.X,S.Y)
plot(lons(:,1), lats(:,1),'rx')
hold off
title('Starting positions all')

figure
hold on
plot(S.X,S.Y)
plot(lons(inArea,1), lats(inArea,1),'rx')
hold off
title('Surviving particles - initial positions')

figure
hold on
plot(S.X,S.Y)
plot(lons(inArea,end), lats(inArea,end),'rx')
hold off
title('Surviving particles - final positions')

figure
hold on
plot(S.X,S.Y)
plot(lons(inArea,:), lats(inArea,:),'.','Color',[0.5 0.5 0.5])
plot(lons(inArea,1), lats(inArea,1),'gx')
plot(lons(inArea,end), lats(inArea,end),'rx')
hold off
title('Surviving particles paths')


Nini = sum(number_ind(:,1));
disp(['Number of trajectories: ',num2str(sum(inArea))])
%

Nparams = length(MPR.k1)*length(MPR.k2)*length(MPR.k3);
disp(['Number of runs: ', num2str(Nparams)])
% preallocate
Nalive = zeros(Nparams,1);
n= 1;
for k1 = MPR.k1
    for k2 = MPR.k2
        for k3 = MPR.k3
            MP(n).k1 = k1;
            MP(n).k2 = k2;
            MP(n).k3 = k3;
            n = n+1;
        end
    end 
end


iniNumPerParticle = number_ind(1,1);
disp(['Initial biomass per particle: ',num2str(iniNumPerParticle)])
parfor n=1:Nparams  % parfor
    disp(['Run ',num2str(n)])
    N = ones(size(PP,1),1)*iniNumPerParticle; % inital number of individuals per particle
    % NOW RUN THE 'MODEL' 
    for t=1:length(nctime)
        N = N + N_change(N,PP(:,t),H(:,t),NS(:,t),dt,MP(n),HShallow);
    end
    Nend = N;
    Nalive(n) = sum(Nend);
    if Nalive(n) > Nini*0.5 && ...
            Nalive(n) < Nini*1.5
        disp('--------------------------------------')
        disp(['Nend: ',num2str(Nalive(n))])
        disp(MP(n))
%         disp(['Initial biomass per particle: ',num2str(iniNumPerParticle)])
%         figure
%         hist(Nend)
%         pause
    end
    
end


save('IBM5_par_4y_10k_2000-3.mat','Nalive','Nini','MP');

RunIDs = Nalive > Nini*0.8 & Nalive < Nini*1.2;



TopMP = MP(RunIDs);
for n=1:length(TopMP)
    disp(TopMP(n));
end

save('IBM5_Top_par_4y_10k_2000-3.mat','Nalive','Nini','TopMP','MP','HShallow');

%
% disp('Plotting ...')

figure
plot(Nalive/Nini,'xb')
ylim([0.5 1.5])


% EUlatlim = [30 46];
% EUlonlim = [-7 36.5];
%  LATS = lats(IDs,end);
%  LONS = lons(IDs,end);
%  hFig = figure(1);
%     set(hFig, 'Position', [10 10 1200 700]); % fixed figure dimensions
%     set(hFig,'Color','w')
%     clf
%     m_proj('mercator' ,'lon',[EUlonlim(1) EUlonlim(2)], ...
%         'lat',[EUlatlim(1) EUlatlim(2)]);    
%     m_grid('xtick',[EUlonlim(1):4:EUlonlim(2)],'ytick',[EUlatlim(1):2:EUlatlim(2)],'linestyle','none')
%     hold on
%     m_plot(LONS,LATS,'or','MarkerFaceColor','r')
%     shading flat
%     caxis([0.5 2.4])
%     colorbar
%     hold off
%     m_gshhs_l('patch',[.7 .7 .7],'edgecolor','black')
%     title('Average travel time [years]','interpreter','none');


% figure
% pcolor(pp)
% set(gca,'ColorScale','log')
% shading flat 
% colorbar


function dN = N_change(N,PP,H,NS,dt,MP,HShallow)
% change in number of ind. at each time step
   %   
   dN1 = N.*dt.* ( MP.k1 .* PP./ (MP.k2 + PP) - MP.k3); 
   % Penalize particles in shallow water with additional mortality
   Shallow = H*0;
   Shallow(H<HShallow) = -0.03; % daily mortality in shallow water - increased predation; -0.05 -> halved population in 14 days
   dN2 = N.*dt.* Shallow;
   % Penalize proximity to the coast - beaching
   Beaching = NS .* (-0.05);
   dN3 = N.*dt.* Beaching;
   dN = dN1 + dN2 + dN3;   
end
