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

# 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
Nk<-200
S<-300

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

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


# Spojitý případ

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

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

soustava2<-function(t,prom,param) {
  return(list(c(param[1]*prom[1]*(1-prom[1]/param[2])-q(prom[1],param[3],param[4]))))
}

reseni2<-ode(podminky,casy,soustava2,parametry,metoda)
lines(reseni2[,1],reseni2[,2],col="cornflowerblue",lwd=3)

soustava3<-function(t,prom,param) {
  return(list(c(param[1]*prom[1]*(1-prom[1]/param[2])-r(prom[1],param[3],param[4]))))
}

reseni3<-ode(podminky,casy,soustava3,parametry,metoda)
lines(reseni3[,1],reseni3[,2],col="gold",lwd=3)



###########################
#                         #
#  Interagující populace  #
#                         #
###########################


# Model konkurekce

r1<-1
r2<-0.6
K1<-1000
K2<-1500
C1<-600
C2<-1000
gamma1<-0.5
gamma2<-0.3
N10<-300
N20<-300

kappa<-function(N,K,C,gamma){
  return((K+C*gamma*N)/(1+gamma*N))
}

podminky<-c(N1=N10,N2=N20)
casy<-0:15
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1])*(1-(prom[1]/kappa(prom[2],param[3],param[5],param[7]))),
                prom[2]*(param[2])*(1-(prom[2]/kappa(prom[1],param[4],param[6],param[8]))))))
}
parametry<-c(r1,r2,K1,K2,C1,C2,gamma1,gamma2)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1100))
lines(reseni[,1],reseni[,3],col="green",lty=2,lwd=3)

# vnořené cykly
plot(c(0,15),c(300,1200),cex=0,xlab="t",ylab="N")
for(i in 1:4){
  for(j in 1:4){
    N10<-300*i
    N20<-300*j
    podminky<-c(N1=N10, N2=N20)
    reseni<-ode(podminky,casy,soustava,parametry,metoda)
    lines(reseni[,1],reseni[,2],lwd=3,col="black")
    lines(reseni[,1],reseni[,3],lwd=3,col="green",lty=2)
  }
}


# Model symbiózy

r1<-1
r2<-0.6
K1<-1000
K2<-1500
C1<-1600
C2<-2000
gamma1<-0.5
gamma2<-0.3
N10<-300
N20<-300

kappa<-function(N,K,C,gamma){
  return((K+C*gamma*N)/(1+gamma*N))
}

podminky<-c(N1=N10, N2=N20)
casy<-0:15
soustava<-function(t,prom,param) {
  return(list(c(prom[1]*(param[1])*(1-(prom[1]/kappa(prom[2],param[3],param[5],param[7]))),
                prom[2]*(param[2])*(1-(prom[2]/kappa(prom[1],param[4],param[6],param[8]))))))
}
parametry<-c(r1,r2,K1,K2,C1,C2,gamma1,gamma2)
metoda<-"lsoda"

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1100))
lines(reseni[,1],reseni[,3],col="green",lty=2,lwd=3)

# vnořené cykly
plot(c(0,15),c(300,2100),cex=0,xlab="t",ylab="N")
for(i in 1:4){
  for(j in 1:4){
    N10<-500*i
    N20<-600*j
    podminky<-c(N1=N10, N2=N20)
    reseni<-ode(podminky,casy,soustava,parametry,metoda)
    lines(reseni[,1],reseni[,2],lwd=3,col="black")
    lines(reseni[,1],reseni[,3],lwd=3,col="green",lty=2)
  }
}