mkinmod <- function(...)
{
  spec <- list(...)
  obs_vars <- names(spec)

  # The returned model will be a list of character vectors, containing
  # differential equations, parameter names and a mapping from model variables
  # to observed variables
  parms <- vector()
  diffs <- vector()
  map <- list()

  # Establish list of differential equations
  for (varname in obs_vars)
  {
    if(is.null(spec[[varname]]$type)) stop("Every argument to mkinmod must be a list containing a type component")
    if(!spec[[varname]]$type %in% c("SFO", "FOMC", "SFORB")) stop("Available types are SFO, FOMC and SFORB only")
    new_parms <- vector()

    # New (sub)compartments (boxes) needed for the model type
    new_boxes <- switch(spec[[varname]]$type,
      SFO = varname,
      FOMC = varname,
      SFORB = paste(varname, c("free", "bound"), sep="_")
    )
    map[[varname]] <- new_boxes
    names(map[[varname]]) <- rep(spec[[varname]]$type, length(new_boxes))

    # Start a new differential equation for each new box
    new_diffs <- paste("d_", new_boxes, " =", sep="")

    # Construct terms for transfer to sink and add if appropriate
    if(is.null(spec[[varname]]$sink)) spec[[varname]]$sink <- TRUE
    if(spec[[varname]]$sink) {
      # Add first-order sink term to first (or only) box for SFO and SFORB models
      if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
        k_compound_sink <- paste("k", new_boxes[[1]], "sink", sep="_")
        sink_term <- paste("-", k_compound_sink, "*", new_boxes[[1]])
        new_diffs[[1]] <- paste(new_diffs[[1]], sink_term)
        new_parms <- k_compound_sink
      }
      if(spec[[varname]]$type == "FOMC") {
        if(match(varname, obs_vars) != 1) {
          stop("Type FOMC is only allowed for the first compartment, which is assumed to be the source compartment")
        }
        # From p. 53 of the FOCUS kinetics report
        fomc_term <- paste("(alpha/beta) * ((time/beta) + 1)^-1 *", new_boxes[[1]])
        new_diffs[[1]] <- paste(new_diffs[[1]], "-", fomc_term)
        new_parms <- c("alpha", "beta")
      }
    }
   
    # Add reversible binding if appropriate
    if(spec[[varname]]$type == "SFORB") {
      k_free_bound <- paste("k", varname, "free", "bound", sep="_")      
      k_bound_free <- paste("k", varname, "bound", "free", sep="_")      
      reversible_binding_terms <- c(
        paste("-", k_free_bound, "*", new_boxes[[1]], "+", k_bound_free, "*", new_boxes[[2]]),
        paste("+", k_free_bound, "*", new_boxes[[1]], "-", k_bound_free, "*", new_boxes[[2]]))
      new_diffs <- paste(new_diffs, reversible_binding_terms)
      new_parms <- c(new_parms, k_free_bound, k_bound_free)
    } 

    # Add observed variable to model
    parms <- c(parms, new_parms)
    names(new_diffs) <- new_boxes
    diffs <- c(diffs, new_diffs)
  }

  # Transfer between compartments
  for (varname in obs_vars) {
    to <- spec[[varname]]$to
    if(!is.null(to)) {
      origin_box <- switch(spec[[varname]]$type,
        SFO = varname,
        FOMC = varname,
        SFORB = paste(varname, "free", sep="_"))
      fraction_left <- NULL
      for (target in to) {
        target_box <- switch(spec[[target]]$type,
          SFO = target,
          SFORB = paste(target, "free", sep="_"))
        if(spec[[varname]]$type %in% c("SFO", "SFORB")) {
          k_from_to <- paste("k", origin_box, target_box, sep="_")
          diffs[[origin_box]] <- paste(diffs[[origin_box]], "-", k_from_to, "*", origin_box)
          diffs[[target_box]] <- paste(diffs[[target_box]], "+", k_from_to, "*", origin_box)
          parms <- c(parms, k_from_to)
        }
        if(spec[[varname]]$type == "FOMC") {
          fraction_to_target = paste("f_to", target, sep="_")
          fraction_not_to_target = paste("(1 - ", fraction_to_target, ")", sep="")
          if(is.null(fraction_left)) {
            fraction_really_to_target = fraction_to_target
            fraction_left = fraction_not_to_target
          } else {
            fraction_really_to_target = paste(fraction_left, " * ", fraction_to_target, sep="")
            fraction_left = paste(fraction_left, " * ", fraction_not_to_target, sep="")
          }
          diffs[[target_box]] <- paste(diffs[[target_box]], "+", fraction_really_to_target, "*", fomc_term)
          parms <- c(parms, fraction_to_target)
        }
      }
    }
  }
  model <- list(diffs = diffs, parms = parms, map = map)
  class(model) <- "mkinmod"
  invisible(model)
}