% The purpose of this code is to simulate theh Nicholson Bailey model for 
% host-parasitoid interaction

clear;
%set the initial data
N =4;
P = 2;

%set the number of iterates
M = 2000;

%set the parameters (the last one set is the one currently being used.)
r = 2;
K = 14.51;
a = 0.2;
c = 1;

r = 1.1;
a = 0.005;
c = 3;
K = 250;
K = 200;

%Leah's parameters:
r = 0.5;
K = 14.47;
c = 1;

r = 2.4;
K = 21.47
K = 14.47;
r = 2.65
K = 24.97;

q = 0.4;
apbar = r*(1-q);
Nbar = K*q;
Pbar = Nbar*(1-exp(-apbar));
a = apbar/Pbar;



% other parameters to try
% r = 1.1;
%a = 0.005;
% c = 3;
%K = 200, 250

%store the initial data
n(1) = N;
p(1) = P;


%iterate
for j=2:M

%evaluate the functions %The beddington model:
lim = exp(r*(1-N/K));
Prob = exp(-a*P);
N_new = N*Prob*lim;
P_new = c*N*(1-Prob);


%store the current values
n(j) = N_new;
p(j) = P_new;

%update
N = N_new;
P = P_new;
end

%make some plots
figure(1)
plot([1:M],n,'b',[1:M],p,'r',[1:M],n,'b*',[1:M],p,'r*')

figure(2)
plot(n,p,'*')