######################################################
#                                                    #
#  Populace pod tlakem nespecializovaného predátora  #
#                                                    #
######################################################

library("deSolve")

# Definice predační funkce p (Nk = kritické N, nad kterým je hodnota p rovna S = saturace)

p<-function(N,Nk,S) {
  if (N <= Nk) {
    return((S/Nk)*N)
  }
  else{
    return(S)
  }
}

q<-function(N,Nk,S) {
  return((S*N)/(N+Nk))
}

r<-function(N,Nk,S) {
  return(2*S/(1+exp(-N/1000))-S)
}

plot(c(0:10000),r(c(0:10000),Nk,S),type="l")


# Diskrétní případ
N0<-500
r<-1
K<-1000
Nkrit<-200
S<-300

podminky<-c("N"=N0)
casy<-0:10
soustava<-function(t,prom,param) {
  return(list(c(prom["N"]+param["r"]*prom["N"]*(1-prom["N"]/param["K"])
                -p(prom["N"],param["Nkrit"],param["S"]))))
}
parametry<-c("r"=r,"K"=K,"Nkrit"=Nkrit,"S"=S)
metoda<-"iteration"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="p",col="brown",lwd=3,which="N")


# Spojitý případ
N0<-10
r<-1
K<-1000
Nkrit<-200
S<-100

podminky<-c("N"=N0)
casy<-0:100
soustava<-function(t,prom,param) {
  return(list(c(param["r"]*prom["N"]*(1-prom["N"]/param["K"])
                -p(prom["N"],param["Nkrit"],param["S"]))))
}
parametry<-c("r"=r,"K"=K,"Nkrit"=Nkrit,"S"=S)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="brown",lwd=3,which="N")
