# Calculating the allometric coefficient for the diameter height relationship, Ap - 16.05.2013. Updated on 15.02.2018.
rm(list = ls())

xPath <- "C:/Users/arng0001/Dropbox/Rcpp/"
# xPath <- "C:/Users/arng0001/Dropbox/Rcpp/"
# xPath <- "/Users/afs80b/Dropbox/Rcpp/"
# xPath <- "/Users/arnepommerening/Dropbox/Rcpp/"


options(digits = 14, width = 100)

# Select data 
# dataFile <- "Clg8_1"
# dataFile <- "Clg8_4"
# dataFile <- "Clg9_1"
# dataFile <- "Clg9_4"
# dataFile <- "Clg10_1"
# dataFile <- "Clg10_2" #*
#dataFile <- "Clg10_3"
# dataFile <- "Clg10_4"
# dataFile <- "Clg10_8"
# dataFile <- "Clg10_15"
# dataFile <- "Clg12_1"
# dataFile <- "Clg12_2"
# dataFile <- "Clg12_4"
# dataFile <- "Clg12_7"
# dataFile <- "Clg12_9"
# dataFile <- "Clg12_10"
# dataFile <- "Gwy4_3"
# dataFile <- "Gwy4_5"
# dataFile <- "Gwy4_10"
# dataFile <- "Gwy1_288"
dataFile <- "Gwy1_486"
# dataFile <- "Gwy1_548"
# dataFile <- "Gwy1_640"
# dataFile <- "GFS1_6"
# dataFile <- "GFS1_7"
# dataFile <- "GFS1_12"
# dataFile <- "Clg1_2000"
# dataFile <- "Clg1_3000"
# dataFile <- "Clg1_4000"
# dataFile <- "Clg1_5000"
# dataFile <- "Clg3_85"
# dataFile <- "Clg3_183"
# dataFile <- "Clg3_267"
# dataFile <- "Clg3_5000"
# dataFile <- "Clg3_5011"
# dataFile <- "Clg3_5100"
# dataFile <- "Clg3_5200"
# dataFile <- "ATW1_1947"



tdata <- read.table(paste(xPath, "StemAnalyses/", dataFile, ".txt", sep = ""), header = TRUE)
names(tdata)
tdata

# Make data set smaller
tdata$d <- tdata$dbh
tdata <- tdata[c("age", "h", "d")]

# Calculate observed relative height and diameter relative growth Rates after Wenk
tdata$ph <- 0
tdata$pd <- 0
for (i in 2 : length(tdata$age)) {
  tdata$ph[i] <- (tdata$h[i] - tdata$h[i - 1]) / tdata$h[i]
  tdata$pd[i] <- (tdata$d[i] - tdata$d[i - 1]) / tdata$d[i]
}

m2 <- log(1 - tdata$pd) / log(1 - tdata$ph)
mean(m2, na.rm = T)

# Calculate observed relative height and diameter RGRs
tdata$pd <- NA
tdata$ph <- NA
for (i in 2 : length(tdata$age)) {
	tdata$pd[i] <- (log(tdata$d[i]) -  log(tdata$d[i - 1])) / (tdata$age[i] - tdata$age[i - 1])
	tdata$ph[i] <- (log(tdata$h[i]) -  log(tdata$h[i - 1])) / (tdata$age[i] - tdata$age[i - 1])
}	


m1 <- tdata$pd / tdata$ph
mean(m1, na.rm = T)


# pdf(file = paste(file.path, "Gwy1_486.pdf", sep = ""))
par(mar = c(2, 2, 0.5, 0.5))
plot(m1, m2, las = 1, ylab = "", xlab = "", cex = .9, col = "black", pch = 19, axes = FALSE, ylim = c(0, 5), xlim = c(0, 5))  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
abline(0, 1, lty = 1, col = "red") 
box(lwd = 2)
# dev.off()

# pdf(file = paste(file.path, "Gwy1_486.pdf", sep = ""))
par(mar = c(2, 2, 0.5, 0.5))
plot(tdata$age, m1, las = 1, ylab = "", xlab = "", cex = .9, col = "black", pch = 19, axes = FALSE, ylim = c(0, 5), xlim = c(10, 90))  
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
# trend <- ksmooth(tdata$age[!is.infinite(m1)], m1[!is.infinite(m1)], kernel = "normal", bandwidth = 20)
trend <- ksmooth(tdata$age[-1], m1[-1], kernel = "normal", bandwidth = 20)
lines(trend, xlab = "", ylab = "", main = "", las = 1, lty = 1, lwd = 2, col = "red")
abline(h = 1, lty = 2) 
box(lwd = 2)
# dev.off()

# Linear regression to estimate the mean allometric coefficient
mregres <- lm(ph ~ pd - 1, data = tdata) # intercept = 0
# mregres <- lm(ph ~ pd, data = tdata) 
summary(mregres)
(b <- cov(tdata$ph[-1], tdata$pd[-1]) / var(tdata$ph[-1])) # slope for intercept not set to 0

# pdf(file = paste(file.path, "Gwy1_486.pdf", sep = ""))
par(mar = c(2, 4, 0.5, 0.5))
plot(tdata$pd, tdata$ph, las = 1, ylab = "", xlab = "", cex = .9, col = "black", pch = 19, axes = FALSE, ylim = c(0, 0.2), xlim = c(0, 0.2))  
curve(summary(mregres)$coefficients[1] * x, from = min(tdata$pd[-1]), to = max(tdata$pd[-1]), lwd = 2, lty = 1, col = "red", add = TRUE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
# abline(0, 1, lty = 1, col = "red") 
box(lwd = 2)
# dev.off()

mregreslog <- lm(log(h) ~ log(d), data = tdata) # Hunt (1990, p. 50)
summary(mregreslog)

# pdf(file = paste(file.path, "Gwy1_486.pdf", sep = ""))
par(mar = c(2, 3, 0.5, 0.5))
plot(log(tdata$d), log(tdata$h), las = 1, ylab = "", xlab = "", cex = .9, col = "black", pch = 19, axes = FALSE, ylim = c(0, 5), xlim = c(0, 5))  
curve(summary(mregreslog)$coefficients[1] + summary(mregreslog)$coefficients[2] * x, from = min(log(tdata$d)), to = max(log(tdata$d)), lwd = 2, lty = 1, col = "red", add = TRUE)
axis(1, lwd = 2, cex.axis = 1.8)
axis(2, las = 1, lwd = 2, cex.axis = 1.8)
# abline(0, 1, lty = 1, col = "red") 
box(lwd = 2)
# dev.off()