# This R script calculates selected parameters of 
# plant root distribution pattern based on measured
# root length densities from different soil depths (cm).

# Author: Kuan Qin
#         Texas A&M AgriLife Research/Extension Ctr at Uvalde
#         1619 Garner Field Road, Uvalde, TX 78801
#         Email: qinkuan@tamu.edu
#         October 1, 2018
# ----------------------------------------------------------------
#  The accuracy of this program or its suitability for use must be 
#  judged by the user. The author disclaims any warranty and/or
#  liability with respect to the consequences of using this program. 
# ----------------------------------------------------------------

# Set your directory to your document location #
# Set soil id to integer numbers #
# Sort soil depth (sd) from each soil id descended, for example, 10 cm, 30 cm, ... , 90 cm #
dta <- read.csv("soil.csv", header = T)
attach(dta) 
d <- split(dta, factor(id))
for (i in 1:length(unique(id))){
	d[[i]]$sdnew[1] <- d[[i]]$sd[1]/2
	for (j in 2:length(d[[i]]$sd)){
		d[[i]]$sdnew[j] <- (d[[i]]$sd[j-1]+d[[i]]$sd[j])/2
	}
	d[[i]]$zr <- d[[i]]$sdnew/200
}
dta_KA <- data.frame()
dta_Lbottom <- data.frame()
dta_Ltop <- data.frame()
dta_nrld <- data.frame()
dta_B <- data.frame()

# Calculate K and alpha (A) #
for (i in 1:length(unique(id))){
	dta.i <-
	data.frame(K=summary(nls(data = d[[i]], rld ~ K*exp(-A*sdnew), start=list(K = 5,A = 0.01)))$parameters[1,1],
	K.SE=summary(nls(data = d[[i]], rld ~ K*exp(-A*sdnew), start=list(K = 5,A = 0.01)))$parameters[1,2],
	A=summary(nls(data = d[[i]], rld ~ K*exp(-A*sdnew), start=list(K = 5,A = 0.01)))$parameters[2,1],
	A.SE=summary(nls(data = d[[i]], rld ~ K*exp(-A*sdnew), start=list(K = 5,A = 0.01)))$parameters[2,2]) 
	dta_KA <- rbind(dta_KA, dta.i)
}

# Calculate Lbottom #
Lm <- 200
for (j in 1:length(unique(id))){
	Lp <- tail(d[[j]], 1)$sd
	dta.j <-
	data.frame(Lbottom=(-1)*dta_KA[j,1]/dta_KA[j,3]*(exp(dta_KA[j,3]*Lm*(-1))-exp(dta_KA[j,3]*Lp*(-1))))
	dta_Lbottom <- rbind(dta_Lbottom, dta.j)	
}

# Calculate Ltop #
for (j in 1:length(unique(id))){
	dt.1 <- d[[j]][1,]$sd*d[[j]][1,]$rld
		for (k in 2:length(unique(d[[j]]$sd))){
		dt.k <- (d[[j]][k,]$sd-d[[j]][k-1,]$sd)*d[[j]][k,]$rld
		dt.1 <- sum(dt.1, dt.k)
	}
	dta_Ltop <- rbind(dta_Ltop, dt.1)
}
colnames(dta_Ltop) <- c("Ltop")
dta_Ltotal <- dta_Lbottom + dta_Ltop

# Calculate normalized root length density for ith soil layer #
for (j in 1:length(unique(id))){
	nrld <-
	Lm * t(t(d[[j]]$rld)) / dta_Ltotal[j,]
	dta_nrld <- rbind(dta_nrld ,nrld)
}
colnames(dta_nrld) <- c("nrld")
dta_nrld <- as.matrix(dta_nrld)
dta$nrld <- dta_nrld
d1 <- split(dta, factor(id))
for (i in 1:length(unique(id))){
	d1[[i]]$zr <- d[[i]]$sdnew/200
}

# Calculate beta(B) #
for (i in 1:length(unique(id))){
	dta.i <-
	data.frame(B=summary(nls(data = d1[[i]], nrld ~ (1+B)*(1-zr)^B, start=list(B=2.6)))$parameters[1,1],
	B.SE=summary(nls(data = d1[[i]], nrld ~ (1+B)*(1-zr)^B, start=list(B=2.6)))$parameters[1,2])
	dta_B <- rbind(dta_B, dta.i)
}

# Calculate z.95 #
z.95 <- (1 - exp(log(0.05)/(dta_B[,1] + 1))) * Lm

# Finally, K, alpha and beta are stored in dta_KA and dta_B #
# Normalized root length density for ith soil layer is stored in dta #
# The absolute depth above which 95% roots are found is stored in z.95 #