%Function to analytically estimate resilience to density dependent kicks, methods modified from Meyer
%et al 2018
function r = resilienceEstimateDD(T,D,H,f,m)
f = 1; %set frequency of disturbance
m = 1; %set proportion killed in affected sites
K = 100; %set carrying capacity
A = 20; % set allee threshold
P = [1:D]; %make a vector of occupied sites that could be affected by disturbance

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

if m*D/H >= 0.799
    r = 2;
    return; 
end 

%Create a range of kick values aligned with the above
krange= m/H:m/H:m*D/H;  %go close to distance to threshold, but not quite (CHECK THIS)

%Create an array to store points on the resilience boundary.
ResilienceCurve=zeros(1,length(krange));

%For each kick, calculate the minimum time needed to fully recover from the
%kick within the basin of attraction.
kcount=0;
for k=krange
    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) (change to allow for different K,A values)
    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(1-k(x0/100) to x0
        end
    ResilienceCurve(kcount)= min(TimeToFlow); %the tau that goes with k
end

%for each kick, calculate the tau estimated by the cumulative density
%function
tau_hcdf = f./(hygepdf(P,T,D,H)+hygecdf(P,T,D,H,'upper'));

%Now, see if the kappa tau estimates fall within the resilience threshold.
%If tau_hcdf-ResilienceCurve is ever negative, do not expect resilience,
%otherwise we do expect resilience, but are unsure

check = tau_hcdf- ResilienceCurve;

if min(check) < 0 
    r=1;
else
    r=0;
end 

end