#final version 20210129

#script for a paper
#difference between hybrid larch and todo-matsu
#data from 2019 (june-september)
#sample trees L5 & A4

library(tidyverse)
library(lubridate)

##############################
Sys.setenv("TZ" = "Asia/Tokyo")

##############################
#data reading
L5 <- read_csv(
   "LDL5_2019.csv",
   col_types = "cddddddd"
 )
A4  <- read_csv(
  "LDA4_2019.csv",
  col_types = "cddddddd"
)
L5$DT <- as.POSIXct(L5$dt, tz = "Asia/Tokyo")
A4$DT <- as.POSIXct(A4$dt, tz = "Asia/Tokyo")
L5$D <- as_date(L5$DT, tz = "Asia/Tokyo")
A4$D <- as_date(A4$DT, tz = "Asia/Tokyo")
#-------
A4soil2019<-read_csv("A4soil_2019.csv", col_types = "cddd")
A4soil2019$DT<-as.POSIXct(A4soil2019$dt, tz = "Asia/Tokyo")
A4soil2019$D <- as_date(A4soil2019$DT, tz = "Asia/Tokyo")
vpd2019<-read_csv("VPD_2019.csv", col_types = "cdddd")
vpd2019$DT<-as.POSIXct(vpd2019$dt, tz = "Asia/Tokyo")
vpd2019$D <- as_date(vpd2019$DT, tz = "Asia/Tokyo")
hourly2019<-read_csv("hourly2019.csv", col_types = "cddddc")
hourly2019$DT<-as.POSIXct(hourly2019$dt, tz = "Asia/Tokyo")
hourly2019$D <- as_date(hourly2019$DT, tz = "Asia/Tokyo")

##############################
#calculating dLD
init_dt <- as.POSIXct("2019-06-12 00:00:00", tz = "Asia/Tokyo")#the initial starting point, as setting for A4 was at 11:50 12 June
init_A4_x <- A4$A4_LDx[A4$DT == init_dt]
init_A4_s <- A4$A4_LDs[A4$DT == init_dt]
init_L5_x <- L5$L5_LDx[L5$DT == init_dt]
init_L5_s <- L5$L5_LDs[L5$DT == init_dt]
init_A4_temp <- A4$A4_temp[A4$DT == init_dt]
init_L5_temp <- L5$L5_temp[L5$DT == init_dt]

A4$dLDx <- A4$A4_LDx - init_A4_x
A4$dLDs <- A4$A4_LDs - init_A4_s
A4$dT <- A4$A4_temp - init_A4_temp
L5$dLDx <- L5$L5_LDx - init_L5_x
L5$dLDs <- L5$L5_LDs - init_L5_s
L5$dT <- L5$L5_temp - init_L5_temp
#dLDs = dLDstem, dLDx = dLDxylem in the paper

##############################
#calculating dPP
#L5
Dhole<-80000 #80mm depth of increment borer for xylem
Lscrew<-105000 #105mm
Lrod<-26000 #26mm
alpha_sw<-27.1 * 10^-6 #Whethermax and Stamm
alpha_screw<-10.7 * 10^-6 #carbon steel rikanenpyo
alpha_rod<-16.0 * 10^-6 #stainless steel japan stainless steel association
offset<-25 #25mm
Rstem<-113 #113mm
Er<-359 #MPa, L5 Erad
Tsw<- 19000 #19 mm sapwood thickness
L5$dLD_tcorr<-L5$dLDx - (L5$dT*(Dhole*alpha_sw+Lrod*alpha_rod-Lscrew*alpha_screw))
L5$dLD_corr<-cos(asin(offset/Rstem))*L5$dLD_tcorr
L5$dPP<-L5$dLD_corr / Tsw * Er #MPa, pressure potential calculation

