###############################################################
#                                                             #
# Ukázka vykreslení směrového pole (bude v nových skriptech)  #
#                                                             #
###############################################################

##########################################################
# Varianta 1: model dvou populací s ovlivněním úživnosti #
##########################################################
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))
}
xlim<-c(300,900)
ylim<-c(600,1400)


################################################
# Varianta 2: model dravec-kořist Gauseho typu #
################################################
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

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

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

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


#########################################
# Pokračování společné pro obě varianty #
#########################################

## Určení úhlu šipky
uhel<-function(vektor) {
  return(Arg(vektor[1]*(1+0i)+vektor[2]*(0+1i)))
}

## Definice počtu šipek vodorovně a svisle
xn<-21
yn<-21

## Prázdný graf
plot(c(0,0),cex=0,xlim=xlim,ylim=ylim,xlab="N1",ylab="N2")

## Velikost obdélníků pro šipky
dx<-(xlim[2]-xlim[1])/xn
dy<-(ylim[2]-ylim[1])/yn

## Vykreslení šipek
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="●")
    }
  }
}

## Vykreslení čar z devíti bodů v grafu (rohy + středy stran)
limitsteps<-1000 # maximální počet iterací
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
    stepsdone<-0
    
    while(sqrt((x-x2zpet)^2+(y-y2zpet)^2)>sqrt(dx^2+dy^2)/100 & abs(x)<Inf & abs(y)<Inf & stepsdone<limitsteps) {
      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
      
      stepsdone<-stepsdone+1
    }
    if(stepsdone==limitsteps) {
      warning(paste0("Maximum number of steps (",limitsteps,") was reached when starting at [",xs,"; ",ys,"]"))
    }
  }
}
box() # rámeček
