diff --git a/src/grib_iterator_class_gaussian.c b/src/grib_iterator_class_gaussian.c
index 7915fd8d2..56afa1661 100644
--- a/src/grib_iterator_class_gaussian.c
+++ b/src/grib_iterator_class_gaussian.c
@@ -90,15 +90,14 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
   grib_iterator_gaussian* self = (grib_iterator_gaussian*)i;
 
   double *lats;
-  double laf;
-  double lal;
-  long trunc;
+    double laf; /* latitude of first point in degrees */
+    double lal; /* latitude of last point in degrees */
+    long trunc; /* number of parallels between a pole and the equator */
   long lai;
   long jScansPositively=0;
   int size;
   double start;
   unsigned long istart=0;
-
   int ret = GRIB_SUCCESS;
 
   const char* latofirst   = grib_arguments_get_name(h,args,self->carg++);
@@ -106,7 +105,6 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
   const char* numtrunc    = grib_arguments_get_name(h,args,self->carg++);
   const char* s_jScansPositively    = grib_arguments_get_name(h,args,self->carg++);
 
-
   if((ret = grib_get_double_internal(h,latofirst,   &laf))) return ret;
   if((ret = grib_get_double_internal(h,latoflast,   &lal))) return ret;
   if((ret = grib_get_long_internal(h,numtrunc,&trunc))) return ret;
@@ -128,14 +126,13 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
 /*
   for(loi=(trunc*2)-1;loi>=0;loi--)
     if(fabs(lats[loi] - lal) < glatPrecision) break;
-
-
   for(j=(trunc*2)-1;j>0;j--) {
     if(fabs(lats[j] - laf) < glatPrecision) break;
   }
 */
 
   binary_search(lats,size-1,start,&istart);
+    Assert(istart < size);
 
   if (jScansPositively) {
     for(lai=0;lai<self->nam;lai++) {
@@ -152,21 +149,25 @@ static int init(grib_iterator* i,grib_handle* h,grib_arguments *args){
   grib_context_free(h->context,lats);
 
   return ret;
-
 }
 
 static void binary_search(double xx[], const unsigned long n, double x, unsigned long *j)
 {
-  /*These routine works only on descending ordered arrays*/
+    /*This routine works only on descending ordered arrays*/
+#define EPSILON 1e-3
+
   unsigned long ju,jm,jl;
   jl=0;
   ju=n;
   while (ju-jl > 1) {
     jm=(ju+jl) >> 1;
-    if (x <= xx[jm]) jl=jm;
+        if (fabs(x-xx[jm]) < EPSILON) {
+            /* found something close enough. We're done */
+            *j=jm;
+            return;
+        }
+        if (x < xx[jm]) jl=jm;
     else ju=jm;
   }
   *j=jl;
 }
-
-