diff --git a/exploration/GR3D_Rdescription/GR3Dfunction.R b/exploration/GR3D_Rdescription/GR3Dfunction.R
index 2c2e43ef9ba1e2e70f875ab1f270d203a1bdf690..be56372edd4e25ead72996485e70cb279604b6ce 100644
--- a/exploration/GR3D_Rdescription/GR3Dfunction.R
+++ b/exploration/GR3D_Rdescription/GR3Dfunction.R
@@ -8,6 +8,15 @@ temperatureEffect = function(tempWater, Tmin, Topt, Tmax){
   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
@@ -27,12 +36,12 @@ vonBertalanffyInverse = function(L, L0, 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, deltaT, sigma, withTempEffect=FALSE, TrefAtSea = c(9.876946, 13.489854, 15.891487, 11.554104) ){
+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 * deltaT))) - sigma * sigma / 2
+    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]
@@ -42,13 +51,13 @@ vonBertalanffyIncrement = function(nStep, L0, Linf, K, deltaT, sigma, withTempEf
   return(L)
 }
 
-vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, deltaT, sigma, tempEffect ){
+vonBertalanffyWithNextIncrement = function(L, L0, Linf, K, timeStepDuration, sigma, tempEffect ){
   if (sigma == 0) {
-    mu = log((Linf - L) * (1 - exp(-K * deltaT)))
+    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration)))
     increment = exp(mu)
   }
   else {
-    mu = log((Linf - L) * (1 - exp(-K * deltaT))) - (sigma * sigma) / 2
+    mu = log((Linf - L) * (1 - exp(-K * timeStepDuration))) - (sigma * sigma) / 2
     increment = exp(rnorm(1, mu, sigma))
   }
   L = L + increment * tempEffect