GR3Dfunction.R

#----------------------------------------------
# Temperature effect
#----------------------------------------------
# temperature effect in GR3D
# from Rosso
temperatureEffect = function(tempWater, Tmin, Topt, Tmax){
  response = (tempWater - Tmin) * (tempWater - Tmax) / ((tempWater - Tmin) * (tempWater - Tmax) - (tempWater - Topt)^2)
  
  response[tempWater <= Tmin | tempWater >= Tmax] = 0
  return(response)
}

thermalRange = function(pct = 0.8, Tmin, Topt, Tmax){
  lower = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax) - pct,
           interval = c(Tmin, Topt))$root
  upper = uniroot(function(x) temperatureEffect(x, Tmin, Topt, Tmax) - pct,
                  interval = c(Topt, Tmax))$root
  return(c(lower = lower,upper = upper))
}

#optimalThermalRange(Tmin = 3, Topt = 17, Tmax =27)

#----------------------------------------------
# Growth simulation
#---------------------------------------------
# von Bertalaffy : L = f(age)
vonBertalanffyGrowth = function(age, L0, Linf, K){
  t0 = log(1 - L0 / Linf) / K
  return(Linf * (1 - exp(-K * (age - t0))))
}

# von Bertalanffy inverse  : age = f(L)
# with L0 rather the usual t0 parameter
vonBertalanffyInverse = function(L, L0, Linf, K){
  t0 = log(1 - L0/Linf)/K
  return(t0 - log(1 - L/Linf)/K)
}


# von Bertalanffy increment
# pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche
vonBertalanffyIncrement = function(nStep, L0, Linf, K, timeStepDuration, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
  tempEffect = temperatureEffect( TrefAtSea , 3, 17, 26)
  L = matrix(nrow = nStep + 1)
  L[1] = L0
  for (i in 1:nStep) {
    mu = log((Linf - L[i]) * (1 - exp(-K * timeStepDuration))) - sigma * sigma / 2
    increment = exp(rnorm(1, mu, sigma))
    if (withTempEffect) {
      increment = increment * tempEffect[((i - 1) %% 4) + 1]
    }
    L[i + 1] = L[i] + increment
  }
  return(L)
}

vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){
  if (sigma == 0) {
    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
    increment = exp(mu)
  }
  else {
    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
    increment = exp(rnorm(length(mu), mu, sigma))
  }
  L = L + increment * tempEffect
  return(L)
}

# -------------------------------------------------------