#########################################
###### Exemplo Juvenis e Adultos ########
#########################################
rm(list=ls())

#Tamanho inicial da Populacao
J0<-0
A0<-1

#Parametros
b<-10
p<-0.1
s<-0.2

matriz.JA=matrix(c(0, b, p, s), nrow=2, ncol=2, byrow=TRUE) #matriz com os paramentros
matriz.JA

#Calculando os autovalores e respectivos autovetores da matriz
eigen(matriz.JA)
eigenva.JA=eigen(matriz.JA)$values[1] # pegando o maior autovalor em modulo
eigenvc.JA=eigen(matriz.JA)$vectors[,1]  # pegando o espectivo autovetor do maior autovalor em modulo


########################################################### 
#### Calculando o numero de individuos iterativamente: ####
###########################################################

#Tempo total a ser analisado
T.max<-1000

#Definicao dos vetores de cada categoria e a populacao total para cada instante de tempo
J=rep(0,(T.max+1))
A=rep(0,(T.max+1))
Ntotal=J+A

#Primeiro valor dos vetores de categorias
J[1]<-J0
A[1]<-A0
Ntotal[1]<-J[1]+A[1]

#Calculando para cada instante de tempo (iterativamente):
for (n in 1:T.max){
  J[n+1]<-A[n]*b
  A[n+1]<-J[n]*p+A[n]*s  
  Ntotal[n+1]<-A[n+1]+J[n+1]
}

#Juntando os dados
tempo=seq(0,(T.max))
dados=data.frame(tempo,J,A, Ntotal)
head(dados)
tail(dados)


#Definido a Eq geral da solucao, completa
EqGeralJA <- function(n){
 vn=0.549719*((1.104988)^n)*c(9.049876,1) + 0.4502481*((-0.9049876)^n)*c(-11.04988,1) 
 vn=as.vector(vn) 
 return(vn)
}


#Definido a Eq geral da solucao, mas levando em conta so o autovalor dominante e seu respectivo autovetor
EqGeralJA.dom <- function(n){
  vn=0.549719*((1.104988)^n)*c(9.049876,1) 
  vn=as.vector(vn) 
  return(vn)
}

# Comparando os valores no tempo 1000 para cada solucao feita:
EqGeralJA(1000) 
EqGeralJA.dom(1000)
dados[1001,]

#definido uma matriz de duas colunas para calcular valores pela eq geral (completa e so com dominante)
Vn=matrix(nrow=1000,ncol=2)
Vn.dom=matrix(nrow=1000,ncol=2)

#Preenchendo valores pela Eq Geral
for (n in 1:T.max){
  Vn[n,]=EqGeralJA(n)
}

#Preenchendo valores pela Eq Geral so com dominante
for (n in 1:T.max){
  Vn.dom[n,]=EqGeralJA.dom(n)
}

#Montando as tabelas:
Tempo=seq(1,(T.max))
Res.EGJA=cbind(Tempo,Vn)
Res.EGJA=as.data.frame(Res.EGJA)
colnames(Res.EGJA)=c("tempo", "J", "A")
tail(Res.EGJA)

Res.EGJA.dom=cbind(Tempo,Vn.dom)
Res.EGJA.dom=as.data.frame(Res.EGJA.dom)
colnames(Res.EGJA.dom)=c("tempo", "J", "A")
tail(Res.EGJA.dom)



# diamante:iterativo, bola vazada: solucao geral, x: solucao geral so com dominante
# vermelho: Juvenis e azul: Adultos
plot(dados$tempo, dados$J, type="b", pch=18,ylim=c(0,10000),xlim=c(0,100),
     col="red", xlab="tempo", ylab="Numero de Individuos", cex=.9) 
lines(dados$tempo, dados$A, pch=18, col="blue", type="b", lty=2,cex=.9)
lines(Res.EGJA$tempo, Res.EGJA$J, pch=1, col="red", type="b", lty=2)
lines(Res.EGJA$tempo, Res.EGJA$A, pch=1, col="blue", type="b", lty=2)
lines(Res.EGJA.dom$tempo, Res.EGJA.dom$J, pch=4, col="red", type="b", lty=2)
lines(Res.EGJA.dom$tempo, Res.EGJA.dom$A, pch=4, col="blue", type="b", lty=2)
legend(1, 10000, legend=c("Juvenis", "Adultos"),col=c("red", "blue"), lty=1:2, cex=0.8)


#Zoom inicio:
plot(dados$tempo, dados$J, type="b", pch=18,ylim=c(0,15),xlim=c(0,10),
     col="red", xlab="tempo", ylab="Numero de Individuos",cex=.9) 
lines(dados$tempo, dados$A, pch=18, col="blue", type="b", lty=2,cex=.9)
lines(Res.EGJA$tempo, Res.EGJA$J, pch=1, col="red", type="b", lty=2)
lines(Res.EGJA$tempo, Res.EGJA$A, pch=1, col="blue", type="b", lty=2)
lines(Res.EGJA.dom$tempo, Res.EGJA.dom$J, pch=4, col="red", type="b", lty=2)
lines(Res.EGJA.dom$tempo, Res.EGJA.dom$A, pch=4, col="blue", type="b", lty=2)


#Zoom meio Juvenis:
plot(dados$tempo, dados$J, type="b", pch=18,ylim=c(100,150),xlim=c(30,35),
     col="red", xlab="tempo", ylab="Numero de Individuos",cex=.9) 
lines(Res.EGJA$tempo, Res.EGJA$J, pch=1, col="red", type="b", lty=2)
lines(Res.EGJA.dom$tempo, Res.EGJA.dom$J, pch=4, col="red", type="b", lty=2)

#Zoom meio Adultos:
plot(dados$tempo, dados$A, type="b", pch=18,ylim=c(10,20),xlim=c(30,35),
     col="blue", xlab="tempo", ylab="Numero de Individuos",lty=2,cex=.9)
lines(Res.EGJA$tempo, Res.EGJA$A, pch=1, col="blue", type="b", lty=2)
lines(Res.EGJA.dom$tempo, Res.EGJA.dom$A, pch=4, col="blue", type="b", lty=2)