##################################
#### Modelo SIR Epidemiologia ####
##################################
#####
##### CASO 1: S0 < a/b --> ponto inicial eh estavel
#####

# Removendo todos os objetos previamente guardados
rm(list=ls())

b=0.6
a=0.3

tempo.Max=100 # definindo tempo maximo a ser rodado
S=rep(0,tempo.Max) # vetor onde guardaremos os valores de Susceptiveis
I=rep(0,tempo.Max) # vetor onde guardaremos os valores de Infectados
R=rep(0,tempo.Max) # vetor onde guardaremos os valores de Recuperados

a/b 
(a-2)/b
S[1]=0.29 # valor S0 < a/b --> ponto estavel
I[1]=0.01 # pequena perturbacao no numero de infectados (eq preve que I=0 eh estavel)
R[1]=0.6

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  S[n+1]=S[n]-b*S[n]*I[n]
  I[n+1]=I[n]+b*S[n]*I[n]-a*I[n]
  R[n+1]=R[n]+a*I[n]
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,S,I,R)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "S", "I", "R")
tail(dados) # fim da tabela com os dados juntos (tempo, S(n), I(n), R(n))

# Plot da populacao com pequena perturbacao no numero de infectados
plot(dados$S~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,1), xlim=c(0,50), main="S(0) < a/b")
lines(dados$I~dados$tempo, col="red",type="o", pch=4)
lines(dados$R~dados$tempo, col="blue",type="o")
legend(35, 0.8, legend=c("Susceptiveis", "Infectados","Recuperados"),col=c("black", "red", "blue"),
       pch=16, cex=0.8)



#####
##### CASO 2: S0 > a/b --> ponto inicial nao eh estavel e vai para um ponto estavel em que S<a/b
#####

# Removendo todos os objetos previamente guardados
rm(list=ls())

b=0.6
a=0.3

tempo.Max=100 # definindo tempo maximo a ser rodado
S=rep(0,tempo.Max) # vetor onde guardaremos os valores de Susceptiveis
I=rep(0,tempo.Max) # vetor onde guardaremos os valores de Infectados
R=rep(0,tempo.Max) # vetor onde guardaremos os valores de Recuperados

a/b 
(a-2)/b
S[1]=0.79 # Agora note que S0 > a/b -> entao vai para outro ponto estavel em que S < a/b
I[1]=0.01
R[1]=0.2

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  S[n+1]=S[n]-b*S[n]*I[n]
  I[n+1]=I[n]+b*S[n]*I[n]-a*I[n]
  R[n+1]=R[n]+a*I[n]
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,S,I,R)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "S", "I", "R")
tail(dados) # fim da tabela com os dados juntos (tempo, S(n), I(n), R(n))

# Plot da populacao com pequena perturbacao nos infectados 
plot(dados$S~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,1), xlim=c(0,50), main="S(0) > a/b")
lines(dados$I~dados$tempo, col="red",type="o", pch=4)
lines(dados$R~dados$tempo, col="blue",type="o")
legend(35, 1, legend=c("Susceptiveis", "Infectados","Recuperados"),col=c("black", "red", "blue"),
       pch=16, cex=0.8)



#####
##### CASO 3: S0 > (a-2)/b e S0 < a/b --> o ponto inicial eh estavel
#####

# Removendo todos os objetos previamente guardados
rm(list=ls())

b=6
a=3

tempo.Max=100 # definindo tempo maximo a ser rodado
S=rep(0,tempo.Max) # vetor onde guardaremos os valores de Susceptiveis
I=rep(0,tempo.Max) # vetor onde guardaremos os valores de Infectados
R=rep(0,tempo.Max) # vetor onde guardaremos os valores de Recuperados

a/b
(a-2)/b 
S[1]=0.2 # Note que S0 > (a-2)/b e S0 < a/b --> o ponto inicial eh estavel
I[1]=0.01
R[1]=0.79

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  S[n+1]=S[n]-b*S[n]*I[n]
  I[n+1]=I[n]+b*S[n]*I[n]-a*I[n]
  R[n+1]=R[n]+a*I[n]
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,S,I,R)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "S", "I", "R")
tail(dados) # fim da tabela com os dados juntos (tempo, S(n), I(n), R(n))

# Plot da populacao com pequena perturbacao nos infectados
plot(dados$S~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,1), xlim=c(0,50), main="S(0) > (a-2)/b E S(0) < a/b")
lines(dados$I~dados$tempo, col="red",type="o", pch=4)
lines(dados$R~dados$tempo, col="blue",type="o")
legend(35, 1, legend=c("Susceptiveis", "Infectados","Recuperados"),col=c("black", "red", "blue"),
       pch=16, cex=0.8)



##################################
####  Parasitoide-Hospedeiro #####
#### Modelo Nicholson-Bailey #####
##################################

#####
##### Sem a perturbacao:
#####

# Removendo todos os objetos previamente guardados
rm(list=ls())

a=1
r=1
lambda = 2

# Ponto de equilibrio 
Hbar = lambda*log(lambda)/(lambda-1)
Pbar = log(lambda)
Hbar 
Pbar

tempo.Max=100 # definindo tempo maximo a ser rodado
H=rep(0,tempo.Max) # vetor onde guardaremos os valores
P=rep(0,tempo.Max) # vetor onde guardaremos os valores

