######################################
#                                    #
#  Model dravec-kořist Gauseho typu  #
#                                    #
######################################

library("deSolve")

## Varianta 1: -delta>c*S
rk<-2
Kk<-1000
S<-4
sigma<-400
delta<--15   # delta = rp (musí být záporná)
c<-3
casy<-c(0:50)/10

-delta>c*S

# predační (trofická) funkce Holling II
p<-function(Nk,S,sigma) {
  return(S*Nk/(Nk+sigma))
}

plot(c(min(casy),max(casy)),c(0,1400),cex=0,main="Model dravec-kořist Gauseho typu",xlab="Čas",ylab="Velikost populace")

for(i in 1:3) {
  N0k<-500*i
  N0p<-500*(4-i)
  
  podminky<-c("Nk"=N0k,"Np"=N0p)
  soustava<-function(t,prom,param) {
    return(list(c(param["rk"]*prom["Nk"]*(1-prom["Nk"]/param["Kk"])-p(prom["Nk"],param["S"],param["sigma"])*prom["Np"],
                  prom["Np"]*param["delta"]+param["c"]*prom["Np"]*p(prom["Nk"],param["S"],param["sigma"]))))
  }
  parametry<-c("rk"=rk,"Kk"=Kk,"delta"=delta,"c"=c,"S"=S,"sigma"=sigma)
  metoda<-"lsoda"
  
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,"time"],reseni[,"Nk"],lwd=3,lty=1,col="black")
  lines(reseni[,"time"],reseni[,"Np"],lwd=3,lty=2,col="green")  

}

## Varianta 2: -delta>c*S a současně -delta>c*S*(Kk-sigma)/(Kk+sigma)
rk<-2
Kk<-1000
S<-2
sigma<-400
delta<--3.5   # delta = rp (musí být záporná)
c<-3
casy<-c(0:400)/10

delta<c*S
-delta>c*S*(Kk-sigma)/(Kk+sigma)

# predační (trofická) funkce Holling II
p<-function(Nk,S,sigma) {
  return(S*Nk/(Nk+sigma))
}

plot(c(min(casy),max(casy)),c(0,1400),cex=0,main="Model dravec-kořist Gauseho typu",xlab="Čas",ylab="Velikost populace")

for(i in 1:3) {
  N0k<-500*i
  N0p<-500*(4-i)
  
  podminky<-c("Nk"=N0k,"Np"=N0p)
  soustava<-function(t,prom,param) {
    return(list(c(param["rk"]*prom["Nk"]*(1-prom["Nk"]/param["Kk"])-p(prom["Nk"],param["S"],param["sigma"])*prom["Np"],
                  prom["Np"]*param["delta"]+param["c"]*prom["Np"]*p(prom["Nk"],param["S"],param["sigma"]))))
  }
  parametry<-c("rk"=rk,"Kk"=Kk,"delta"=delta,"c"=c,"S"=S,"sigma"=sigma)
  metoda<-"lsoda"
  
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,"time"],reseni[,"Nk"],lwd=3,lty=1,col="black")
  lines(reseni[,"time"],reseni[,"Np"],lwd=3,lty=2,col="green")  
  
}


## Varianta 3: -delta>c*S a současně -delta<c*S*(Kk-sigma)/(Kk+sigma)
rk<-2
Kk<-1000
S<-2
sigma<-400
delta<--2.5   # delta = rp (musí být záporná)
c<-3
casy<-c(0:400)/10

-delta<c*S
-delta<c*S*(Kk-sigma)/(Kk+sigma)

# predační (trofická) funkce Holling II
p<-function(Nk,S,sigma) {
  return(S*Nk/(Nk+sigma))
}

plot(c(min(casy),max(casy)),c(0,1400),cex=0,main="Model dravec-kořist Gauseho typu",xlab="Čas",ylab="Velikost populace")

for(i in 1:3) {
  N0k<-500*i
  N0p<-500*(4-i)
  
  podminky<-c("Nk"=N0k,"Np"=N0p)
  soustava<-function(t,prom,param) {
    return(list(c(param["rk"]*prom["Nk"]*(1-prom["Nk"]/param["Kk"])-p(prom["Nk"],param["S"],param["sigma"])*prom["Np"],
                  prom["Np"]*param["delta"]+param["c"]*prom["Np"]*p(prom["Nk"],param["S"],param["sigma"]))))
  }
  parametry<-c("rk"=rk,"Kk"=Kk,"delta"=delta,"c"=c,"S"=S,"sigma"=sigma)
  metoda<-"lsoda"
  
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,"time"],reseni[,"Nk"],lwd=3,lty=1,col="black")
  lines(reseni[,"time"],reseni[,"Np"],lwd=3,lty=2,col="green")  
  
}

# Směrové pole
sklon<-function(N1,N2) {
  d1<-rk*N1*(1-N1/Kk)-p(N1,S,sigma)*N2
  d2<-N2*delta+c*N2*p(N1,S,sigma)
  return(c(d1,d2))
}

uhel<-function(vektor) {
  return(Arg(vektor[1]*(1+0i)+vektor[2]*(0+1i)))
}

xlim<-c(0,1000)
ylim<-c(0,1200)

xn<-21
yn<-21

plot(c(0,0),cex=0,xlim=xlim,ylim=ylim,main="Směrové pole pro model dvou populací",
     xlab=expression(N[1]),ylab=expression(N[2]))

dx<-(xlim[2]-xlim[1])/xn
dy<-(ylim[2]-ylim[1])/yn

