pp_time = function (r1,K,n01,c,a,Th,d,n02,nostep) {
#This is a function simulating predator-prey dynamics with the type III functional response

#Prey:
#r1 - intrinsic growth rate 
#K - prey carrying capacity
#n01 - initial population size

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

#nostep - number of steps for simulation

if (c*a > d*a*Th)
   R0=sqrt(d/(c*a-d*a*Th))
else {
   print('Nonrealistic values: Th or d too large')
   print(d*Th)
   return()
}



  temp1 = rep(NA,nostep+1)
  temp2 = rep(NA,nostep+1)
  temp1[1] = n01
  temp2[1] = n02
  for (t in 1:nostep) {
     nprime1 = temp1[t] + r1*temp1[t]*(K-temp1[t])/K - temp2[t]*(a*temp1[t]^2/(1+a*Th*temp1[t]^2)) ;
     nprime2 = temp2[t] + c*temp2[t]*(a*temp1[t]^2/(1+a*Th*temp1[t]^2)) - d*temp2[t];
     if (nprime1 < 0 ) {
        nprime1 = 0;
     }
     if (nprime2 < 0) {
        nprime2 = 0;
     }
     temp1[t+1] = nprime1;
     temp2[t+1] = nprime2;
  }
  xx=cbind(temp1,temp2)
  list(xx=xx,R0=R0)
}