H[1]=Hbar  # comeco no pto de eq 1
P[1]=Pbar

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  H[n+1]=lambda*H[n]*exp(-a*P[n])
  P[n+1]=r*H[n]*(1-exp(-a*P[n]))
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,H,P)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "H", "P")
tail(dados) # fim da tabela com os dados juntos (tempo, H(n), P[n])

# Plot da populacao com saturacao comecando no ponto de eq 1 sem perturbacao
plot(dados$H~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,5), xlim=c(0,30), main="Ponto eq 1 sem perturbacao")
lines(dados$P~dados$tempo, col="red",type="o", pch=16)
legend(10, 4, legend=c("Hospedeiro", "Parasitoide"),col=c("black", "red"), pch=16, cex=0.8)



#####
##### Fazendo a perturbacao:
#####

rm(list=ls())

a=1
r=1
lambda = 2

# Ponto de equilibrio 
Hbar = lambda*log(lambda)/(lambda-1)
Pbar = log(lambda)
Hbar 
Pbar

tempo.Max=100 # definindo tempo maximo a ser rodado
H=rep(0,tempo.Max) # vetor onde guardaremos os valores
P=rep(0,tempo.Max) # vetor onde guardaremos os valores

pert=1
H[1]=Hbar+pert  # comeco no pto de eq 1
P[1]=Pbar+pert

for(n in 1:(tempo.Max-1)){ #calculando iterativamente
  H[n+1]=lambda*H[n]*exp(-a*P[n])
  P[n+1]=r*H[n]*(1-exp(-a*P[n]))
}

tempo=seq(0,(tempo.Max-1))
dados=cbind(tempo,H,P)

dados=as.data.frame(dados)
colnames(dados)=c("tempo", "H", "P")
tail(dados) # fim da tabela com os dados juntos (tempo, H(n), P[n])

# Plot da populacao com saturacao comecando no ponto de eq 1 sem perturbacao
plot(dados$H~dados$tempo, type="o",pch=16, xlab="tempo", ylab="x(n)", 
     ylim=c(0,50), xlim=c(0,100), main="Ponto eq 1 com perturbacao")
lines(dados$P~dados$tempo, col="red",type="o", pch=16)
legend(10, 50, legend=c("Hospedeiro", "Parasitoide"),col=c("black", "red"), pch=16, cex=0.8)







##############################################
####    2 especies -  Modelo Nao Linear  #####
####  Parasitoide-Hospedeiro (May 1978)  #####
####       Periodo da oscilacao          #####
##############################################

####### Caso que a extincao nao eh estavel e aumentando o valor de f

rm(list=ls())

T.max=200

P=matrix(0,T.max,1)
H=matrix(0,T.max,1)

# Primeiro valor dos vetores das spp
P0=3
H0=5

a=1

#Calculando para cada  instante de tempo
Par.Hosp<-function(par.a,par.f,H.0,P.0,t.max){
  P=matrix(0,T.max,1)
  H=matrix(0,T.max,1)
  P[1]=P.0
  H[1]=H.0
  a=par.a
  f=par.f
  for (n in 1:t.max){
    H[n+1]= f*H[n]/(1+a*P[n])
    P[n+1]= H[n]*(1-1/(1+a*P[n]))
  }
  tempo=seq(0,t.max)
  dados=data.frame(tempo,H,P)
  return(dados)
}

#Juntando os dados
dados.f1.1=Par.Hosp(a,1.1,P0,H0,T.max)

dados.f5=Par.Hosp(a,5,P0,H0,T.max)

dados.f10=Par.Hosp(a,10,P0,H0,T.max)


#Valores baixos de f
plot(dados.f1.1$tempo, dados.f1.1$H, type="b", pch=16,xlim=c(0,100), ylim=c(0,10),col="#DCDCDC", 
     xlab="tempo", ylab="Numero de Individuos")
lines(dados.f1.1$tempo, dados.f1.1$P, type="b", pch=16)
legend(60, 10, legend=c("Hospedeiro", "Parasitoide"),col=c("#DCDCDC", "#000000"), pch=16, cex=0.8)


# Variando o valor de f para hospedeiros
plot(dados.f1.1$tempo, dados.f1.1$H, type="b", pch=16,xlim=c(0,20), ylim=c(0,35),col="#DCDCDC", 
     xlab="tempo", ylab="Numero de Hospedeiros", main="Hospeideiros") 
lines(dados.f5$tempo, dados.f5$H, pch=16, type="b", lty=2, col="#A9A9A9")
lines(dados.f10$tempo, dados.f10$H, pch=16, type="b", lty=2, col="#778899")
legend(18, 35, legend=c("f=2", "f=5", "f=10"),
       col=c("#DCDCDC", "#A9A9A9", "#778899"), pch=16, cex=0.8)



# Variando o valor de f para parasitoides:
plot(dados.f1.1$tempo, dados.f1.1$P, type="b", pch=16,xlim=c(0,20), ylim=c(0,35),col="#DCDCDC", 
     xlab="tempo", ylab="Numero de Parasitoides", main="Parasitoides") 
lines(dados.f5$tempo, dados.f5$P, pch=16, type="b", lty=2, col="#A9A9A9")
lines(dados.f10$tempo, dados.f10$P, pch=16, type="b", lty=2, col="#778899")
legend(18, 35, legend=c("f=2", "f=5", "f=10"),
       col=c("#DCDCDC", "#A9A9A9", "#778899"), pch=16, cex=0.8)