for(xi in 0:(xn-1)) {
  for(yi in 0:(yn-1)) {
    vektor<-sklon(xlim[1]+(xi+0.5)*dx,ylim[1]+(yi+0.5)*dy)
    if(sqrt(vektor[1]^2+vektor[2]^2)>10) {
      text(xlim[1]+(xi+0.5)*dx,ylim[1]+(yi+0.5)*dy,labels="→",srt=180/pi*uhel(vektor))
    } else {
      text(xlim[1]+(xi+0.5)*dx,ylim[1]+(yi+0.5)*dy,labels="●")
    }
  }
}


for(xs in c(xlim[1],(xlim[1]+xlim[2])/2,xlim[2])) {
  for(ys in c(ylim[1],(ylim[1]+ylim[2])/2,ylim[2])) {
    x<-xs
    y<-ys
    x2zpet<-x1zpet<-xstep<-Inf
    y2zpet<-y1zpet<-ystep<-Inf
    
    i<-0
    while(i<10000 & sqrt((x-x2zpet)^2+(y-y2zpet)^2)>sqrt(dx^2+dy^2)/100 & x>=xlim[1] & x<=xlim[2] & y>=ylim[1] & y<=ylim[2]) {
      vektor<-sklon(x,y)
      
      if(min(vektor==0)==TRUE) {
        points(x,y,pch=19,col="green")
      } else {
        xstep<-dx/4*vektor[1]/sqrt(vektor[1]^2+vektor[2]^2)
        ystep<-dy/4*vektor[2]/sqrt(vektor[1]^2+vektor[2]^2)
        lines(c(x,x+xstep),c(y,y+ystep),lwd=3,col="green")
      }
      
      x2zpet<-x1zpet; x1zpet<-x
      y2zpet<-y1zpet; y1zpet<-y
      
      x<-x+xstep
      y<-y+ystep
      
      i<-i+1
    }
  }
}




## Vyšetření řešení modelu (první varianta pomocí balíku nleqslv)
# Np*(delta+c*p(Nk,S,sigma)) = 0
# 1. Np = 0 (predátor nikdy nebyl přítomen)
# 2. delta = - c*p(Nk, S, sigma)

library("nleqslv")

# Ustálená velikost populace kořisti (z rovnice pro Np)
fce.k<-function(Nk) {
  delta+c*S*Nk/(Nk+sigma)
}
Nklim<-nleqslv(500,fce.k,control=list(ftol=0.0001))$x

# Ustálená velikost populace predátora (dosazením do rovnice pro Nk)
fce.p<-function(Np) {
  rk*Nklim*(1-Nklim/Kk)-Np*S*Nklim/(Nklim+sigma)
}
Nplim<-nleqslv(400,fce.p,control=list(ftol=0.0001))$x

## Vyšetření řešení modelu (druhá varianta pomocí balíku utils)
# Np*(delta+c*p(Nk,S,sigma)) = 0
# 1. Np = 0 (predátor nikdy nebyl přítomen)
# 2. delta = - c*p(Nk, S, sigma)

#install.packages("utils")
library("utils")

# Ustálená velikost populace kořisti (z rovnice pro Np)
fce.k<-function(Nk,delta,c,S,sigma) {
  delta+c*S*Nk/(Nk+sigma)
}
Nklim<-uniroot(fce.k,c(0,1000),tol=0.1,delta=delta,c=c,S=S,sigma=sigma)$root

# Ustálená velikost populace predátora (dosazením do rovnice pro Nk)
fce.p<-function(Np,Nklim,rk,Kk,S,sigma) {
  rk*Nklim*(1-Nklim/Kk)-Np*S*Nklim/(Nklim+sigma)
}
Nplim<-uniroot(fce.p,c(0,1000),tol=0.1,rk=rk,Kk=Kk,S=S,sigma=sigma,Nklim=Nklim)$root

lines(c(min(casy),max(casy)),rep(Nklim,2),col="red",lty=1)
lines(c(min(casy),max(casy)),rep(Nplim,2),col="red",lty=2)




## Varianty chování modelu
# 1) -delta > c * S -> populace dravce vymře
# 2) -delta < c * S -> populace koexistují
# 3) -delta > c * S * (Kk - sigma) / (Kk + sigma) -> populace se ustálí
# 3) -delta < c * S * (Kk - sigma) / (Kk + sigma) -> populace kolísají



#######################################
#                                     #
#  Model dravec-kořist Leslieho typu  #
#                                     #
#######################################

library("deSolve")

rk<-1
rp<-0.6
Kk<-1000
Kp<-1500
alfa<-0.0001
gama<-0.3
casy<-c(0:100)/10

plot(c(min(casy),max(casy)),c(0,2200),cex=0,main="Model dravec-kořist Leslieho typu",xlab="Čas",ylab="Velikost populace")

for(i in 1:3) {
  N0k<-700*i
  N0p<-700*(4-i)
  
  podminky<-c("Nk"=N0k,"Np"=N0p)
  soustava<-function(t,prom,param) {
    return(list(c(prom["Nk"]*(param["rk"]-param["rk"]*prom["Nk"]/param["Kk"]-param["alfa"]*prom["Np"]),
                  prom["Np"]*(param["rp"]-param["rp"]*prom["Np"]/(param["Kp"]+param["gama"]*prom["Nk"])))))
  }
  parametry<-c("rk"=rk,"Kk"=Kk,"rp"=rp,"Kp"=Kp,"alfa"=alfa,"gama"=gama)
  metoda<-"lsoda"
  
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  lines(reseni[,"time"],reseni[,"Nk"],lwd=3,lty=1,col="black")
  lines(reseni[,"time"],reseni[,"Np"],lwd=3,lty=2,col="green")  
}