#A4
Dhole<-25000 #25mm depth of increment borer for xylem
Lscrew<-70000 #70mm
Lrod<-45000 #45mm
alpha_sw<-27.1 * 10^-6 #Whethermax and Stamm
alpha_screw<-10.7 * 10^-6 #carbon steel rikanenpyo
alpha_rod<-16.0 * 10^-6 #stainless steel japan stainless steel association
offset<-25 #25mm
Rstem<-88 #88mm
Er<-174 #MPa, A4 Erad 
Tsw<- 17000 #17 mm sapwood thickness
A4$dLD_tcorr<-A4$dLDx - (A4$dT*(Dhole*alpha_sw+Lrod*alpha_rod-Lscrew*alpha_screw))
A4$dLD_corr<-cos(asin(offset/Rstem))*A4$dLD_tcorr
A4$dPP<-A4$dLD_corr / Tsw * Er #MPa, pressure potential calculation

##############################
#overall range of dPP
L5_analy <- L5 %>% 
	filter(!is.na(dPP))	%>% 
	filter(D >= as_date("2019-06-12") & D <= as_date("2019-09-04")) 
max(L5_analy$dPP)#0.7759078
min(L5_analy$dPP)#-0.784395
	
A4_analy <- A4 %>% 
	filter(!is.na(dPP))	%>% 
	filter(D >= as_date("2019-06-12") & D <= as_date("2019-09-04")) 
max(A4_analy$dPP)#0.07166741
min(A4_analy$dPP)#-1.072432

##############################
#change into dailydata
makedaily <- function(x){
	x %>%
	select(DT, D, dLDx, dLDs, dPP) %>%
	filter(!is.na(dLDx) & !is.na(dLDs)) %>%
	filter(D >= as_date("2019-06-12") & D <= as_date("2019-09-04") ) %>%
	mutate(dLDb = dLDs-dLDx)%>% #dLDb = dBCX in the paper
	group_by(D) %>% summarise(
		dLDx.av = mean(dLDx),
		dLDs.av = mean(dLDs),
		dLDb.av = mean(dLDb),
		dLDx.rg = max(dLDx)-min(dLDx),
		dLDs.rg = max(dLDs)-min(dLDs),
		dLDb.rg = max(dLDb)-min(dLDb),
		dPP.av = mean(dPP),
		dPP.rg = max(dPP)-min(dPP)
	)
}
L5_daily <- makedaily(L5)
A4_daily <- makedaily(A4)

write.csv(L5_daily, "L5_daily.csv")
write.csv(A4_daily, "A4_daily.csv")

#daily data of climatic conditions
vpd_daily <- vpd2019 %>% group_by(D) %>% summarise(
	vpd.n = sum(!is.na(vpd_C2)),
	vpd.max = max(vpd_C2),
	vpd.min = min(vpd_C2),
	vpd.rg = max(vpd_C2) - min(vpd_C2),
	vpd.av = mean(vpd_C2)
) %>% filter(D >= as_date("2019-06-12") & D <= as_date("2019-09-04") ) 

prec_daily <- hourly2019 %>% group_by(D) %>% summarise(
	prec.sum = sum(prec)
)%>% filter(D >= as_date("2019-06-12") & D <= as_date("2019-09-04") ) 

soil_daily <- A4soil2019 %>% group_by(D) %>% summarise(
	vwc.av = mean(VWC_A4soil)
)%>% filter(D >= as_date("2019-06-12") & D <= as_date("2019-09-04") ) 

##############################
#combining daily data
#make complete set for daily data from 12 June to 4 September
dailyp1 <- left_join(vpd_daily, prec_daily, by = "D")
dailyp2 <- left_join(dailyp1, soil_daily, by = "D")
dailypL5 <- L5_daily %>%
	mutate(
		L5.bark.rg = dLDb.rg,
		L5.dPP.rg = dPP.rg
	) %>%
	select(D, L5.bark.rg, L5.dPP.rg)
dailypA4 <- A4_daily %>%
	mutate(
		A4.bark.rg = dLDb.rg,
		A4.dPP.rg = dPP.rg
	) %>%
	select(D, A4.bark.rg, A4.dPP.rg)
dailyp3 <- left_join(dailyp2, dailypL5, by = "D")
daily <- left_join(dailyp3, dailypA4, by = "D")

write.csv(daily, "daily_data.csv")

plot(daily$vpd.max, daily$vpd.rg)#confirmation for the relation max-range

##############################
#plot

