#########################################################
######## Modelos em tempo continuo nao-lineares ########
########################################################


######################################
######## Competicao (de novo) ########
######################################

#### Definindo a funcao de competicao:
Competicao <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dN1 = r1*N1*((K1 - N1 - beta1*N2)/K1)
    dN2 = r2*N2*((K2 - N2 - beta2*N1)/K2)
    return(list(c(dN1, dN2)))
  })
}

###############################################################
##### Caso 2 sa aula passada (sp 2 sempre exclui a sp 1): #####
##############################################################

#Parametros do modelo
Pars <- c(r1 = 0.05, r2=0.03, beta1 = 2, beta2 = 0.25, K1=40, K2=50)
#Tamanho inicial das populacoes
State <- c(N1 = 10, N2 =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 300, by = 1)

#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = Competicao, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out

par(mfrow=c(1,2)) #Plotando 2 resultados na mesma tela

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Plotando a densidade de uma sp pela outra (escala de cinza de acordo com o tempo, 
# onde mais claro eh o inicio)
cinzas=  gray.colors(dim(out)[1], start = 0.8, end = 0, gamma = 0.7, alpha = NULL)

plot(1, type="n", xlab="", ylab="", xlim=c(0, 100), ylim=c(0, 210))
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2)
points(out$N2~out$N1, col=cinzas,xlim=c(0,100), ylim= c(0,210), 
       ylab="N2", xlab="N1", pch=16)

dev.off()


#####################################################
######## Competicao com predacao intraguilda ########
#####################################################

#### Definindo a funcao de competicao com predacao intraguilda:
CompeticaoPredIntra <- function (Time, State, Pars) {
  with(as.list(c(State, Pars)), {
    dN1 = r1*N1*((K1 - N1 - beta1*N2)/K1) + gama*N1*N2
    dN2 = r2*N2*((K2 - N2 - beta2*N1)/K2) - delta*N1*N2
    return(list(c(dN1, dN2)))
  })
}

########################################################################
#####  Olhando pro caso que no modelo de competica original   ##########
##########         a gente tinha persistencia so da sp 2      ##########
########################################################################

Pars <- c(r1 = 0.05, r2=0.03, beta1 = 2, beta2 = 0.25, K1=40, K2=50,gama=0.002, delta=0.0002)
#Tamanho inicial das populacoes
State <- c(N1 = 10, N2 =10)
#Tempo a ser rodado a solucao numerica
Time <- seq(0, 500, by = 1)

#Saidas gravadas no arquivo out
out <- as.data.frame(ode(func = CompeticaoPredIntra, y = State, parms = Pars, times = Time))
colnames(out)=c("tempo", "N1", "N2")
out


par(mfrow=c(1,2)) #Plotando 2 resultados na mesma tela

# Plotando os tamanhos de cada especie pelo tempo
plot(out$N1~out$tempo , type = "l", xlab = "tempo", ylab = "Numero de individuos", 
     col=1, lty=1, ylim=c(0,50))
lines(out$N2~out$tempo , type = "l", col=2, lty=2)
legend("topright", c("N1", "N2"), lty = c(1,2), col = c(1,2), box.lwd = 0)

# Plotando a densidade de uma sp pela outra (escala de cinza de acordo com o tempo, 
# onde mais claro eh o inicio)
cinzas=  gray.colors(dim(out)[1], start = 0.8, end = 0, gamma = 0.7, alpha = NULL)

plot(1, type="n", xlab="", ylab="", xlim=c(0, 100), ylim=c(0, 210))
legend("topright", c("dN1/dt=0", "dN2/dt=0"), lty = c(1,2), box.lwd = 0)
abline(a=Pars[5]/Pars[3],b=-(Pars[5]/Pars[3])/Pars[5],lty=1, col="#C0C0C0")
abline(a=Pars[5]/(Pars[3]-Pars[7]*Pars[5]/Pars[1]),b=-Pars[5]/(Pars[3]-Pars[7]*Pars[5]/Pars[1])/Pars[5],lty=1)
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/Pars[4]),lty=2, col="#C0C0C0")
abline(a=Pars[6],b=-Pars[6]/(Pars[6]/(Pars[4]+Pars[8]*Pars[6]/Pars[2])),lty=2)
points(out$N2~out$N1, col=cinzas,xlim=c(0,100), ylim= c(0,210), 
       ylab="N2", xlab="N1", pch=16)

dev.off()