%This script produces figure 1

%Panel a
%CALCULATE FLOW-KICK TRAJECTORIES FOR A POPULATION AFFECTED BY
%A PROPORTIONAL KICK
%The population's undisturbed dynamics are given by the ODE function dxdtAllee.m

figure
t=tiledlayout(2,2);
nexttile
ax = gca;
ax.TitleHorizontalAlignment = 'left';
title('a) Proportional Kicks')
hold on

%Trajectory A
kappa=-0.12; %kick proportion and direction
tau=0.25; %flow time
N=64; % number of flow-kick iterations
postkick=100; %initial population size
count=0;


for i=1:N
    count=count+1;
    [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow 
    postkick=X(end)* (1+kappa); %kick
    p1= plot(T,X,'Color',[255/256 128/256 0],'LineWidth',1);
end

%Trajectory B
kappa=-0.6; %kick size and direction
tau=10/8; %flow time
N=15;
postkick=100; %initial position
count=0;

for i=1:N
    count=count+1;
    [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow
    postkick=X(end)*(1+kappa); %kick
    p2 = plot(T,X,'Color',[0.4940 0.1840 0.5560],'LineWidth',1);
end

%Trajectory C
kappa=-0.8; %kick size and direction
tau=10/6; %flow time
N=15;
postkick=100; %initial position
count=0;

for i=1:N
    count=count+1;
    [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow
    postkick=X(end)*(1+kappa); %kick
    p3 = plot(T,X,'Color',[0.3010 0.7450 0.9330],'LineWidth',1);
end

axis([0 12 0 105])
xlabel('Time (years)')
ylabel('Population size')

% Panel b
% CALCULATES FLOW-KICK TRAJECTORIES FOR A POPULATION AFFECTED BY
%A DENSITY DEPENDENT KICK
%The population's undisturbed dynamics are given by the ODE function dxdtAllee.m

%Trajectory A
kappa=-0.12; %kick proportion and direction
tau=0.25; %flow time
N=64; % number of flow-kick iterations
postkick=100; %initial population size
count=0;

nexttile(3)
ax = gca;
ax.TitleHorizontalAlignment = 'left';
title('b) Density-Dependent Kicks')
hold on

for i=1:N
    count=count+1;
    [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow 
    postkick=X(end)* (1+(kappa*(X(end)/100))); %kick
    plot(T,X,'Color',[255/256 128/256 0],'LineWidth',1)
end

%Trajectory B
kappa=-0.6; %kick size and direction
tau=10/8; %flow time
N=15;
postkick=100; %initial position
count=0;

for i=1:N
    count=count+1;
    [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow
    postkick=X(end)* (1+(kappa*(X(end)/100))); %kick
    plot(T,X,'Color',[0.4940 0.1840 0.5560],'LineWidth',1)
end

%Trajectory C
kappa=-0.8; %kick size and direction
tau=10/6; %flow time
N=15;
postkick=100; %initial position
count=0;

for i=1:N
    count=count+1;
    [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow
    postkick=X(end)* (1+(kappa*(X(end)/100))); %kick
    plot(T,X,'Color',[0.3010 0.7450 0.9330],'LineWidth',1)
end

axis([0 12 0 105])
xlabel('Time (years)')
ylabel('Population size')

%
%CALCULATE THE FLOW-KICK RESILIENCE BOUNDARY FOR 2 KICK FUNCTIONAL
%FORMS - Modified from Meyer et al 2018
%For each kick size, the script approximates the minimum time it takes to
%population to recover by that magnitude. 

%Create the inverse of the vector field function, i.e. dt/dx
fun=@(x) 1./(x.*(1-(x/100)).*((x/20)-1)); 

%Create a range of kick values 
%for proportional or density-dependent kicks
krange_p_d= 0.01:0.01:0.79;  %go close to distance to threshold, but not quite

%Create an array to store points on the resilience boundary.
%proportional
ResilienceCurve_p=zeros(2,length(krange_p_d));
%density-dependent
ResilienceCurve_d=zeros(2,length(krange_p_d));

%For each kick, calculate the minimum time needed to fully recover from the
%kick within the basin of attraction.

%proportional
kcount=0;
for k=krange_p_d
    kcount=kcount+1;
    x0range=ceil(10*(20.1/(1-k)))/10:0.1:99.9; %identify a set of x values at which to start flowing (min computed over this set)
    TimeToFlow=zeros(1,length(x0range));
    x0count=0;
        for x0=x0range
            x0count=x0count+1;
            TimeToFlow(x0count)=integral(fun,x0*(1-k),x0);%calculate time to flow from x0+k to x0
        end
    ResilienceCurve_p(1,kcount)= min(TimeToFlow); %the tau that goes with k
    ResilienceCurve_p(2,kcount)= k; % the k
end

%density-dependent
kcount=0;
for k=krange_p_d
    kcount=kcount+1;
    x0range=ceil(10*(20.1/(1-k)))/10:0.1:99.9; %identify a set of x values at which to start flowing (min computed over this set) 
    TimeToFlow=zeros(1,length(x0range));
    x0count=0;
        for x0=x0range
            x0count=x0count+1;
            TimeToFlow(x0count)=integral(fun,x0*(1-k*(x0/100)),x0);%calculate time to flow from x0+k to x0
        end
    ResilienceCurve_d(1,kcount)= min(TimeToFlow); %the tau that goes with k
    ResilienceCurve_d(2,kcount)= k; % the k
end
%% Panel c
%Plots the resilience boundaries for proportional and density dependent kicks

nexttile(2)
ax = gca;
ax.TitleHorizontalAlignment = 'left';
title('c) Resilience Boundaries')
hold on
p4 = patch([ResilienceCurve_d(2,:) .79*ones(size(ResilienceCurve_d(2,:))) 0.79 1 1],[ResilienceCurve_d(1,:) fliplr(ResilienceCurve_d(1,:)) 6 6 0], 'r','FaceAlpha',.3, 'EdgeColor','none');
p5 = patch([ResilienceCurve_d(2,:) fliplr(ResilienceCurve_p(2,:))],[ResilienceCurve_d(1,:) fliplr(ResilienceCurve_p(1,:))], 'y','FaceAlpha',.3, 'EdgeColor','none');
p6 = patch([ResilienceCurve_p(2,:) .79 0 min(ylim)*ones(size(ResilienceCurve_p(2,:)))],[ResilienceCurve_p(1,:) 6 6 fliplr(ResilienceCurve_p(1,:))], 'b','FaceAlpha',.3, 'EdgeColor','none');
ylabel('flow time (\tau)')
p7 = plot(.12, .25, '.','Color',[255/256 128/256 0], 'MarkerSize', 18);
p8 = plot(.6, 1.25, '.','Color',[0.4940 0.1840 0.5560], 'MarkerSize', 18);
p9 = plot(.8, 1.66, '.','Color',[0.3010 0.7450 0.9330], 'MarkerSize', 18);
p10 = plot(.12,.25, 'LineStyle', 'none');
p11 = plot(.6, 1.25,'LineStyle', 'none');
hold off
xlabel('kick size (\kappa)')
axis([0 1 0 6])
lgd = legend([p1,p2,p3,p4,p5,p6,p7,p8,p9,p10,p11],{'        \kappa = 0.12, \tau = 0.25','        \kappa = 0.6, \tau = 1.25','        \kappa = 0.8, \tau = 1.66','Not Resilient (both kick types)' ,'Not Resilient (P) and Resilient (DD)','Resilient (both kick types)','','','','',''}, 'FontSize',12);
lgd.Layout.Tile = 4;
lgd.NumColumns = 2;