function  [y,R0] = pp_lines(r1,K,c,a,Th,d)
%This is a function calculating equilibrium lines for predator-prey dynamics 

%Prey:
%r1 - intrinsic growth rate 
%K - prey carrying capacity

%Predator:
%c - prey conversion rate
%a - encounter rate
%Th - handling time
%d - predator mortality

nobin = 100;

if (c*a > d*a*Th)
   R0=sqrt(d/(c*a-d*a*Th));
else
   'Nonrealistic values: Th or d too large'
   R0=0;
end

incr1 = max(K,R0)/nobin;
x = [];
y1 = [];
y2 = [];
for t=0:nobin
   x1 = incr1*t ;
   if (x1>0)
      y1t = (r1/a)*((1-x1/K)*(1+a*Th*x1^2))/x1;
   else
      y1t=0;
   end
   
   x = [x x1];
   y1=[y1 y1t];
   %y2=[y2 y2t];
end
%y=[x; y1; y2];
y=[x; y1];
   
   
   