Description Usage Arguments Author(s) Examples
Creates a Model14 object from different sources Have a look at the methods for details.
1 2 3 | GeneralModel_14(t, A, ivList, initialValF, inputFluxes, inputFc,
Fc, di = -0.0001209681, lambda = -0.0001209681, solverfunc = deSolve.lsoda.wrapper,
pass = FALSE, ...)
|
t |
|
A |
|
ivList |
|
initialValF |
|
inputFluxes |
|
inputFc |
|
Fc |
|
di |
|
lambda |
|
solverfunc |
The function used by to actually solve the ODE system. This can be |
pass |
if TRUE Forces the constructor to create the model even if it is invalid |
... |
Carlos A. Sierra, Markus Mueller
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))
|
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.