Vezy_et_al_2018_Fig.A1

# Script to reproduce the figure A1: Tuzet et al. (2003) conductance model parameterization using 6 Coffea
# sprouts sap fluxes and leaf water potential. (1) Sf and Ψf parameters fit using data Ψmin and Ψmax data from (Dauzat et al., 2001) ;
# (2) Root-to-leaf plant conductivity (K) and (3) G0 and G1 parameters fitting.

# Importing packages ------------------------------------------------------

library(RColorBrewer)
CustomPalette= rev(brewer.pal(4, "RdYlBu")) # Not very nice for screen but Colorblind friendly (difficult for so many colors).
Linecolors= list(Measurement= CustomPalette[1],
                 Simulated= CustomPalette[3], 
                 CoffeePart= CustomPalette[4], 
                 Identity= "#636363", errors_bars= "#bdbdbd")



# Dummy data for example -------------------------------------------------

DataFluxes= data.frame(DeltaPsi= rnorm(100,mean = 0.5,0.5),
                       SapFF= rnorm(100,0.8,0.5),
                       Acagamfpsiv= rnorm(100,0.013,0.01))
DataFluxes$DeltaPsi[DataFluxes$DeltaPsi<0]= 0
DataFluxes$SapFF[DataFluxes$SapFF<0]= 0
DataFluxes$Acagamfpsiv[DataFluxes$Acagamfpsiv<0]= 0
DataFluxes$gs=0.0033+1.809*(DataFluxes$Acagamfpsiv*rnorm(100,1,0.1))
# Warning !! : replace with your own data !


# Finding G0 and G1 -------------------------------------------------------

Fit= lm(gs~Acagamfpsiv, data=DataFluxes)$coefficients

# Function to add or change alpha value of colours:
add.alpha <- function(col, alpha=1){
    # Function to add or change alpha value of colours
    # Source: https://gist.github.com/mages/5339689 ; Markus Gesmann
    if(missing(col))
        stop("Please provide a vector of colours.")
    apply(sapply(col, col2rgb)/255, 2,
          function(x)
              rgb(x[1], x[2], x[3], alpha=alpha))
}

# PsiF/Sf:

plotTuzet= function(){
    layout(matrix(c(1,2,3,3), nrow = 2, byrow = T), widths= c(1,2))
    par(oma= c(4,0.5,1.5,0.5), cex=1.5, bty = 'n',mar= c(1.5, 4, 1, 2))
    Psilmin= -2.3
    Psilmax= -1.3
    PsiF= median(c(Psilmin, Psilmax))
    Sf= 8
    PsiV_Tuz= seq((-3),0, 0.1)
    Fpsiv_Tuz= (1+exp(Sf*PsiF))/(1+exp(Sf*(PsiF-PsiV_Tuz)))
    
    plot(Fpsiv_Tuz~PsiV_Tuz, type='l', lwd=2, xlab="", ylab="", col=Linecolors[["Simulated"]])
    mtext(side=3, text=expression(Psi[v[max]]), at=-1.3)
    abline(v=-1.3, lwd=2)
    mtext(side=3, text=expression(Psi[v[min]]), at=-2.3)
    abline(v=-2.3, lwd=2)
    mtext(side=3, text=expression(Psi[italic(f)]), at=-1.8)
    abline(v=-1.8, lwd=2, lty=4, col='red')
    mtext(expression(f(Psi[italic(v)])), side = 2, line = 2.5, cex=1.5)
    mtext(expression(Psi[italic(v)]~(MPa)), side = 1, line = 3, cex=1.5)
    text(expression("" %->% Sf %->% ""), x= -1.99, y=0.3, srt= 74)
    text(x = -2.9, y = 1, labels = "1")
    par(mar= c(1.5, 4, 1, 0.5))
    # K:
    # NB: run 1-Finding_K.R to import table
    DataFluxes= DataFluxes[!is.na(DataFluxes$DeltaPsi),]
    
    plot(DataFluxes$DeltaPsi~DataFluxes$SapFF, ylab= "", xaxt= 'n', xlab="",
         pch= 21, col= Linecolors[["Measurement"]],xlim=c(0,4.5),
         bg= add.alpha(col= Linecolors[["Measurement"]], 0.5))
    mtext(expression(Delta~Psi~(MPa)), side = 1, line = 3, cex=1.5)
    mtext(expression(Sap~Flux~(mmol[H2O]~m[leaves]^-2~.s^-1)), side = 2, line = 2.5, cex=1.5)
    lines(I(0.7946*DataFluxes$DeltaPsi)~DataFluxes$DeltaPsi, col= Linecolors[["Simulated"]], lwd= 2)
    # abline(0,0.7946, col= Linecolors[["Simulated"]], lwd= 2)
    axis(side = 1, at= seq(0,4.5,0.5))
    text(expression(""%->% K %->% ""), x=  2, y= 2*0.7946*0.7946*1.35, srt= 24.5)
    text(x = 0, y = 2.5, labels = "2")
    
    # G0/G1:
    par(mar= c(0.5, 4, 1, 0.5))
    plot(DataFluxes$gs~DataFluxes$Acagamfpsiv, ylab= "", xaxt= 'n', xlab="",
         pch= 21, col= Linecolors[["Measurement"]],xlim= c(0,0.025),
         bg= add.alpha(col= Linecolors[["Measurement"]], 0.5))
    mtext(expression(A[CO2]%*%f(Psi[italic(v)])/(C[s]-Gamma)), side = 1, line = 3, cex=1.5)
    mtext(expression(G[s]~(mol[CO2]~m[leaves]^-2~.s^-1)), side = 2, line = 2.5, cex=1.5)
    # abline(Overall_G0_G1, col= Linecolors[["Simulated"]], lwd= 2)
    lines(I(Fit[1]+Fit[2]*DataFluxes$Acagamfpsiv)~DataFluxes$Acagamfpsiv, col= Linecolors[["Simulated"]], lwd= 2)
    axis(side = 1, at= seq(min(DataFluxes$Acagamfpsiv), 0.025, 0.005))
    text(expression(g0 %->% ""), x= -0.0004, y= Fit[1])
    text(expression(""%->% g1 %->% ""), x= 0.005, y= Fit[1]+Fit[2]*0.005*1.3, srt= 10.5)
    text(expression(G[s] == g0 + g1 %.% frac(A,c[i]-Gamma) %.% f(Psi[V])), x= 0.0025, y= 0.04)
    # text(expression(G[s] == 0.0033 + 1.809 %.% frac(A,c[i]-Gamma) %.% f(Psi[V])), x= 0.0025, y= 0.04)
    
    text(x = 0, y = 0.06, labels = "3")
    # legend("topleft", legend= paste0("Adjustment:","\n","g0= ", round(Fit3[1],5), "\n","g1= ", round(Fit3[2],3)),
    #        lwd=2, col=Linecolors[["Simulated"]], bty="n")
}

pdf(file= "Fig.A1_Tuzet_Parameterisation.pdf", width = 16, height = 9)
plotTuzet()
dev.off()

png(filename = "Fig.A1_Tuzet_Parameterisation.png", width = 16, height = 9, units = "in", res= 100)
plotTuzet()
dev.off()