diff --git a/Trajectories/cost_function.py b/Trajectories/cost_function.py
index b9bfa1abfc701a82ce1c04a7e134c470c359034d..b795e0dd487e39bee385274345b6cbb54094bd5e 100755
--- a/Trajectories/cost_function.py
+++ b/Trajectories/cost_function.py
@@ -103,12 +103,78 @@ def node_to_prop(node_list, e_overestim, SHPC, orientation):
     return properties
 
 
-def search_beg(inst_eddies, max_delta):
-    return min(max_delta + 1, len(inst_eddies))
+def search_beg(inst_eddies, max_delta, avg_fix, e_overestim):
+    """Return an index ip in inst_eddies. inst_eddies is a list of
+    inter-date indices of instantaneous eddies. len(inst_eddies) must
+    be >= 1. max_delta must be >= 0. If avg_fix then inst_eddies[:ip]
+    is a fixed number of instantaneous eddies, if the list is long
+    enought. If not avg_fix then ip = bisect.bisect_right(inst_eddies,
+    date(inst_eddies[0]) + max_delta, key = date), although the
+    computation is not implemented that way.
 
+    """
+    ip = min(max_delta + 1, len(inst_eddies))
+
+    if not avg_fix and ip >= 2:
+        d_max = (
+            util_eddies.node_to_date_eddy(
+                inst_eddies[0], e_overestim, only_date=True
+            )
+            + max_delta
+        )
+        # {date(elem) > d_max for elem in inst_eddies[ip:]}
+
+        while (
+            ip >= 2
+            and util_eddies.node_to_date_eddy(
+                inst_eddies[ip - 1], e_overestim, only_date=True
+            )
+            > d_max
+        ):
+            ip -= 1
+
+        # {date(elem) <= d_max for elem in inst_eddies[:ip] and
+        # date(elem) > d_max for elem in inst_eddies[ip:]}
+
+    # {1 <= ip <= min(max_delta + 1, len(inst_eddies))}
+    return ip
+
+
+def search_end(inst_eddies, max_delta, avg_fix, e_overestim):
+    """Return an index ip in inst_eddies. inst_eddies is a list of
+    inter-date indices of instantaneous eddies. len(inst_eddies) must
+    be >= 1. max_delta must be >= 0. If avg_fix then inst_eddies[ip:]
+    is a fixed number of instantaneous eddies, if the list is long
+    enought. If not avg_fix then ip = bisect.bisect_left(inst_eddies,
+    date(inst_eddies[-1]) - max_delta, key = date), although the
+    computation is not implemented that way.
+
+    """
+    ip = max(len(inst_eddies) - max_delta - 1, 0)
+
+    if not avg_fix and ip <= len(inst_eddies) - 2:
+        d_min = (
+            util_eddies.node_to_date_eddy(
+                inst_eddies[-1], e_overestim, only_date=True
+            )
+            - max_delta
+        )
+        # {date(elem) < d_min for elem in inst_eddies[:ip]}
+
+        while (
+            ip <= len(inst_eddies) - 2
+            and util_eddies.node_to_date_eddy(
+                inst_eddies[ip], e_overestim, only_date=True
+            )
+            < d_min
+        ):
+            ip += 1
 
-def search_end(inst_eddies, max_delta):
-    return max(len(inst_eddies) - max_delta - 1, 0)
+        # {date(elem) < d_min for elem in inst_eddies[:ip] and
+        # date(elem) >= d_min for elem in inst_eddies[ip:]}
+
+    # {max(len(inst_eddies) - max_delta - 1, 0) <= ip <= len(inst_eddies) - 1}
+    return ip
 
 
 t0 = time.perf_counter()
@@ -129,6 +195,11 @@ parser.add_argument(
 parser.add_argument(
     "--debug", help="save properties to output file", action="store_true"
 )
+parser.add_argument(
+    "--avg_fix",
+    help="average over a fixed number of instantaneous eddies",
+    action="store_true",
+)
 args = parser.parse_args()
 
 # Set some values needed for the cost function:
@@ -185,7 +256,9 @@ print("Iterating on vertices...")
 for n in g.vertices():
     if n.in_degree() != 0:
         # Define properties for beginning of the segment:
-        ip_beg = search_beg(g.vp.inst_eddies[n], max_delta)
+        ip_beg = search_beg(
+            g.vp.inst_eddies[n], max_delta, args.avg_fix, g.gp.e_overestim
+        )
         properties = node_to_prop(
             g.vp.inst_eddies[n][:ip_beg],
             g.gp.e_overestim,
@@ -203,7 +276,9 @@ for n in g.vertices():
         if ip_beg < len(g.vp.inst_eddies[n]):
             # We have to read more from the shapefiles and redefine
             # properties.
-            ip_end = search_end(g.vp.inst_eddies[n], max_delta)
+            ip_end = search_end(
+                g.vp.inst_eddies[n], max_delta, args.avg_fix, g.gp.e_overestim
+            )
 
             if ip_beg <= ip_end:
                 # We cannot use part of properties from the beginning