##############################
#Fig1
par(mfrow=c(5,1))
par(mar=c(1,3,1,3))
s <- as.POSIXct("2019-06-01 00:00:00", tz = "Asia/Tokyo")
g <- as.POSIXct("2019-09-05 00:00:00", tz = "Asia/Tokyo")
plot(vpd2019$DT, vpd2019$temp_C2, type="n", xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(s,g), ylim=c(0, 40))
setwd(datadir)
source("abline_daily_2019.R")
abline(h=c(-10,0,10,20,30), col="grey")
lines(vpd2019$DT, vpd2019$temp_C2, col="orange")
axis(2, at=c(0,10,20,30,40), col="orange")
mtext("Crown temp, ºC", side=2, line=2, cex=0.7, col="orange")
par(new=T)
plot(hourly2019$DT, hourly2019$prec, type="h", axes=F, xlab="", ylab="", xlim=c(s,g), ylim=c(0, 40), col="blue")
axis(4, at=c(0,10,20,30,40), col="blue")
mtext("Precipitation, mm", side=4, line=2, cex=0.7, col="blue")

plot(vpd2019$DT, vpd2019$vpd_C2, type="n", xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(s,g), ylim=c(0, 4))
source("abline_daily_2019.R")
abline(h=c(0,1,2,3,4), col="grey")
lines(vpd2019$DT, vpd2019$vpd_C2, col="green")
axis(2, at=c(0,1,2,3,4), col="green")
mtext("VPD, kPa", side=2, line=2, cex=0.7, col="green")
par(new=T)
plot(A4soil2019$DT, A4soil2019$VWC_A4soil, type="l", axes=F, xlab="", ylab="", xlim=c(s,g), ylim=c(0, 1), col="black")
axis(4, at=c(0,0.5,1), col="black", line=0)
mtext("Soil VWC", side=4, line=2, cex=0.7, col="black")

plot(L5$DT, L5$dLD_corr, type="n", xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(s,g), ylim=c(-100, 1100))
source("abline_daily_2019.R")
abline(h=c(0, 500, 1000), col="grey")
lines(L5$DT, L5$dLDx, col="red")
lines(L5$DT, L5$dLDs, col="pink")
lines(A4$DT, A4$dLDx, col="blue")
lines(A4$DT, A4$dLDs, col="light blue")
axis(2, at=c(0,500,1000), col="black")
mtext("dLD, µm", side=2, line=2, cex=0.7, col="black")

plot(L5$DT, L5$dLD_corr, type="n", xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(s,g), ylim=c(-100, 100))
source("abline_daily_2019.R")
abline(h=c(-100, 0, 100), col="grey")
lines(L5$DT, L5$dLD_corr, col="red")
lines(A4$DT, A4$dLD_corr, col="blue")
axis(2, at=c(-100, 0, 100), col="black")
mtext("dLD_corr, µm", side=2, line=2, cex=0.7, col="black")

plot(L5$DT, L5$dPP, type="n", xaxt="n", yaxt="n", xlab="", ylab="", xlim=c(s,g), ylim=c(-1, 1))
source("abline_daily_2019.R")
abline(h=c(-1, 0, 1), col="grey")
lines(L5$DT, L5$dPP, col="red")
lines(A4$DT, A4$dPP, col="blue")
axis(2, at=c(-1, 0, 1), col="black")
mtext("dPP, MPa", side=2, line=2, cex=0.7, col="black")

##############################
#Fig2
par(mfrow = c(1,1))
boxplot(L5_daily$dLDs.rg, L5_daily$dLDx.rg, A4_daily$dLDs.rg, A4_daily$dLDx.rg, xaxt="n")
mean(L5_daily$dLDs.rg)
mean(L5_daily$dLDx.rg)
mean(A4_daily$dLDs.rg)
mean(A4_daily$dLDx.rg)

