#Proposta A

setwd("C:/Users/Mauro Brum-Jr/Documents/Arquivos_R")

#A funçao calcula parâmetros derivados da curva de pressão volume foliar em pesquisas ecofisiologicas
#  e agronomicas. A entrada water.potential representa o vetor de um data.frame com os valores de 
#  potencial hídrico. A entrada leaf.mass representa o vetor de um data.frame com os valores de 
#  perda de massa foliar. A entrada leaf.dry.mass, representa um vetor que representa a massa 
#  seca da folha de cada apos a finalizaçao das medidas. Nessa medidas o usuario usa a mesma folha
#  para fazer as medidas de potencial hídrico e concomitantemente a massa de água da folha. 

#Autor: Mauro Brum Monteiro Junior
#e-mail: maurobrumjr@gmail.com
#data: 18/05/2014
#Versão 1

data.curve=read.table("folha3.txt", header=TRUE, sep="\t")
pv.curve(data.curve$water.potential,data.curve$leaf.mass,data.curve$leaf.dry.mass)

pv.curve<-function(wp,lmass,ldmass)
{
  water.potential=wp
  leaf.mass=lmass
  leaf.dry.mass=ldmass

### Verifica classe de vetores, pois todos os vetores devem ser numericos. 
if(!is.numeric(water.potential) & !is.integer(water.potential)){
  stop("Vetor 'water.potential' nao e 'numeric' ou 'integer'\n",call.=FALSE);
}
if(!is.numeric(leaf.mass) & !is.integer(leaf.mass)){
  stop("Vetor 'leaf.mass' nao e 'numeric' ou 'integer'\n",call.=FALSE);
}
if(!is.numeric(leaf.dry.mass) & !is.integer(leaf.dry.mass)){
  stop("Vetor 'leaf.dry.mass' nao e 'numeric' ou 'integer'\n",call.=FALSE) 
}

### A funçao irá remover NA automaticamente.
qNA1<-length(water.potential)-length(na.omit(water.potential))
qNA2<-length(leaf.mass)-length(na.omit(leaf.mass))
qNA3<-length(leaf.dry.mass)-length(na.omit(leaf.dry.mass))
#para o calculo da quantidade de NA eu usei o comprimento de vetor menos o comprimento do
#vetor gerado pela omissão de NA. A diferença e a quantidade de NA. Se a diferença for diferente
#de zero é porque tem NA.

if(qNA1!=0){
  water.potential<-as.vector(na.exclude(water.potential) )
  cat("\n\t Foi removido",qNA1,"NA´s no vetor water.potential\n
       \n\t Se foi removido o NA apenas da variável water.potential colocar NA da linha
       \n\t respectiva da coluna leaf.mass e leaf.dry.mass do data.frame original\n\n")}
if(qNA2!=0){
  leaf.mass<-as.vector(na.exclude(leaf.mass)) 
  cat("\n\t Foi removido",qNA2,"NA´s no vetor leaf.mass\n
       \n\t Se foi removido o NA apenas da variável leaf.mass colocar NA da linha
       \n\t respectiva da coluna water.potential e e leaf.dry.mass do data.frame original\n\n")}  
if(qNA3!=0){
  leaf.dry.mass<-as.vector(na.exclude(leaf.dry.mass))
  cat("\n\t Foi removido",qNA3,"NA´s no vetor leaf.dry.mass\n
      \n\t Se foi removido o NA apenas da variável leaf.mass colocar NA da linha
      \n\t respectiva da coluna water.potential e e leaf.mass do data.frame original\n\n")}

#Acima uma sequecia de comando usando a funcao if para remover NA se as diferenças calculadas 
#acima forem diferentes de 0. 

if(qNA1==0 & qNA2==0 & qNA3==0){
  cat("\n\n\t Nenhum NA foi removido nos vetores\n\n\n")}
#Acima um comando para retornas a mensagem que nenhum NA foi removido quando as diferenças calculadas
#acima forem iguais a 0.

### Calcular e criar um vetor com dados do conteúdo relativo de água
  ldm<-data.curve$leaf.dry.mass[1] #fixa um valor do vetor da massa seca foliar (leaf.dry.mass)
#Calculo baseado na fórmula RWC=(massa folha-massa seca/massa turgida-massa seca)*100  
  rwc<-((leaf.mass-ldm)/(max(leaf.mass)-ldm))*100 
#Calculo para determinar a perda de água relativa de água  
  water.loss<- 100-rwc
#Calculo para determinar o inverso   
  inverse.w.potential<- -1/water.potential
#Calculo da massa de água contida na folha  
  w.mass<-leaf.mass-ldm
#Criando um data.frame com os vetores criado até o momento, pois o usuário pode recuperar os data.frame  
  data.curve<-data.frame(water.potential,leaf.mass, leaf.dry.mass, 
                         inverse.w.potential, water.loss, rwc, w.mass)

###Parameter asymptotic exponential
#Estimadores A,B,C por meio de estimadores self.starting SSasymp(). Colsulte pagina 675 Crawley, 2007.
  modelo<-nls(inverse.w.potential~SSasymp(water.loss,a,b,c))
###Gráfico  
  plot(inverse.w.potential~water.loss, data=data.curve, pch=20, main="blue dot=Po red dot=TLP or PPT",
       cex.main=0.8)
  xv<-seq(min(water.loss), max(water.loss)) 
  yv<-predict(modelo, list(water.loss=xv))
  lines(xv,yv, col="red", lty=3, lwd=2 )
#acima criei um uma sequencia de valor xv para estimar y em funçao do modelo  
 
#Usar o valor [a] do modelo encontrado anteriormente como o ponto em que a curva começa a formar
  #uma reta para encontrar o PPT e calcular a equaçao da reta.

### Encontrar uma equaçao da reta. O objetivo é encontrar uma equaçao da linear com os pontos que
# estão após o ponto de inflexão em que a curva comeca a virar uma reta. Esse ponto é o ponto previsto
# como o ponto em que as celular da folha perde a pressão da parede celular. Pp=0, depois desse ponto o 
# componente mais importante do potencial hidrico da folha e o potencial osmotico.

# O modelo nls gera os coeficiente a,b e c. O valor a é o resultado estimado do ponto de de assintota
# e determinará o ponto de perda de turgor. 
  coef.model<-coef(modelo)
# Encontrar o desvio do valor a estimado pelo modelo nls pela funçao profile()
  get.sd<-profile(modelo)
#o objeto é uma lista contendo 3 arquivos. Tenho interesse apenas no abjeto que está atribuido 
#a essa lista summary, isto é o desvio padrão do valor A
  o.summary<-attr(get.sd, "summary")
  o.parameter<-as.data.frame(o.summary[[11]])
  sd.a<-o.parameter[1,2]
#Somei o devido padrão ao coeficiente a para determinar o ponto de perdar de turgor do eixo Y.
a.value<-as.vector(coef.model[1]+sd.a) #equivale ao ponto de assintota predito pelo modelo.

#Encontra o valor de X predito pelo modelo de assintota exponencial
#Tentei usar a funçao predic(). Mas ela não estimou o valor x no local correto da curva. Então
#usei os valores de X estimado no na parte do gráfico (ver comando line()) para encontrar o valor x.. 
  pred.x.y<-data.frame(xv,yv)
  table.xy<-pred.x.y[pred.x.y$y<=a.value, c("xv","yv")]
  x.point<-table.xy[1,]
  x.value<-x.point$xv
#coloca um ponto vermelho no gráfico onde está localizado o ponto de perda de  turgor PPT.
  points(x.value,a.value, pch=20, cex=2, col="Red")
  
  
#Todos os calculos daqui para frente estao relacionados com os dados que estão antes e depois do PPT.  
#Dados abaixo do ponto de perda de turgor
  under.ppt<-data.curve[data.curve$water.loss>=x.value,
                        c("water.potential","leaf.mass","inverse.w.potential","rwc","water.loss", "w.mass")]
#Acrecentarei os valores ppt previsto pelo modelo os valores x e y dos dados originais. 
#Isso foi feito para considerar os valores de ppt previsto na parte linear apos o PPT.
w.lm.x<-c(x.value,under.ppt$water.loss)
iwp.lm.y<-c(a.value, under.ppt$inverse.w.potential)
#O intercepto representa o ponto em que o potencial osmotico esta hidrataçao máxima. Será usado para
#definir um parametro Po.
y.point<-as.vector(coef(lm(iwp.lm.y~w.lm.x))[1])
#representa o coeficiente angular da reta
a.point<-as.vector(coef(lm(iwp.lm.y~w.lm.x))[2])
#mostra o ponto Po no gráfico.
points(0,y.point,pch=20,col="blue", cex=2)
#estabelce a linha predida pelo modelo linear.
abline(lm(iwp.lm.y~w.lm.x))


##############Calculos parametros derivados da curva###########
  
  ###potencial osmótico na hidrataçao máxima (MPa). O inverso do intercepto y da relaçao linear.
  Po<-as.vector(-1/y.point)
  ###potencial de perda de turgor (MPa). O inverso do a valor determinado pelo modelo nls
  PPT<-as.vector((-1/a.value))
  ###Calculo do potencial osmótico (MPa). Calculado 
  data.curve$p.sol<-p.sol<-(-1/(y.point+a.point*data.curve$water.loss)) 
  ###Calculo do potencial de pressão (MPa)
  data.curve$p.press<-p.press<-(data.curve$water.potential)-(data.curve$p.sol)
  
  
  ###Separar os valores que estão acima do PPT
  above.ppt<-data.curve[data.curve$water.loss<=x.value, 
                        c("water.potential","leaf.mass","rwc","water.loss", 
                          "inverse.w.potential","p.sol", "p.press", "w.mass")]
  
  #PLATEAU EFFECT & ESTIMATION OF SATURATED WATER CONTENT (SWC)
  #Para estimar o conteudo saturado de água é necessário fazer uma modelo linear da massa de água 
  #em funçao do potencial hídrico. Usaremos os coeficentes dessa relaçao em outros parâmetros. 
  wm.wp<-lm(w.mass~water.potential, data=above.ppt)
  b.wm.wp<-as.vector(coef(lm(w.mass~water.potential, data=above.ppt))[1])
  a.wm.wp<-as.vector(coef(lm(w.mass~water.potential, data=above.ppt))[2])
  #estimando a proporçao saturada de água baseada no intercepto da funçao w.mass~water.potential
  #tanto para os valore acima e abaixo do PPT.
  above.ppt$swc.e<-above.ppt$w.mass/b.wm.wp
  under.ppt$swc.e<-under.ppt$w.mass/b.wm.wp
  
  ###Calculo do coeficiente de elásticidade (MPa)
  #Por definiçao o coeficiente de elasticidade é o slope da inclinaçao da funçao entre 
  #o potencial hídrico e o conteúdo saturado de água dos dados acima do ponto de perda de turgor.
  elast<-as.vector(coef(lm(p.press~swc.e, data=above.ppt))[2])
  #retoma o coeficiente de elasticidade (MPa)
  
  ###Satured water contend (g/g). Calculado pelo intercepto da relaçao 
  #linear funçao w.mass~water.potential 
  SatWcontent<-as.vector(b.wm.wp/ldm)
  
  ###Fraçao simplástica de água: quantidade de água no citoplásma (g)
  SympWFrac<-as.vector((((-y.point)/a.point)/100)*y.point)
  
  ###RWCppt: Conteudo relativo de água no PPT (%)
  RWCtlp<-(under.ppt$swc.e[1])*100
  
  ###Capacitancia antes do ponto de perda de turgor.
  #Corresponde ao coeficiente angular da relaçao swc.e~water.potential abaixo do ponto de perda
  #de turgor (1/MPa)
  Cft<-as.vector(coef(lm(swc.e~water.potential, data=above.ppt))[2])
  ###Capacitancia depois do ponto de perda de turgor (1/MPa)
  Ct<-as.vector(coef(lm(swc.e~water.potential,data=under.ppt))[2])
  ###Capacitancia absoluta (mol/kg.MPa).
  Cabs<-as.vector(Cft*SatWcontent/18*1000)
  #ps: 1 mol de H2O = 18 g
  ###Water extracted between full turgor and turgor loss point (mol/kg)
  Wt<-(SatWcontent*(1-RWCtlp/100))/18*1000
  ###Water extracted between gravitational potential and turgor loss point
  Wgt<-((a.wm.wp*(-0.01)*ldm+b.wm.wp)-(RWCtlp/100*b.wm.wp))/ldm/18*1000
  ###Concentraçao de solutos osmóticamente ativos
  Ns<-((SympWFrac*Po/(8.314462*294.26))*1000)*(-1)
  
  ###Gerando o data.frame
  par.pvcurve<-c("Po","PPT", "elast", "SatWcontent", "SympWFrac", 
                 "RWCppt","Cft","Ct","Cabs","Wt","Wgt", "Ns")
  res.pvcurve<-round(c(Po,PPT,elast,SatWcontent,SympWFrac,RWCtlp,Cft,Ct,Cabs,Wt,Wgt,Ns),4)
  unidade<-c("MPa","MPa","MPa","g/g","g","%","1/MPa","1/MPa",
             "mol/kg.MPa","mol/kg", "mol/kg","osmol*1000")
  result.1<-data.frame(par.pvcurve,res.pvcurve,unidade)
  
  return(result.1)
}