## SimMat: a dynamic model for the prediction of inoculum production of L. maculans/L. biglobosa 
## species complex. Parameterisation and evaluation in French conditions.
## project of the ENDURE summer school Volterra 2016

## author: Stephanie
## creation date: 11/10/2016
## last update: 13/10/2016

## set directory:
setwd("~/courses/ENDURE")
#dataset
weather<-read.table("database_pseudothecia_ascospores_dull.csv", header=TRUE, sep=",", dec=".")
weather<-weather[c(367:731),] # there are NA, I chose to avoid them for now, plus the table is huge!

############################################################
# CREATION OF THE MODEL
############################################################

## 1. parameters
Tmax=22 #°C     (Salam et al., 2003)
Tmin= 6 #°C       (Salam et al., submitted)
rmin= 4 #mm     (Salam et al., 2003)
nd = 7 #days    (Salam et al., 2003)
NFD = 65 #days  (Aubertot et al., 2005)
SFD = 24 #days  (Aubertot et al., 2005)
k = 0.046 #mm-1 

SimMat= function(weather,Tmax, Tmin, rmin, nd, NFD, SFD, k)
  #weather is a table with at least a variable called "Temperature" and a variable called "Rainfall"
  {
    # 2. favourable days: cf power point for the equation
    # for day d>nd
  # I create a vector with only 0. Its length is equal to the number of days in the weather table
    DFM<-rep(0, dim(weather)[1])
      for (d in (nd+1):dim(weather)[1]) # I just consider the nd+1 days because of the rain factor wich requires the nd last days
      { rain_last_days<-sum(weather[c((d-nd):(d-1)), which(names(weather)=="Rainfall")])
      if(weather[d, which(names(weather)=="Temperature")]>Tmin & 
           weather[d, which(names(weather)=="Temperature")]<Tmax & rain_last_days>rmin){DFM[d]=1}
      }

    # 3.cumulative favorable days from the harvest (beginning of the simulation according to my code)
    CDFM<-rep(0, dim(weather)[1])
    for (d in (nd+1):dim(weather)[1])
    {CDFM[d]<-sum(DFM[1:d])}

    # 4. probability of the fungus maturation and spores release
    #4.1. the maximum proportion of ripe pseudotecia ready to sporulate from the harvest until the day d
    prop.full.pseudo.cumul<-pnorm(CDFM, mean=NFD, sd=SFD)
   
    prop.empty.pseudo<-rep(0,dim(weather)[1])
    prop.full.pseudo<-rep(0,dim(weather)[1])
    
    for (d in (nd+1):dim(weather)[1])
    {# 4.2. the ripe pseudotecia which has sporulated already at the day d
      prop.empty.pseudo[d]<-prop.empty.pseudo[d-1]+
        (prop.full.pseudo.cumul[d]-
           prop.empty.pseudo[d-1])*(1-exp(-k*weather[d, which(names(weather)=="Rainfall")]))
      # 4.3. the ripe pseudotecia which are ready to sporulate during the next rainfall
      prop.full.pseudo[d]<-(prop.full.pseudo.cumul[d]-prop.empty.pseudo[d])/(1-prop.empty.pseudo[d])  
    }
      
# combine everything in one table:
proba.dev.fungus<-data.frame(DFM, CDFM, prop.full.pseudo.cumul, prop.empty.pseudo, prop.full.pseudo)
return(proba.dev.fungus)
}

essai<-SimMat(weather,Tmax, Tmin, rmin, nd, NFD, SFD, k)
plot(essai[,5]~c(1:dim(weather)[1]), xlab="number of day since the harvest", ylab="proportion of full pseudotecia ready to sporulate")