##############################
#Fig3
sd <- date("2019-06-11")
gd <- date("2019-09-05")
par(mfrow=c(5,1))
par(mar=c(1,4,1,4))
plot(daily$D, daily$vwc.av, typ="n", col="black", xlim=c(sd, gd), ylim=c(0.2,0.6), xaxt="n", yaxt="n", xlab="", ylab="")
source("abline_daily_2019_indailydata.R")
abline(h=c(0.2,0.3,0.4,0.5,0.6), col="grey")
lines(daily$D, daily$vwc.av)
axis(4, at=c(0.2,0.3,0.4,0.5,0.6))
mtext("Daily mean soil VWC", cex=0.7, side=4, line=2)
par(new=T)
plot(daily$D, daily$prec.sum, type="l", col="blue", xlim=c(sd, gd), ylim=c(0,40), axes=F, xaxt="n", yaxt="n", xlab="", ylab="")
axis(2, at=c(0,10,20,30,40), col="blue")
mtext("Daily precipitation, mm", cex=0.7, side=2, line=2, col="blue")

plot(daily$D, daily$vpd.rg, type="n", xlim=c(sd, gd), ylim=c(0,3), xaxt="n", yaxt="n", xlab="", ylab="")
source("abline_daily_2019_indailydata.R")
abline(h=0, lwd=2, col="grey")
abline(h=c(1,2,3), col="grey")
lines(daily$D, daily$vpd.rg, lty=1, col="green")
axis(2, at=c(0,1,2,3), col="green")
mtext("Daily range of VPD, kPa", side=2, line=2, col="green", cex=0.7)
par(new=T)
plot(daily$D, daily$L5.bark.rg, xlim=c(sd, gd), ylim=c(0,300), type="n", col="red", axes=F, xlab="", ylab="")
lines(daily$D, daily$L5.bark.rg, col="pink")
lines(daily$D, daily$A4.bark.rg, col="light blue")
axis(4, at=c(0,100,200,300), col="black")
mtext("Daily range of dBCX, µm", side=4, line=2, col="black", cex=0.7)

plot(daily$D, daily$L5.dPP.rg / daily$vpd.rg, type="n", col="red", xlim=c(sd, gd), ylim=c(0,400), xaxt="n", yaxt="n", xlab="", ylab="")
source("abline_daily_2019_indailydata.R")
abline(h=0, lwd=2, col="grey")
abline(h=c(100, 200, 300, 400), col="grey")
lines(daily$D, daily$L5.bark.rg / daily$vpd.rg, col="pink")
lines(daily$D, daily$A4.bark.rg / daily$vpd.rg, col="light blue")
axis(2, at=c(0, 100, 200, 300, 400), col="black")
mtext("dBCX range / VPD range", side=2, line=2, col="black", cex=0.7)

plot(daily$D, daily$vpd.rg, type="n", xlim=c(sd, gd), ylim=c(0,3), xaxt="n", yaxt="n", xlab="", ylab="")
source("abline_daily_2019_indailydata.R")
abline(h=0, lwd=2, col="grey")
abline(h=c(1,2,3), col="grey")
lines(daily$D, daily$vpd.rg, lty=1, col="green")
axis(2, at=c(0,1,2,3), col="green")
mtext("Daily range of VPD, kPa", side=2, line=2, col="green", cex=0.7)
par(new=T)
plot(daily$D, daily$L5.dPP.rg, xlim=c(sd, gd), ylim=c(0,1.5), type="n", col="red", axes=F, xlab="", ylab="")
lines(daily$D, daily$L5.dPP.rg, col="red")
lines(daily$D, daily$A4.dPP.rg, col="blue")
axis(4, at=c(0,0.5,1,1.5), col="black")
mtext("Daily range of dPP, MPa", side=4, line=2, col="black", cex=0.7)

plot(daily$D, daily$L5.dPP.rg / daily$vpd.rg, type="n", col="red", xlim=c(sd, gd), ylim=c(0,2), xaxt="n", yaxt="n", xlab="", ylab="")
source("abline_daily_2019_indailydata.R")
abline(h=0, lwd=2, col="grey")
abline(h=c(0.5,1,1.5,2), col="grey")
lines(daily$D, daily$L5.dPP.rg / daily$vpd.rg, col="red")
lines(daily$D, daily$A4.dPP.rg / daily$vpd.rg, col="blue")
axis(2, at=c(0,0.5,1,1.5,2), col="black")
mtext("dPP range / VPD range", side=2, line=2, col="black", cex=0.7)

