development/RK4.R

In gholtgrieve/gassyPants: Functions Useful for Working With Gasses Dissolved in Water

#'@title
#'@description

#'@usage

#'
#'@param
#'@return
#'@note
#'@references
#'@author
#'@examples

#'@export


  #Function to do fourth-order Runge-Kutta integration
  #Calculated the value of state variable at next time step given value at this
  #time step, environmental forcing variables, and parameter values
  #CTS 16 May 2011 based on version from SRC

  #Note that in SRC version ('RKfix.R', there was more flexibility for different
  #time step sizes

  #Inputs
  # ode - Character string giving name of a function to calculate dy/dt
  #       This should return an output named "rate" that contains dy/dt;
  #       it may optionally return other outputs as well, which are ignored by RK
  # y0  - Value of state variable y at beginning of time step
  # E   - Vector containing environmental forcing variables required by ode
  # p   - Vector containing parameter values required by ode

  #Output
  # y1  - value of state variable y at end of time step

  RK4 <- function(y0, params, env) {
    TAlk <- params$TAlk
    PAR <- env$PAR
    temperature <- env$temperature
    salinity <- env$salinity
    pressure <- env$pressure

    K1 <- metabolism(TCO2_init = y0, params = params, env=env)
    k1 <- K1$TCO2

    K2 <- metabolism(TCO2_init = y0+(k1/2), params = params, env=env)
    k2 <- K2$TCO2

    K3 <- metabolism(TCO2_init = y0+(k2/2), params = params, env=env)
    k3 <- K3$TCO2

    K4 <- metabolism( TCO2_init = y0+k3, params = params, env=env)
    k4 <- K4$TCO2

    Y1 <- y0 + (k1/6) + (k2/3) + (k3/3) + (k4/6)
    CO2Sys <- get.CO2Sys.TAlkmolL_TCO2molL(TAlkmolL = TALk, TCO2molL = Y1,
                                           temperature = temperature, salinity = salinity,
                                           pressure = pressure)

    return(CO2Sys)
 }