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

library("deSolve")

# Model konkurekce (ovlivnění úživnosti)
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))
}

# Graf funkce úživnosti pro hodnoty N od 0 do 1000
plot(0:1000,kappa(0:1000,K=1000,C=1600,gamma=0.01),type="l",lwd=3,col="navyblue",xlab="N",ylab="kappa(N)")

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

reseni<-ode(podminky,casy,soustava,parametry,metoda)
plot(reseni,type="l",col="black",lwd=3,which="N1",ylim=c(0,1100),main="Velikosti populací")
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",main="Velikosti populací")
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 (ovlivnění úživnosti)

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))
}

# Graf funkce úživnosti pro hodnoty N od 0 do 1000
plot(0:1000,kappa(0:1000,K=1000,C=1600,gamma=0.01),type="l",lwd=3,col="navyblue",xlab="N",ylab="kappa(N)")


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

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

# vnořené cykly
plot(c(0,15),c(300,2400),cex=0,xlab="t",ylab="N",main="Velikosti populací")
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)
  }
}


# Model kompetice (ovlivnění koeficientu růstu)
plot(0,0,cex=0,xlim=c(0,20),ylim=c(0,2200),main="Velikosti populací",xlab="t",ylab="N(t)")

for(i in c(1:3)) {
  N10<-700*i
  N20<-700*(4-i)
  r1<-1
  r2<-0.6
  K1<-1000
  K2<-1500
  beta12<--0.0005
  beta21<--0.0002
  
  
  podminky<-c("N1"=N10,"N2"=N20)
  casy<-0:20
  soustava<-function(t,prom,param) {
    return(list(c(prom["N1"]*(param["r1"]-param["r1"]/param["K1"]*prom["N1"]+param["beta12"]*prom["N2"]),
                  prom["N2"]*(param["r2"]-param["r2"]/param["K2"]*prom["N2"]+param["beta21"]*prom["N1"]))))
  }
  parametry<-c("r1"=r1,"r2"=r2,"K1"=K1,"K2"=K2,"beta12"=beta12,"beta21"=beta21)
  metoda<-"lsoda"
  
  reseni<-ode(podminky,casy,soustava,parametry,metoda)
  
  lines(reseni[,"time"],reseni[,"N1"],lwd=3,col="black",lty=1)
  lines(reseni[,"time"],reseni[,"N2"],lwd=3,col="green",lty=2)
}


# Směrové pole
#############################
# Vykreslení směrového pole #
#############################

r1<-1.0
r2<-0.6
K1<-1000
K2<-1500
γ12<-0.5
γ21<-0.3
C12<-600
C21<-1000


sklon<-function(N1,N2) {
  d1<-(r1*N1*(1-(N1*(1+γ12*N2))/(K1+C12*γ12*N2)))
  d2<-(r2*N2*(1-(N2*(1+γ21*N1))/(K2+C21*γ21*N1)))
  return(c(d1,d2))
}

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

xlim<-c(300,900)
ylim<-c(600,1400)

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
    
    while(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
    }
  }
}





