ODE forcings

## =============================================================================
## R-Code to solve example 3.3.2 from the book
## K. Soetaert, J. Cash and F. Mazzia, 2012.  
## Solving differential equations in R. UseR, Springer, 248 pp. 
## http://www.springer.com/statistics/computational+statistics/book/978-3-642-28069-6.
##
## An ordinary differential equation model with forcing functions
## The ozone model  - needs file "light.rda"
## implemented by Karline Soetaert
## =============================================================================

library (deSolve)

load(file = "light.rda")   # contains data.frame Light
head(Light, n = 4)

irradiance <- approxfun(Light)
irradiance(seq(from = 0, to = 1, by = 0.25))
k3 <- 1e-11; k2 <- 1e10; k1a <- 1e-30
k1b <- 1; sigma <- 1e11
yini <- c(O = 0, NO = 1.3e8, NO2 = 5e11, O3 = 8e11)

chemistry <- function(t, y, parms) {
  with(as.list(y), {

    radiation <- irradiance(t)
    k1  <- k1a + k1b*radiation

    dO   <-  k1*NO2 - k2*O
    dNO  <-  k1*NO2        - k3*NO*O3 + sigma
    dNO2 <- -k1*NO2        + k3*NO*O3
    dO3  <-           k2*O - k3*NO*O3
    list(c(dO, dNO, dNO2, dO3), radiation = radiation)
  })
}

times <- seq(from = 8/24, to = 5, by = 0.01)
out <- ode(func = chemistry, parms = NULL, y = yini, 
          times = times, method = "bdf")

plot(out, type = "l", lwd = 2, las = 1)