#====================================================================================
#Zorros y conejos
#===================================================================================
#install.packages('deSolve')
#install.packages('ggplot2')
library("deSolve")
library("ggplot2")

#Considera a los parametros como variables exogenas que no pueden ser modificadas
parameters<-c(vida.media.zorros= 4,#[año]
              vida.media.conejos= 2.5, #[año]
              tasa.natalidad.conejos= 2, #0.25, #[numero de conejos/conejo]
              poblacion.sostenible.conejos= 500 #[numero de conejos]
              )

#Para definir tus condiciones iniciales utiiza tus variables de estado
InitialConditions <- c(poblacion.conejos= 500, #[número de conejos]= poblacion.sostenible.conejos
                       poblacion.zorros= 40 #[número de zorros]
                       )

times <- seq(0, #periodo inicial [año]
             50, #periodo final [año]
             0.1)#incremento de tiempo [año]

intg.method<-c("euler")


zorros.conejos <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    
    poblacion.relativa<-poblacion.conejos/poblacion.sostenible.conejos #[Dmnl] 
  
    #variables auxiliares endógenas
    tasa.natalidad.zorros<-approx(c(0.00,0.50,1.00,1.50,2.0,1e6),
                                c(-0.12,0.12,0.25,0.27,0.3,0.3),xout = poblacion.relativa)$y # [1/año]
   

   
     caza.de.conejos.por.zorro <-approx(c(0.0,1.0,2.0,1e6),
                                        c(0.0,20,43,43), xout =poblacion.relativa)$y # [conejos/año/zorros]

     caza.de.conejos<-poblacion.zorros*caza.de.conejos.por.zorro
    
     pulso<- 0.0 #ifelse(t>10 & t<= 11,1,0)*100

    #variables de flujo
     flujo.nacimientos.conejos<-poblacion.conejos*tasa.natalidad.conejos+pulso #[conejos/año]
     flujo.muertes.conejos<-poblacion.conejos/vida.media.conejos+caza.de.conejos #[conejos/año]
     flujo.nacimientos.zorros<-poblacion.zorros*tasa.natalidad.zorros #[zorros/año]
     flujo.muertes.zorros<-poblacion.zorros/vida.media.zorros #[zorros/año]
     
    #variables de estado
     dpoblacion.conejos<-flujo.nacimientos.conejos-flujo.muertes.conejos
     dpoblacion.zorros<-flujo.nacimientos.zorros-flujo.muertes.zorros
    
     list(c(dpoblacion.conejos,dpoblacion.zorros),
          flujo.nacimientos.conejos=flujo.nacimientos.conejos,
          flujo.nacimientos.zorros=flujo.nacimientos.zorros,
          poblacion.relativa=poblacion.relativa,
          caza.de.conejos=caza.de.conejos,
          tasa.natalidad.zorros=tasa.natalidad.zorros,
          pulso = pulso,
          poblacion.relativa = poblacion.relativa
          
          
           )
  } )
}

#simulación del modelo
out <- ode(y = InitialConditions,
           times = times,
           func = zorros.conejos,
           parms = parameters,
           method =intg.method )

plot(out)

head(out)

class(out)

out1<-data.frame(out)

#Graficar el comportamiento
plot(out, which=c("dpoblacion.conejos","flujo.nacimientos.conejos"),xlab = "time", ylab =c("dpoblacion.zorros","flujo.nacimientos.zorros"))

#Graficar el comportamiento
plot(out, which=c("population","renewable.resources"),xlab = "time", ylab =c("","renewable resource units"))

#======================================================================================================================================

#Discutir como graficar diferentes corridas y diferentes variables 
#Graficar el comportamiento
plot(out, which=c("population","renewable.resources",
                  "resource.availability.dependent.lifetime"),xlab = "time", ylab =c("people","renewable resource units","lifetime"))


#ggplot2
library(ggplot2)
  ggplot(out1,aes(x=time, y=renewable.resources))+geom_line()+scale_y_continuous(sec.axis = out1$population )

head(out)