## Carregando pacotes e funcoes necessarias ##
library(deSolve)
library(rootSolve)
source("eq_funcoes.r")

### Sistema Lotka-Volterra de competidores ###

## Coexistência das duas espécies, 1/beta < k1/k2 < alfa ##
plota.LV(n=c(n1=1,n2=1),r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,beta=0.1,time=150)

## Calculo dos valores no equilibrio
(lv.n1 <- lv.neq(r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,beta=0.1))

## Perturbando no equilibrio
## Grafico ##
plota.LV(n=lv.n1,r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,
         beta=0.1,time=150, perturb=c(-1,1), t.perturb=50)

## Coexistência das duas espécies, 1/beta > k1/k2 > alfa ##
## Calculo dos valores no equilibrio
(lv.n2 <- lv.neq(r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9))
## Grafico ##
plota.LV(n=lv.n2,r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9,time=150)

## Perturbando no equilibrio ##
plota.LV(n=lv.n2,r1=0.2,K1=150,r2=0.2,K2=100,
         alfa=1.8,beta=0.9,time=200, perturb=c(-1,1), t.perturb=50)

## Experimentos do Bruno e Liedson
plota.LV(n=lv.n2,r1=0.9,K1=150,r2=0.9,K2=100,
         alfa=1.8,beta=0.9,time=200, perturb=c(-40,0), t.perturb=50)

## Invertendo a perturbação ##
plota.LV(n=lv.n2,r1=0.2,K1=150,r2=0.2,K2=100,
         alfa=1.8,beta=0.9,time=200, perturb=c(1,-1), t.perturb=50)

## Calculo dos autovalores da Matriz da Comunidade ##

##Primeira combinação de parametros: 1/beta > k1/k2 > alfa  ##
## Matriz da comunidade
j1 <- jacob.lv(r1=0.2,K1=150,r2=0.2,K2=100,alfa=0.2,beta=0.1)
##Autovalores
eigen(j1, only.values=TRUE)$values

## Segunda combinação de parametros: 1/beta < k1/k2 < alfa  ## 
## Matriz da comunidade
j2 <- jacob.lv(r1=0.2,K1=150,r2=0.2,K2=100,alfa=1.8,beta=0.9)
##Autovalores
eigen(j2, only.values=TRUE)$values

### ESTUDO DE R. MAY ###
##Conectância e força de interação fixa, aumento do n de especies ##
(sim.1 <- may(S=20,C=0.3, f=0.2, nsim=100))
(sim.1 <- rbind(sim.1,may(S=40,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=60,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=80,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=90,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=100,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=110,C=0.3, f=0.2, nsim=100)))
(sim.1 <- rbind(sim.1,may(S=120,C=0.3, f=0.2, nsim=100)))
## Grafico
plot(p.estab~S, data=sim.1, xlab="N de espécies",ylab="Proporção matrizes estáveis")

##Riqueza e força de interação fixa, aumento de Conectancia ##
(sim.2 <- may(S=120,C=0.3, f=0.2))
(sim.2 <- rbind(sim.2,may(S=120,C=0.28, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.26, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.24, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.22, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.20, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.18, f=0.2)))
(sim.2 <- rbind(sim.2,may(S=120,C=0.16, f=0.2)))
##grafico
plot(p.estab~C, data=sim.2, xlab="Conectancia",ylab="Proporção matrizes estáveis")

### Matrizes de varios tamanhos e conectancias ###
## Sorteio de matrizes de 20 a 200 especies, 0.1 a 0.5 de conectncia
sim.3 <- yodzis(Smin=20,Smax=200,Cmin=0.1,Cmax=0.5, f=0.2, nsim=200)
plot(C~S,data=sim.3,subset=max.aut<0, ylab="Conectância", xlab="N de espécies")
## limite teorico
curve(1/(0.2^2*x), add=T, col="red", 20,200)