--- STD/helpers.hh	2020-04-30 09:19:51.575537041 +0200
+++ MOD/helpers.hh	2020-05-02 14:36:21.035245866 +0200
@@ -1,5 +1,8 @@
+// helpers.hh edited by Lorenzo Sala
+
 #include "nuc_eos.hh"
 #include <cstdlib>
+#include <iostream>
 
 namespace nuc_eos {
 
@@ -368,7 +371,7 @@
   const double dltp = log(1.2);
   const double dltm = log(0.8);
 
-  const double leps0_prec = leps0*prec;
+  double leps0_prec = fabs(leps0*prec);
 
   // temporary local vars
   double lt, lt1, lt2;
@@ -380,8 +383,31 @@
   double delx,dely,delz;
   int idx[8];
   
+  // LSMOD (Modification made by Lorenzo Sala)
+  // LSMOD: The following lines calculate eps in 
+  //        f2a = eps(rho,Tmin, Ye) and f1a = eps(rho,Tmax,Ye)
+  get_interp_spots(lr,ltmax,ye,&delx,&dely,&delz,idx);
+  nuc_eos_C_linterp_one(idx,delx,dely,delz,&f1a,iv);
+  get_interp_spots(lr,ltmin,ye,&delx,&dely,&delz,idx);
+  nuc_eos_C_linterp_one(idx,delx,dely,delz,&f2a,iv);
 
   // prepare
+  // check if your energy is actually tabulated at this rho and ye.
+  // f2a is the energy evaluated at ltmin, so it is the minimum energy tabulated 
+  // at this rho ad ye.
+  // If leps0 <= f2a, then ltout is likely to be the minimum temperature tabulated.
+  if(leps0 <= f2a) { // + 1.0E-6
+    *ltout = ltmin;
+    return;
+  }
+
+  /* // If leps0 >= f1a, then ltout is likely to be the maximum temperature tabulated.
+  if(leps0 >= f1a) { // + 1.0E-6
+    *ltout = ltmax;
+    return;
+  } */
+
+  // otherwise, proceed finding extrema for applying bisection method.
   lt = lt0;
   lt1 = MIN(lt0 + dlt0p,ltmax);
   lt2 = MAX(lt0 + dlt0m,ltmin);
@@ -416,6 +442,10 @@
 
     bcount++;
     if(CCTK_BUILTIN_EXPECT(bcount >= maxbcount, false)) {
+#if DEBUG
+      fprintf(stderr,"bcount out of range it %d, lr: %15.6E, lt1: %15.6E, lt2: %15.6E, f1a: %18.11E, f2a: %18.11E leps0: %18.11E, ye: %15.6E\n",
+	      bcount,lr,lt1,lt2,f1a,f2a,leps0,ye);
+#endif
       *keyerrt = 667;
       return;
     }
@@ -510,7 +540,12 @@
   const int irho = MIN(MAX(1 + (int)(( lr - logrho[0] - 1.0e-12) * drhoi),1),nrho-1); 
   const int iye = MIN(MAX(1 + (int)(( ye - yes[0] - 1.0e-12) * dyei),1),nye-1); 
 
-  while(it < itmax) {
+  /* ******* if temp low for high density, switch directly to bisection. 
+             Verifying Newton-Raphson result evaluating the derivative. 
+             The variable shouldgotobisection will be modified accordingly
+             to the value of derivative of eps wrt temp ******* */
+  bool shouldgotobisection = false; // LSMOD
+  while(it < itmax && shouldgotobisection == false) {
     it++;
 
     // step 2: check if the two bounding values of the temperature
@@ -572,6 +607,13 @@
     dlepsdlti = (logtemp[itemp]-logtemp[itemp-1])/(epst2-epst1);
     ldt = -(leps - leps0) * dlepsdlti * fac;
 
+    //LSMOD: too large a dlt means that the energy dependence on the temperature
+    //       is weak ==> We'd better try bisection.
+    //       Factor 1/12.0 come from tests by LSMOD
+    //       This is done in order to limit the "velocity" of T variation
+    //       given by Newton-Raphson.    
+    if(ldt > (ltmax-ltmin) / 12.0 ) shouldgotobisection = true;
+
     ltn = MIN(MAX(lt + ldt,ltmin),ltmax);
     lt = ltn;
 
@@ -605,9 +647,9 @@
   bisection(lr,lt0,ye,leps0,prec,ltout,1,keyerrt);
 #if DEBUG
   if(*keyerrt==667) {
-    fprintf(stderr,"This is worse. Bisection failed!\n");
-    abort();
-  }      
+     fprintf(stderr,"This is worse. Bisection failed!\n");
+     abort();
+  }
 #endif