#====================================================================================
#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( 
              )

#Para definir tus condiciones iniciales utiiza tus variables de estado
InitialConditions <- c(
                       )

times <- seq(
  
            )

intg.method<-c("euler") #usaremos el método euler de integración


zorros.conejos <- function(t, state, parameters) {
  with(as.list(c(state,parameters)), {
    
    poblacion.relativa<-
  
    #variables auxiliares endógenas
     tasa.natalidad.zorros<-
     caza.de.conejos.por.zorro <-
     caza.de.conejos<-
    
     pulso<- 0.0 # ifelse(t>10 & t<= 11,1,0)*100  #Profero/a: explicar la función de "pulso" durante el proceso de modelación 

    #variables de flujo
    #especifica las variables de flujo
    
    
     
    #variables de estado
    #especifica las variables de estado. la "d" en población.conejos y población zorros señala que estás construyendo una ecuación diferencial
     dpoblacion.conejos<-
     dpoblacion.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 comportamiento
plot(out, which=c("dpoblacion.conejos","flujo.nacimientos.conejos"),xlab = "time", ylab =c("dpoblacion.zorros","flujo.nacimientos.zorros"))

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

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

#=====Profesor/a=====
#Discuta como graficar dir¡ferentes corridas y diferentes variables 
#Graficar comportamiento
plot(out, which=c("population","renewable.resources",
                  "resource.availability.dependent.lifetime"),xlab = "time", ylab =c("people","renewable resource units","lifetime"))

#Revisar cómo usar la ligrería ggplot2 para visualizar estas gráficas
library(ggplot2)
  ggplot(out1,aes(x=time, y=renewable.resources))+geom_line()+scale_y_continuous(sec.axis = out1$population )

head(out)