#Clear working environment
rm(list=ls())
#Load the needed library
library(deSolve)

#Define the differential equation
SOM_Diffeq <- function(t,y,params){
  L = y[1]; POM = y[2]; Mic = y[3]; MAOM = y[4];
  with(
    as.list(params),
    {
      dL <- S - r1*L - m1*L*Mic
      dPOM <- p1*r1*L - m2*POM*Mic - l1*POM
      dMic <- e1*m1*L*Mic + e2*m2*POM*Mic - d*Mic
      dMAOM <- p2*d*Mic - l2*MAOM
      res <- c(dL,dPOM,dMic,dMAOM)
      list(res)
    }
  )
}

#Define a time vector over which to simulate
time_vector=seq(from=0,to=200,by=1)

#Define parameter values
default_params=c(S=10,r1=.1,m1=.1,p1=.5,m2=.1,l1=.1,e1=1,e2=1,d=.05,p2=.1,l2=1)

#Put the simulation into some sort of object that can hold it
SOM_output1=rk(#Vector of starting values
                c(L=1,POM=1,Mic=1,MAOM=1),
                #Time vector
                time_vector,
                #Differential equation to simulate
                SOM_Diffeq,
                #parameters
                default_params,
                #My preferred method of using the rk function
                method="rk45dp7")

#Results take the form of a table with a number of rows equal to the length of the
#time vector and number of columsn equal to the number of players + 1. Column 1 is time,
#column 2 is the first player, column 3 is the second player, etc.

#Plot the results in the default manner
plot(SOM_output1)

#Perhaps we want to see the effect of changing a parameter, specifically the rate at
#which microbes metabolize litter: m1

#Define a range of m1 values to use
m1_range=seq(from=.05,to=1,by=.05)

#Create an array to hold the results from simulating many times
#This array has a number of rows equal to the number of m values used and a
#number of columns equal to the number of players
m1_effect=array(NA,dim=c(length(m1_range),4))

#Loop through and simulate for each m value
for(i in 1:length(m1_range)){
 m1_effect[i,]=rk(#Vector of starting values
                 c(L=1,POM=1,Mic=1,MAOM=1),
                  #Time vector
                  time_vector,
                  #Differential equation to simulate
                  SOM_Diffeq,
                  #Instead of using the same paramters every time, use a different value of
                  #m1 every time
                  c(S=10,r1=.1,m1=m1_range[i],p1=.5,m2=.1,l1=.1,e1=1,e2=1,d=.05,p2=.1,l2=1),
                  #My preferred method of using the rk function then select the players
                  #values at the end of the simulation to store in m1_effect
                  method="rk45dp7")[length(time_vector),2:5]
}

#Plot the results
#Tell R to plot 2 rows and 2 columns of graphs at a time
par(mfrow=c(2,2),mar=c(4,4,1,1))

#Plot the effect of m1 on the long-term amount of litter, POM, Mic, and MAOM
plot(
      #vector of x-axis values
      m1_range,
      #vector of y-axis values
      m1_effect[,1],
      #axis lavels
      xlab="microbial metabolism of litter: m1",ylab="Litter: L",
      #set the y window on the plot
      ylim=c(.9*min(m1_effect[,1]),1.1*max(m1_effect[,1])))

plot(m1_range,m1_effect[,2],xlab="microbial metabolism of litter: m1",ylab="POM",ylim=c(.9*min(m1_effect[,2]),1.1*max(m1_effect[,2])))
plot(m1_range,m1_effect[,3],xlab="microbial metabolism of litter: m1",ylab="Microbial biomass: Mic",ylim=c(.9*min(m1_effect[,3]),1.1*max(m1_effect[,3])))
plot(m1_range,m1_effect[,4],xlab="microbial metabolism of litter: m1",ylab="MAOM",ylim=c(.9*min(m1_effect[,4]),1.1*max(m1_effect[,4])))