R/MeanAge.R

In SoilR: Models of Soil Organic Matter Decomposition

MeanAge=function
(IdotT,	
OdotLin,
sol,	
times 	
){
rho=function(a,tk){
	startTime=tk-a
	startVal=sol(startTime)
	if (startTime<0){
	        print("if")
		print(paste("a=",a,"tk=",tk))
		return(startVal)
	}
	times=c(startTime,tk,2*tk)
	s=solver(times,OdotLin,startVal)
	res=s[2]/sol(tk)
	return(res)
}
E=function(tk){
	E_integrand=function(Y,a){a*rho(a,tk)}
	res=solver(c(0,tk*0.99),E_integrand,0)[2]
return(res)
}
res=sapply(times,E)
return(res)
}
MeanAge2=function
(IdotT,	
OdotLin,
sol,	
times 	
){
   maxage=max(times)-min(times)
fineTimes=(seq(sqrt(min(times)),sqrt(max(times)),sqrt(maxage)/10000))^2   
startval=1
OofT=splinefun(fineTimes,solver(fineTimes,OdotLin,startval))
rho=function(a,tk){
	s=OofT(tk)*IdotT(tk-a)/OofT(tk-a)
	res=s/sol(tk)
	return(res)
}
E=function(tk){
	rho_tk=function(a){rho(a,tk)}
	E_integrand=function(Y,a){a*rho(a,tk)}
	res=solver(c(tk/1000,tk),E_integrand,0)[2]
return(res)
}
res=mcmapply(E,times,mc.cores=16)
return(res)
}
MeanAge3=function
(IdotT,	
OdotLin,
sol,	
times 	
){
so=solver(times,OdotLin,1)
deltaT=times[[2]]-times[[1]]
tstart=times[[2]] 
startvalue=sol(tstart)
if(startvalue==0){stop("cannot compute the solution since we are in a fixed point")}
Id=function(t){IdotT(t)}
l=100
pdf(file="runit.MeanAge.Densities.pdf",paper="a4r")
OofT=splinefun(times,so)
s=function(a,tk){
   OofT(tk)*Id(tk-a)/OofT(tk-a)
}
rho=function(a,tk){
	res=s(a,tk)/sol(tk)
	return(res)
}
E=function(tk){
        l=100
        ages=seq(0,tk-times[[2]],length=l)
	s_tk=function(a){s(a,tk)}
	ss_tk=splinefun(ages,mapply(s_tk,ages))
	one=integrate(ss_tk,lower=0,upper=tk-times[[2]])$value/sol(tk)
	rho_0=(1-one)/(deltaT)
        Integrand=function(a){a*rho(a,tk)} 
        spIntegrand=splinefun(ages,sapply(ages,Integrand))
	res=integrate(spIntegrand,lower=0,upper=max(ages))$value+rho_0*(tk-(deltaT)/2)
	return(res)
}
times[[1]]<-times[[2]]
dev.off()
res=mclapply(times,E,mc.cores=16)
return(res)
}