Merge branch 'patrick' into Camille_Bis

.gitignore

 <<<<<<< HEAD
-output
-.settings
-.classpath
-.project
-
+output
+.settings
+.classpath
+.project
+
 /bin/
 /target/
 =======
@@ -111,6 +111,9 @@ org.*
 /outcmaesxmean.dat
 /outcmaesxrecentbest.dat
 /exploration/scriptR/BruchPrg.R.orig
+.Rproj.user
+/exploration/scriptR/.RData
+/exploration/scriptR/*.Rproj
 /exploration/scriptR/.Rhistory
 /exploration/scriptR/scriptR.Rproj
 /exploration/scriptR/.Rproj.user/

data/input/fishTryRealBV_CC.xml

@@ -249,7 +249,7 @@
 
 				<species.PopulateBasinNetworkSeveralTimesAccordingToBasinSize>
 					<synchronisationMode>ASYNCHRONOUS</synchronisationMode>
-					<nbFishPerSI>2500</nbFishPerSI>
+					<nbFishPerSI>100000</nbFishPerSI>
 					<initialLength>2.0</initialLength>
 					<etaPopulate>40.0</etaPopulate>
 					<timesOfPopulate>5</timesOfPopulate>
@@ -265,8 +265,8 @@
 					<tempMinGrow>3.0</tempMinGrow>
 					<tempMaxGrow>26.0</tempMaxGrow>
 					<tempOptGrow>17.0</tempOptGrow>
-					<kOptForFemale>0.29</kOptForFemale>
-					<kOptForMale>0.21</kOptForMale>
+					<kOptForFemale>0.3236</kOptForFemale>
+					<kOptForMale>0.2141</kOptForMale>
 					<sigmaDeltaLVonBert>0.2</sigmaDeltaLVonBert>
 				</species.Grow>
 

exploration/scriptR/BruchPrg.R

@@ -60,64 +60,63 @@ temperatureEffect= function(temp, Tmin, Topt, Tmax){
 # ------------------------------------------------
 # temperature effect on growth
 # ------------------------------------------------
-tempData=read.csv("/home/patrick/Documents/workspace/GR3D/data/input/reality/SeasonTempBVFacAtlant1801_2100_newCRU_RCP85.csv", sep=";")
-sel = tempData$NOM=="Garonne" & tempData$Year>=2008 & tempData$Year<=2018
-Tref=colMeans(tempData[sel, c("Winter", "Spring", "summer", "Autumn")])
+tempData = read.csv("/home/patrick/Documents/workspace/GR3D/data/input/reality/SeasonTempBVFacAtlant1801_2100_newCRU_RCP85.csv", sep =";")
+sel = tempData$NOM == "Garonne" & tempData$Yea >= 2008 & tempData$Year <= 2018
+Tref = colMeans(tempData[sel, c("Winter", "Spring", "summer", "Autumn")])
 TrefAtSea = (12 + Tref) / 2
 
-temperature=seq(0,30,.1)
+temperature = seq(0,30,.1)
 # temperature effect on growth
 tempEffects = temperatureEffect(temperature,  3, 17, 26)
-plot(temperature,tempEffects, type='l')
+plot(temperature,tempEffects, type = 'l')
 
-points(TrefAtSea,  temperatureEffect(TrefAtSea, 3, 17, 26), col="red", pch =20)
-text(TrefAtSea, temperatureEffect(TrefAtSea, 3, 17, 26),  c("Winter", "Spring", "Summer", "Autumn"), pos=4)
+points(TrefAtSea,  temperatureEffect(TrefAtSea, 3, 17, 26), col = "red", pch = 20)
+text(TrefAtSea, temperatureEffect(TrefAtSea, 3, 17, 26),  c("Winter", "Spring", "Summer", "Autumn"), pos = 4)
 
-points(Tref,  temperatureEffect(Tref, 3, 17, 26), col="red")
-text(Tref, temperatureEffect(Tref, 3, 17, 26),  c("Winter", "Spring", "Summer", "Autumn"), pos=2, cex=.5)
+points(Tref,  temperatureEffect(Tref, 3, 17, 26), col = "red")
+text(Tref, temperatureEffect(Tref, 3, 17, 26),  c("Winter", "Spring", "Summer", "Autumn"), pos = 2, cex = .5)
 
-plot(tempData$Year[sel], tempData$Winter[sel], type='l')
+plot(tempData$Year[sel], tempData$Winter[sel], type = 'l')
 
 # ----------------------------------------------
 # growth simulation
 # ----------------------------------------------
 vonBertalaffyGrowth = function(age, L0, Linf, K){
-  t0=log(1-L0/Linf)/K
-  return(Linf*(1-exp(-K*(age-t0))))
+  t0 = log(1 - L0 / Linf) / K
+  return(Linf * (1 - exp(-K * (age - t0))))
 }
 
 vonBertalaffyInverse = function(L, L0, Linf, K){
-  t0=log(1-L0/Linf)/K
-  return (t0 - log(1-L/Linf)/K)
+  t0 = log(1 - L0/Linf)/K
+  return(t0 - log(1 - L/Linf)/K)
 }
 
-KvonBertalaffy = function(age, L, L0, Linf){
-  t0=log(1-L0/Linf)/K
-  return (-log(1-L/Linf)/(age-t0))
+KvonBertalaffy = function(age, L, L0, Linf) {
+  return(-log((Linf - L ) / (Linf - L0)) / age)
 }
 
 #vonBertalaffyInverse(L=40, L0, Linf, koptMale)
 
-Zfct=function(age, Ltarget,  L0, Linf, K, Tref ) {
-  return (sapply(age,function(a) (vonBertalaffyGrowth(a, L0, Linf, K) * mean(temperatureEffect(Tref, 3, 17, 26)) - Ltarget)^2))
+Zfct = function(age, Ltarget,  L0, Linf, K, Tref ) {
+  return(sapply(age,function(a) (vonBertalaffyGrowth(a, L0, Linf, K) * mean(temperatureEffect(Tref, 3, 17, 26)) - Ltarget)^2))
 }
 
-Pauly= function(age, t0, Linf, K, D){
-  return(Linf/10*((1-exp(-K*D*(age-t0)))^(1/D)))
+Pauly = function(age, t0, Linf, K, D){
+  return(Linf / 10 * ((1 - exp(-K * D * (age - t0)))^(1 / D)))
 }
 
 # pas cohérent avec la temperature effet sur le coeff de croissance mais ca marche
 vonBertalaffyIncrement = function(nStep, L0, Linf, K, deltaT, 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
+  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 * deltaT))) - sigma * sigma / 2
     increment = exp(rnorm(1, mu, sigma))
-    if (withTempEffect){
-      increment = increment * tempEffect[((i-1) %% 4)+1]
+    if (withTempEffect) {
+      increment = increment * tempEffect[((i - 1) %% 4) + 1]
     }
-    L[i+1]=L[i]+increment
+    L[i + 1] = L[i] + increment
   }
   return(L)
 }
@@ -126,17 +125,17 @@ vonBertalaffyIncrement = function(nStep, L0, Linf, K, deltaT, sigma, withTempEff
 # parametres
 # ----------------------------------------------------
 L0 = 2
-Linf = 80
+Linf = 70
 #kopt = 0.3900707
 #koptMale = 0.2
 #koptFemale = 0.3
-timestep = 0.25
+timestep = 0.25 # time step of the simumation
 sigma = 0.2
 
 tempEffect = mean(temperatureEffect( TrefAtSea, 3, 17, 26))
 
 #c(mean(temperatureEffect( Tref, 3, 17, 26)), mean(temperatureEffect( TrefAtSea, 3, 17, 26)))
-age=seq(0,8,.25)
+age = seq(0,8,.25)
 
 
 # ###########################################################################"
@@ -203,8 +202,9 @@ extractTemp<-as.data.frame(tempData %>%
 )
 
 plot(extractTemp$Lat, extractTemp$meanSpring, pch=20)
+abline(v=45.25) # Garonne
 abline(v=49.25) # vire
-abline(v=53.75) # elbe
+abline(v=52.25) # Rhine
 
 
 abline(h=11)
@@ -212,7 +212,7 @@ abline(h=12.5)
 # ======================================================================
 # exploration of growth for male and female
 # ======================================================================
-correction=mean(temperatureEffect(Tref, 3, 17, 26))
+correction=mean(temperatureEffect(TrefAtSea, 3, 17, 26))
 
 age=seq(0,10,.25)
 present = vonBertalaffyGrowth(age, 2, 60, 0.3900707 * correction)

exploration/scriptR/PatrickPrg.R

 # see BruchPrg for upload of data
 
-# ======================================================================
-#  temp en fonction de la latitude
-# ====================================================================
+basins = read.csv("/home/patrick/Documents/workspace/GR3D/data/input/reality/basins.csv", sep = ";")
+tempData = read.csv("/home/patrick/Documents/workspace/GR3D/data/input/reality/SeasonTempBVFacAtlant1801_2100_newCRU_RCP85.csv", sep = ";")
+
+riverName = "Garonne"
+riverName = "Rhine"
+riverName = "Meuse"
+
+surfaceBasin = basins[basins$nomBV == riverName, 'surface']
+
+sel = tempData$NOM == riverName & tempData$Year >= 1901 & tempData$Year <= 1910
+(TrefInRiver = as.data.frame(t(colMeans(tempData[sel, c("Winter", "Spring", "summer", "Autumn")]))))
+
+
+# =============================== temp en fonction de la latitude ====
 library(dplyr)
 
-extractTemp<-as.data.frame(tempData %>%
-                             filter(Year>1901 & Year <= 1910) %>%
-                             group_by(NOM,Lat)%>%
-                             summarize(meanSpring=mean(Spring)) 
-                           %>% arrange (Lat)
+extractTemp <- as.data.frame(tempData %>%
+                               filter(Year > 1991 & Year <= 2009) %>%
+                               group_by(NOM,Lat) %>%
+                               summarize(meanSpring = mean(Spring)) 
+                             %>% arrange(Lat)
 )
 
-plot(extractTemp$Lat, extractTemp$meanSpring, pch=20)
-abline(v=49.25) # vire
-abline(v=53.75) # elbe
+plot(extractTemp$Lat, extractTemp$meanSpring, pch = 20)
+abline(v = 45.25) # Garonne
+abline(v = 49.25) # vire
+abline(v = 52.25) # Rhine
+
+extractTemp %>% filter(NOM == riverName)
 
+abline(h = 11)
+abline(h = 12.5)
 
-abline(h=11)
-abline(h=12.5)
+levels(extractTemp$NOM)
 
 
-# ========================================================================
-# effect temperature sur la reproduction
-# ========================================================================
+# ==== effect temperature sur la reproduction ====================================
 
-thermalBevertonHolt = function (stock, temperature, juvenileTolerance) {
+
+temperatureEffect = function(tempWater, TminRep, ToptRep, TmaxRep){
+  #  if (tempWater<=TminRep | tempWater >= TmaxRep)
+  #    return(0)
+  #  else
+  response = (tempWater - TminRep) * (tempWater - TmaxRep) / ((tempWater - TminRep) * (tempWater - TmaxRep) - (tempWater - ToptRep)^2)
+  
+  response[tempWater <= TminRep | tempWater >= TmaxRep] = 0
+  return(response)
+  
+}
+
+thermalBevertonHolt = function(stock, temperature, juvenileTolerance, surface = 84810) {
   survOptRep = 0.0017
-  surface =84810
   delta_t = 0.33
   lambda = 4.1E-4
-  a=135000 # ~fecondité
+  a = 135000 # ~fecondité
   tempEffect = temperatureEffect(temperature, juvenileTolerance[1], juvenileTolerance[2], juvenileTolerance[3] )
   b =  -log(survOptRep * tempEffect)/delta_t
   c = lambda/surface
   
   eta = 2.4
-  theta =1.9
+  theta = 1.9
   #theta = 1.15
   S95 = eta * surface
   S50 = S95 / theta
   
-  AlleeEffect = 1 / (1+ exp(-log(19) * (S - S50) / (S95 - S50))) 
+  AlleeEffect = 1 / (1 + exp(-log(19) * (stock - S50) / (S95 - S50))) 
   
   alpha = (b * exp(-b * delta_t)) / (c * (1 - exp(-b * delta_t)))
   alpha[tempEffect == 0] = 0
-  beta = b /(a * c * (1-exp(-b * delta_t)))
+  beta = b / (a * c * (1 - exp(-b * delta_t)))
   
   recruit = alpha * stock * AlleeEffect / (beta + stock * AlleeEffect)
-  return (recruit)
+  return(recruit)
+}
+
+stockAtCrash = function(temperature, juvenileTolerance, surface = 84810) {
+  survOptRep = 0.0017
+  delta_t = 0.33
+  lambda = 4.1E-4
+  a = 135000 # ~fecondité
+  tempEffect = temperatureEffect(temperature, juvenileTolerance[1], juvenileTolerance[2], juvenileTolerance[3])
+  b =  -log(survOptRep * tempEffect) / delta_t
+  b[b == Inf] = 0
+  c = lambda / surface
+  
+  eta = 2.4
+  theta = 1.9
+  #theta = 1.15
+  S95 = eta * surface
+  S50 = S95 / theta
+  
+  beta = b / (a * c * (1 - exp(-b * delta_t)))
+  
+  res = S50 + (S95 - S50) * log(beta * log(19) / (S95 - S50)) / log(19)
+  res[is.nan(res)] = 0
+  return(res)
 }
 
-thermalBevertonHoltWithoutAllee = function (stock, temperature, juvenileTolerance) {
+
+thermalBevertonHoltWithoutAllee = function(stock, temperature, juvenileTolerance, surface = 84810)  {
   survOptRep = 0.0017
-  surface =84810
   delta_t = 0.33
   lambda = 4.1E-4
-  a=135000 # ~fecondité
+  a = 135000 # ~fecondité
   tempEffect = temperatureEffect(temperature, juvenileTolerance[1], juvenileTolerance[2], juvenileTolerance[3] )
   b =  -log(survOptRep * tempEffect)/delta_t
   c = lambda/surface
   
-  
   alpha = (b * exp(-b * delta_t)) / (c * (1 - exp(-b * delta_t)))
   alpha[tempEffect == 0] = 0
-  beta = b /(a * c * (1-exp(-b * delta_t)))
+  beta = b / (a * c * (1 - exp(-b * delta_t)))
   
   recruit = alpha * stock  / (beta + stock )
-  return (recruit)
+  return(recruit)
 }
 
 
 modifiedBevertonHolt = function(stock){
-  alpha = 6.4e6
+  alpha = 6.4e6 # for the Garonne River !!!
   beta = .172e6
-  d=19.2
+  d = 19.2
   
-  return  (alpha * (S^d) / ((beta^d) + (S^d)))
+  return(alpha * (stock^d) / ((beta^d) + (stock^d)))
+}
+
+survivingRecruit = function(stock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface){
+  survStock = stock * temperatureEffect(temperature, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3])
+  recruit = thermalBevertonHolt(survStock, temperature, juvenileTolerance, surface)
+  survRecruit =  recruit * exp(-sigmaZ)
+  return(survRecruit)
 }
 
-temperature =18
 
-juvenileTolerance = c(9.8, 20, 26)
+
 spawnerTolerance = c(10, 20, 23)
-Z=0.4
+juvenileTolerance = c(9.8, 20, 26)
 
-S=seq(0, 1e6, 100)
-plot(S, thermalBevertonHolt(S, Tref[2], juvenileTolerance), type ='l', ylab ='R', ylim = c(0, 8e6) )
-lines(S, thermalBevertonHolt(S, Tref[2]+2, juvenileTolerance), col='red')
-lines(S, thermalBevertonHolt(S, Tref[2]-2, juvenileTolerance), col='blue')
 
+spawnerTolerance = c(10.7, 17, 25.7) #  spawning temperature range Paumier et al 2019
+juvenileTolerance = c(10.8, 20.8, 29.8) #tolerance temperature range Jatteau 2017
+
+Tmin = 8
+spawnerTolerance = c(Tmin, 20, 23) 
+juvenileTolerance = c(Tmin, 16, 23) 
+
+sigmaZ = 0.4 * 5 # Z * lifespan of female
+propFemale = 0.5
+
+# ================================ stock-recruitment relationship at 1900-1910 SPRING temperature ====
+
+ymax = 1e7
+S = seq(0, 1e6, 100)
+# with spawner mortality according to temperature
+survS = S * temperatureEffect(TrefInRiver$Spring, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3])
+plot(S, thermalBevertonHolt(survS, TrefInRiver$Spring, juvenileTolerance, surface = surfaceBasin), type = 'l',
+     xlab = 'Stock (female)', ylab = 'R (female and male)', ylim = c(0, ymax) )
+lines(S, thermalBevertonHolt(survS, TrefInRiver$Spring + 2, juvenileTolerance, surface = surfaceBasin), col = 'red')
+lines(S, thermalBevertonHolt(survS, TrefInRiver$Spring - 2, juvenileTolerance, surface = surfaceBasin), col = 'blue')
+# replacement line with prop of female in recruit
+lines(c(0, ymax * exp(-sigmaZ) * propFemale), c(0, ymax), col = 'green', lty = 2)
+
+
+# =============================== comparison with Rougier et al stock-recruitment relationship for the Garonne River ======
+
+tempe = TrefInRiver$Spring
+tempe = as.numeric(tempData %>%
+  filter(Year > 1991 & Year <= 2009 & NOM == riverName) %>%
+  summarize(meanSpring = mean(Spring), meanSummer = mean(summer)) )
+tempe = juvenileTolerance[2]
+tempe = 20.2 # Rougier 2014
+
+Stot = seq(0, 2e6, 1e4)
+ymax = 1e7
+
+survS = Stot * temperatureEffect(tempe, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3])
+plot(Stot, thermalBevertonHolt(stock = survS * propFemale, temperature = tempe, juvenileTolerance = juvenileTolerance, 
+                               surface = surfaceBasin),
+     type = 'l', xlab = 'Total stock (2 x female)', ylab = 'R (female and male)', ylim = c(0, ymax) )
+# lines(Stot, thermalBevertonHolt(stock = survS, temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin), lty = 3)
+lines(Stot, thermalBevertonHoltWithoutAllee(Stot * propFemale, tempe, juvenileTolerance,  surface = surfaceBasin), col = 'blue')
+## lines(Stot, thermalBevertonHolt(Stot, tempe, juvenileTolerance), lty=2 )
+lines(Stot, modifiedBevertonHolt(Stot), col = 'green') # for the Garonne River
+# replacement line
+lines(c(0, ymax * exp(-sigmaZ)), c(0, ymax), col = 'green', lty = 2)
+
+# ===================================== stock at equilibrium =============================================================
+stocksAtEquilibrium = function(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface, propFemale = 0.5) {
+
+   # !!!!! stock and recruit consider only females !!!!!!
+  objFct = function(stock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface){
+    survRecruit = survivingRecruit(stock, temperature, juvenileTolerance, spawnerTolerance, sigmaZ, surface) * propFemale 
+    return((survRecruit - stock)^2)
+  }
+  
+  stockatequilibrium = 0
+  stockattrap = 0
+  
+  # stock at crash before the spawner mortality
+  tempEffSpawner =  temperatureEffect(temperature, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3])
+  stockatcrash = ifelse(tempEffSpawner > 0, stockAtCrash(temperature, juvenileTolerance, surface) / tempEffSpawner, 0)
+  
+  
+  if (stockatcrash > 0) {
+    recruitatcrash = thermalBevertonHolt(stockatcrash * tempEffSpawner,
+                                         temperature, juvenileTolerance, surface) * propFemale 
+    sigmaZatcrash = -log(stockatcrash / recruitatcrash)
+    
+    if (sigmaZ < sigmaZatcrash) {
+      stockattrap =  optimise(objFct, c(0,  stockatcrash), temperature = temperature, juvenileTolerance = juvenileTolerance, spawnerTolerance = spawnerTolerance, 
+                              sigmaZ = sigmaZ, surface = surface)$minimum
+      
+      stockatequilibrium = optimise(objFct, (c(1, 100) * stockatcrash), temperature = temperature, juvenileTolerance = juvenileTolerance, spawnerTolerance = spawnerTolerance, 
+                                      sigmaZ = sigmaZ, surface = surface)$minimum
+    }
+  }
+  
+  return(c(stockatequilibrium, stockattrap, stockatcrash))
+}
 
-Stot=seq(0, 0.3e6, 100)*2
-plot(Stot, thermalBevertonHolt(Stot/2, juvenileTolerance[2], juvenileTolerance), type ='l', ylab ='R', ylim = c(0, 8e6) )
-lines(Stot, thermalBevertonHoltWithoutAllee(Stot/2, juvenileTolerance[2], juvenileTolerance), col ='blue')
-## lines(Stot, thermalBevertonHolt(Stot, juvenileTolerance[2], juvenileTolerance), lty=2 )
-lines(Stot, modifiedBevertonHolt(Stot), col = 'green')
+temperature = seq(5,30,.1)
 
+stocks = as.data.frame(t(sapply(temperature, function(temp) (stocksAtEquilibrium(temp, juvenileTolerance, spawnerTolerance, sigmaZ, surface = surfaceBasin)))))
+names(stocks) = c('atEquilibrium', 'atTrap', 'atCrash')
+stocks$temperature = temperature
 
+sel = TRUE
+plot(stocks$temperature[sel], stocks$atEquilibrium[sel], type = 'l', xlab = 'T (°C)', ylab = 'Stock (# female)' , main = riverName)
+lines(stocks$temperature[sel], stocks$atTrap[sel], col = 'red')
+legend('topleft', legend = c('at equilibrium', 'at trap'), col = c('black', 'red'), lty = 1)
+max(stocks$temperature[stocks$atEquilibrium == 0 & stocks$temperature < 20])
 
-thermBH_without = sapply(temperature, function(temp) (thermalBevertonHolt(stock = 4e5, temp)))
-plot(temperature, thermBH_without, type='line,', ylab ='recruit')
+# verif
+tempe = 14.4
+stocks[stocks$temperature == tempe, ]
 
-thermBH_with = sapply(temperature, function(temp) (thermalBevertonHolt(stock = 4e5 *  
-                                                                         temperatureEffect(temp, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]), temp, juvenileTolerance)))
-lines(temperature, thermBH_with, col='red')
 
 
-recruitProduction = function(temperature, juvenileTolerance, spawnerTolerance, Z, stock) {
-  survivingStock = stock * temperatureEffect(temperature, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3])
-  recruit = sapply(1:length(temperature), function(i) (thermalBevertonHolt(survivingStock[i], temperature[i], juvenileTolerance)))
-  survivingRecruit =  recruit * exp(-Z*5)
-  return(survivingRecruit)
-}
+S = seq(0, 5e5, 1e3) # only females
 
+recruit = thermalBevertonHolt(S * temperatureEffect(tempe, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]),
+                              temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin) * propFemale
 
-fct = function(temperature, juvenileTolerance, spawnerTolerance, Z, stock){
-  return( (recruitProduction(temperature, juvenileTolerance, spawnerTolerance, Z, stock) -stock)^2)
-}
+plot(S, recruit, type = 'l', xlab = 'Stock (# female)', ylab = 'Recruit (# of female)', main = paste(riverName, surfaceBasin))
+lines(c(0, max(recruit) * exp(-sigmaZ)), c(0, max(recruit)), col = 'green', lty = 2)
+sAtEqui = stocks$atEquilibrium[stocks$temperature == tempe]
+sAtCrash = stocks$atCrash[stocks$temperature == tempe]
+sAtTrap = stocks$atTrap[stocks$temperature == tempe]
+points(sAtEqui, thermalBevertonHolt(sAtEqui * temperatureEffect(tempe, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]),
+                                    temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin) * propFemale,
+       pch = 20) 
+points(sAtTrap, thermalBevertonHolt(sAtTrap * temperatureEffect(tempe, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]),
+                                    temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin) * propFemale, 
+       col = 'red', pch = 20 )
+points(sAtCrash, thermalBevertonHolt(sAtCrash * temperatureEffect(tempe, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]),
+                                    temperature = tempe, juvenileTolerance = juvenileTolerance, surface = surfaceBasin) * propFemale)
 
 
-stock =4e5
-temperature =seq(0,30,.1)
-plot(temperature, recruitProduction(temperature, juvenileTolerance, spawnerTolerance, Z, stock), type='l')
-plot(temperature, fct(temperature, juvenileTolerance, spawnerTolerance, Z, stock), type='l')
+# ==== recruitment according to temperature ===========================
+# but need to fix a reference stock !
+refStock = 1.5e6 / 2 # female !!!
+temperature = seq(5,30,.1)
 
+thermBH_with = sapply(temperature, function(temp) (thermalBevertonHolt(stock = refStock *  
+                                                                         temperatureEffect(temp, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3]), temp, juvenileTolerance)))
+plot(temperature, thermBH_with, type = "l", ylab = 'recruit')
+
+# verify impact of Allee effect 
+thermBH_without = sapply(temperature, function(temp) (thermalBevertonHolt(stock = refStock, temp, juvenileTolerance)))
+lines(temperature, thermBH_without, col = 'red')
 
-optimise(fct, interval = c(5, 20),  stock=stock, Z=Z, juvenileTolerance = c(9, 18, 26), spawnerTolerance = c(9, 18, 26) )$minimum
 
-plot(temperature, thermBH_with * exp(-Z*5), ylab = 'surviving spawner' , type = 'line')
+# ===== with a fixed reference stock
+recruitProduction = function(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, refStock) {
+  survivingStock = refStock * temperatureEffect(temperature, spawnerTolerance[1], spawnerTolerance[2], spawnerTolerance[3])
+  recruit = sapply(1:length(temperature), function(i) (thermalBevertonHolt(survivingStock[i], temperature[i], juvenileTolerance)))
+  survivingRecruit =  recruit * exp(-sigmaZ) # surviving spawners
+  return(survivingRecruit)
+}
 
+refStock = 2e5 
+temperature = seq(5,30,.1)
+plot(temperature, recruitProduction(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, refStock = refStock), type = 'l',
+     ylab = "survival spawner")
+abline(h = refStock)
 
-abline(h=4e5)
-#lines(temperature, thermBH_with * exp(-(Z/2)*5), col='green')
-cbind(temperature, thermBH_with *exp(-Z*5) -4e5)
+fct = function(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, refStock){
+   return( (recruitProduction(temperature, juvenileTolerance, spawnerTolerance, sigmaZ, refStock) - refStock)^2)
+ }
+ optimise(fct, interval = c(5, 20),  refStock=refStock, sigmaZ=sigmaZ, juvenileTolerance , spawnerTolerance )$minimum
 
 
+# ===================
+# OLD
+# ==================
 # temperature effect on spawner survival 
 plot(temperature, temperatureEffect(temperature, 10, 20, 23), type='l', type = 'line')