54: A general constructor

Description Usage Arguments Author(s) Examples

Description

Creates a Model14 object from different sources Have a look at the methods for details.

Usage

1
2
3
GeneralModel_14(t, A, ivList, initialValF, inputFluxes, inputFc, 
    Fc, di = -0.0001209681, lambda = -0.0001209681, solverfunc = deSolve.lsoda.wrapper, 
    pass = FALSE, ...)

Arguments

t
A
ivList
initialValF
inputFluxes
inputFc
Fc
di
lambda
solverfunc

The function used by to actually solve the ODE system. This can be deSolve.lsoda.wrapper or any other user provided function with the same interface.

pass

if TRUE Forces the constructor to create the model even if it is invalid

...

Author(s)

Carlos A. Sierra <[email protected]>, Markus Mueller <[email protected]>

Examples

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
t_start=1960
t_end=2010
tn=220
timestep=(t_end-t_start)/tn 
t=seq(t_start,t_end,timestep) 
n=3
At=new(Class="BoundLinDecompOp",
  t_start,
  t_end,
  function(t0){
        matrix(nrow=n,ncol=n,byrow=TRUE,
          c(-1,    0.1,    0, 
             0.5  , -0.4,    0,   
             0,    0.2,   -0.1)
        )
  }
) 
 
c0=c(100, 100, 100)
F0=ConstFc(c(0,10,10),"Delta14C")
#constant inputrate
inputFluxes=new(
  "TimeMap",
  t_start,
  t_end,
  function(t0){matrix(nrow=n,ncol=1,c(10,10,10))}
) 
# we have a dataframe representing the C_14 fraction 
# note that the time unit is in years and the fraction is given in
# the Absolute Fraction Modern format.
# This means that all the other data provided are assumed to have the same value
# This is especially true for the decay constants to be specified later
inputFc=BoundFc(C14Atm_NH,format="Delta14C")
# add the C14 decay to the matrix which is done by a diagonal matrix which does not vary over time
# we assume a half life th=5730 years
th=5730
k=log(0.5)/th #note that k is negative and has the unit y^-1

mod=GeneralModel_14(t=t,A=At,ivList=c0,initialValF=F0,inputFluxes=inputFluxes,inputFc=inputFc,di=k)
#start plots
par(mfrow=c(3,2))
   lt1=1;  lt2=2; lt3=3 
   col1=1;  col2=2; col3=3
   # plot the C and C14 curves
   Ct=getC(mod)
   plot(t,Ct[,1],type="l",lty=lt1,col=col1, ylim=c(0,200),
        ylab="C stocks (arbitrary units)",xlab="Time") 
   lines(t,Ct[,2],type="l",lty=lt2,col=col2) 
   lines(t,Ct[,3],type="l",lty=lt3,col=col3,lwd=2) 
   legend(
      "topright",
      c("C in pool 1",
        "C in pool 2",
        "C in pool 3"
      ),
      lty=c(lt1,lt2,lt3),
      col=c(col1,col2,col3)
   )
   C14t=getC14(mod)
   plot(t,C14t[,1]/1000,type="l",lty=lt1,col=col1, ylim=c(0,100),
        ylab="14C stocks (arbitrary units)",xlab="Time") 
   lines(t,C14t[,2]/1000,type="l",lty=lt2,col=col2) 
   lines(t,C14t[,3]/1000,type="l",lty=lt3,col=col3) 
   legend(
      "topright",
      c("14C in pool 1",
        "14C in pool 2",
        "14C in pool 3"
      ),
      lty=c(lt1,lt2,lt3),
      col=c(col1,col2,col3)
   )
   #now plot the C14 Fraction in the atmosphere and compute the C14/C fraction of in the Soil 

   FC14=getF14(mod)
   plot(C14Atm_NH, type="l",xlim=c(1960,2010))
   lines(t,FC14[,1],lty=lt1,col=col1) 
   lines(t,FC14[,2],lt2,type="l",lty=lt2,col=col2) 
   lines(t,FC14[,3],type="l",lty=lt3,col=col3) 
   legend("topleft",c(
                      expression(F[1]),
                      expression(F[2]),
                      expression(F[3])
                    )
   ,lty=c(lt1,lt2,lt3),col=c(col1,col2,col3))


#now compute the release flux
   Rt=getReleaseFlux(mod)
   plot(
     t,
     Rt[,1],
     type="l",
     lty=lt1,
     col=col1,
     ylab="C Release Flux (arbitrary units)",
     xlab="Time",
     ylim=c(0,50)
   ) 
   lines(t,Rt[,2],lt2,type="l",lty=lt2,col=col2) 
   lines(t,Rt[,3],type="l",lty=lt3,col=col3) 
   legend("topleft",c("RF1","RF2","RF3"),lty=c(lt1,lt2,lt3),col=c(col1,col2,col3))
   #now compute the c14 release flux
   R14t=getReleaseFlux14(mod)/1000
   plot(
     t,
     R14t[,1],
     type="l",
     lty=lt1,
     col=col1,
     ylab="C14 Release Flux (arbitrary units)",
     xlab="Time"
   ) 
   lines(t,R14t[,2],lt2,type="l",lty=lt2,col=col2) 
   lines(t,R14t[,3],type="l",lty=lt3,col=col3) 
   legend("topleft",c(
                      expression(RF[14]^1),
                      expression(RF[14]^2),
                      expression(RF[14]^3)
                    )
   ,lty=c(lt1,lt2,lt3),col=c(col1,col2,col3))

R14m=getF14R(mod)
C14m=getF14C(mod)
plot(C14Atm_NH, type="l",xlim=c(1960,2010),col=4)
lines(t,C14m) 
lines(t,R14m,col=2) 
legend(
  "topright",
  c("Atmosphere","Mean SOM-14C","Mean Release 14C"),
  lty=rep(1,3),
  col=c(4,1,2),
  bty="n"
)
par(mfrow=c(1,1))

SoilR documentation built on May 1, 2019, 8:06 p.m.