R/thermal_response_functions_basic.R

In lorecatta/DENVclimate:

briere <- function(T, c, T0, Tm) {

  N <- length(T)

  out <- rep(0, N)

  for (i in 1:N) {

    if(T[i] >= T0 & T[i] <= Tm) {

      out[i] <- c*T[i]*(T[i]-T0)*sqrt(Tm-T[i])

    }

  }

  out

}

quadratic <- function(T, c, T0, Tm) {

  N <- length(T)

  out <- rep(0, N)

  for (i in 1:N) {

    if(T[i] >= T0 & T[i] <= Tm) {

      out[i] <- -c*(T[i]-T0)*(T[i]-Tm)

    }

  }

  out

}

R_0 <- function(T, m, T_inf, c_a, T0_a, Tm_a, c_b, T0_b, Tm_b, c_c, T0_c, Tm_c,
                c_lf, T0_lf, Tm_lf, c_PDR, T0_PDR, Tm_PDR) {

  a <- briere(T,c_a,T0_a,Tm_a);
  b <- briere(T,c_b,T0_b,Tm_b);
  c <- briere(T,c_c,T0_c,Tm_c);
  lf <- quadratic(T,c_lf,T0_lf,Tm_lf);
  mu <- 1 / (lf+1);
  PDR <- briere(T,c_PDR,T0_PDR,Tm_PDR);
  eip <- 1 / PDR;

  message("===================================================================")
  message("a = ", a)
  message("b = ", b)
  message("c = ", c)
  message("lf = ", lf)
  message("mu = ", mu)
  message("PDR = ", PDR)
  message("===================================================================")

  (a^2*m*b*c*T_inf*(exp(-mu*eip)))/mu

}

prior_values <- function(n) {

  c_a <- rgamma(n, shape = 1, rate = 10)
  T0_a <- runif(n, min = 0, max = 24)
  Tm_a <- runif(n, min = 25, max = 45)
  c_b <- rgamma(n, shape = 1, rate = 10)
  T0_b <- runif(n, min = 0, max = 24)
  Tm_b <- runif(n, min = 25, max = 45)
  c_c <- rgamma(n, shape = 1, rate = 10)
  T0_c <- runif(n, min = 0, max = 24)
  Tm_c <- runif(n, min = 25, max = 45)
  c_lf <- rgamma(n, shape = 1, rate = 1)
  T0_lf <- runif(n, min = 0, max = 24)
  Tm_lf <- runif(n, min = 25, max = 45)
  c_PDR <- rgamma(n, shape = 1, rate = 10)
  T0_PDR <- runif(n, min = 0, max = 24)
  Tm_PDR <- runif(n, min = 25, max = 45)

  list(c_a = c_a,
       T0_a = T0_a,
       Tm_a = Tm_a,
       c_b = c_b,
       T0_b = T0_b,
       Tm_b = Tm_b,
       c_c = c_c,
       T0_c = T0_c,
       Tm_c = Tm_c,
       c_lf = c_lf,
       T0_lf = T0_lf,
       Tm_lf = Tm_lf,
       c_PDR = c_PDR,
       T0_PDR = T0_PDR,
       Tm_PDR = Tm_PDR)

}