nu = 51
sig = 0.0011
K=23490
mu=3.9

### calcul de Emid

#fonction pour l'intégrale double
f <- function(t,s){
  a <- (exp(s)-1)/s
  b <- (sig**2*t/6 + sig**2)*exp(-mu*t)*a
  c <-  nu*exp(-nu*t)*b
  return(c)
}

#fonction pour l'intégrale simple
g <- function(t){
  a <- nu*exp(-nu*t)
  b <- sig**2*t/3*(1-exp(-mu*t))
  return(a*b)
}


emid <- function(h){

  tmin <- h
  tmax <- 100
  smax <- function(x) mu*x
  smin <- 0

  I <- integral2(f,tmin,tmax,smin,smax)$Q

  J <- integrate(g,h,100)$value

  a <- (1 - exp(-nu*h)*(nu/(nu+h)))*exp(nu*h)

  return(a*(I+J))
}


#### calcul de E final

# fonction pour l'intégrale triple
k <- function(t,s,u){
  a <- (exp(u)-1)/u
  b <- 1/s
  c <- nu*exp(-nu*t)*exp(-mu*t)
  return(a*b*c)
}

efin <- function(h){

  tmin <- h
  tmax <- 100
  smax <- function(x) mu*x
  smin <- 0
  umax <- function(x,y) y
  umin <- 0

  I <- integral3(k,tmin,tmax,smin,smax,umin,umax)

  a <- (1-exp(-mu*h)*mu/(mu+h))
  b <- exp(mu*h)/a
  J <- nu**2*(nu*h+0.5)

  e <- emid(h) + (J/K**2)*b*I
  return(e)
}


#Graphique
c <- 0.1 * 1:140
d <- lapply(c,efin)

plot(c,d,xlab = 'h en secondes',ylab = 'E(h)',main = 'Courbe de E(h)')