% Runs simulation for density dependent kicks
kappa= 1; %mortality proportion in affected sites
tau= 1; %interval between disturbances
L = 100; % number of sites in landscape
disturbed =[1:1:100]; % number of sites affected by disturbance
occupied = [1:1:100]; % number of sites occupied by population
T = 100; % number of flow-kick iterations
S = 100; % number of simulation runs
crashcountmatrix_dd = nan(length(occupied),length(disturbed));
disturbedcount = 1;

for D = disturbed
occupiedcount = 1;
for N = occupied
    curparam = [N,D];
    disp(curparam)
    crashcount=0;
    for j=1:S
    postkick=100; %initial population size
    count=0;
        for i=1:100
            count=count+1;
            k = hygernd(L,D,N); %draws a random number of sites to be affected from the hypergeometric distribution
            [T,X]=ode45(@dxdtAllee,[(count-1)*tau, count*tau],postkick); %flow 
            postkick=X(end)*(1-(kappa*(k/N)*(X(end)/100))); %kick
        end 
        if postkick < 1
            disp('crash')
            crashcount = crashcount+1;
        end 
    end
crashcountmatrix_dd(occupiedcount, disturbedcount) = crashcount;
occupiedcount = occupiedcount+1;
end
disturbedcount = disturbedcount+1;
end

save output.mat
save(strcat(['sddkcountoutput_high_res.mat']))