##############################
#Fig4
#separated with terms
T1 <- as_date("2019-06-30")
T2 <- as_date("2019-07-06")
T3 <- as_date("2019-07-16")
T4 <- as_date("2019-08-02")
T5 <- as_date("2019-08-11")
T6 <- as_date("2019-08-20")
t1 <- daily %>% filter(D>=T1 & D<(T1+ddays(6)))
t2 <- daily %>% filter(D>=T2 & D<(T2+ddays(6)))
t3 <- daily %>% filter(D>=T3 & D<(T3+ddays(6)))
t4 <- daily %>% filter(D>=T4 & D<(T4+ddays(6)))
t5 <- daily %>% filter(D>=T5 & D<(T5+ddays(6)))
t6 <- daily %>% filter(D>=T6 & D<(T6+ddays(6)))
pre <- rbind(t1, t2, t3)
post <- rbind(t5, t6)

par(mfrow=c(1,1))
par(mar=c(3,3,3,3))
plot(t1$vpd.rg, t1$L5.dPP.rg, type="n", xlim=c(0,2.5), ylim = c(0,1), xaxt="n", yaxt="n", xlab="", ylab="")
abline(v=c(0,1,2), col="grey")
abline(h=c(0,0.2,0.4,0.6,0.8,1), col="grey")
points(daily$vpd.rg, daily$L5.dPP.rg, col="pink", pch=1)
points(daily$vpd.rg, daily$A4.dPP.rg, col="light blue", pch=1)
points(t1$vpd.rg, t1$L5.dPP.rg, col="red", pch=1)
points(t2$vpd.rg, t2$L5.dPP.rg, col="red", pch=4)
points(t3$vpd.rg, t3$L5.dPP.rg, col="red", pch=2)
points(t4$vpd.rg, t4$L5.dPP.rg, col="red", pch=7)
points(t5$vpd.rg, t5$L5.dPP.rg, col="red", pch=15)
points(t6$vpd.rg, t6$L5.dPP.rg, col="red", pch=16)
points(t1$vpd.rg, t1$A4.dPP.rg, col="blue", pch=1)
points(t2$vpd.rg, t2$A4.dPP.rg, col="blue", pch=4)
points(t3$vpd.rg, t3$A4.dPP.rg, col="blue", pch=2)
points(t4$vpd.rg, t4$A4.dPP.rg, col="blue", pch=7)
points(t5$vpd.rg, t5$A4.dPP.rg, col="blue", pch=15)
points(t6$vpd.rg, t6$A4.dPP.rg, col="blue", pch=16)
axis(1, at=c(0,1,2,3))
axis(2, at=c(0,0.4,0.8))
mtext("Daily range of VPD, kPa", side=1, line=2)
mtext("Daily range of dPP, MPa", side=2, line=2)

text(2.1, 0.9, "L5")
text(2.2, 0.9, "A4")
points(2.1,0.86, pch=1, col="pink" )
points(2.2,0.86, pch=1, col="light blue")
text(2.3,0.86, "All data", adj=0)
points(2.1,0.82, pch=1, col="red" )
points(2.2,0.82, pch=1, col="blue")
text(2.3,0.82, "T1", adj=0)
points(2.1,0.78, pch=4, col="red" )
points(2.2,0.78, pch=4, col="blue")
text(2.3,0.78, "T2", adj=0)
points(2.1,0.74, pch=2, col="red" )
points(2.2,0.74, pch=2, col="blue")
text(2.3,0.74, "T3", adj=0)
points(2.1,0.70, pch=7, col="red" )
points(2.2,0.70, pch=7, col="blue")
text(2.3,0.70, "T4", adj=0)
points(2.1,0.66, pch=15, col="red" )
points(2.2,0.66, pch=15, col="blue")
text(2.3,0.66, "T5", adj=0)
points(2.1,0.62, pch=16, col="red" )
points(2.2,0.62, pch=16, col="blue")
text(2.3,0.62, "T6", adj=0)

length(daily$vpd.rg)#85
cor(daily$vpd.rg, daily$L5.dPP.rg, use="complete.obs", method = "pearson")#0.4172233
length(daily$vpd.rg)#85
cor(daily$vpd.rg, daily$A4.dPP.rg, use="complete.obs", method = "pearson")#0.555724