pp_plot_para = function(r1,K,n01,c,a,Th,d,n02,paratochange,inival,endval,stepval,animation) {
#This function plots simulated dynamics and isolines together (predator-prey)

#Prey:
# 1 r1 - intrinsic growth rate 
# 2 K - prey carrying capacity
# 3 n01 - initial population sizes 

#Predator:
# 4 c - prey conversion rate
# 5 a - encounter rate
# 6 Th - handling time
# 7 d - predator mortality
# 8 n02 - initial population size

#further parameters:
# paratochange: number of the para to be changed (see above)
# inival: initial parameter value
# endval: final parameter value
# stepval: parameter change step

#animation - whether to animate

para=NULL
para[1] = r1; 
para[2] = K; 
para[4] = c; 
para[5] = a;
para[6] = Th;
para[7] = d;
paraname = c("Prey growth rate","Prey carrying capacity","  ","Prey conversion rate","Encounter rate","Handling time","Predator mortality")

xx=NULL;
yy=NULL;

for (value in seq(inival,endval,by=stepval)) {
   para[paratochange] = value;
   aa=pp_lines(para[1],para[2],para[4],para[5],para[6],para[7])
   xx =c(xx, aa$y[,1], aa$R0);
   yy =c(yy, aa$y[,2]);

   if (aa$R0 == 0)
      break
   
}
maxx = max(xx);
maxy = max(yy);

   
for (value in seq(inival,endval,by=stepval)) {
   para[paratochange] = value;
   aa=pp_lines(para[1],para[2],para[4],para[5],para[6],para[7]);
   if (aa$R0 == 0)
      break
   text = paste(paraname[paratochange], " = ",as.character(value), ", R0 = ",format(aa$R0,digits=3), sep="")
   plot(aa$y[,1],aa$y[,2],type="l",lty=1,lwd=2,col="blue",xlim=c(0,maxx),ylim=c(0,maxy),
     xlab="Prey density",ylab="Predator density",main=text)
   lines(c(aa$R0, aa$R0),c(0, maxy),lty=1,lwd=2,col="green");
   if (animation > 0)
     Sys.sleep(1)
   
   
}
}

pp_lines = function(r1,K,c,a,Th,d) {
  nobin = 100

  if (c*a > d*a*Th)
    R0=sqrt(d/(c*a-d*a*Th))
  else {
    R0 = 0
  }

  incr1 = max(K,R0)/nobin;
  x = NULL;
  y1 = NULL;
  y2 = NULL;
  for (t in 1:nobin) {
    x1 = incr1*t ;
    if (x1>0)
      y1t = (r1/a)*((1-x1/K)*(1+a*Th*x1^2))/x1
    else
      y1t=0
   
   
    x = c(x, x1)
    y1=c(y1, y1t)
   
  }
  y=cbind(x, y1)

  list(y=y,R0=R0)
}