diff --git a/Trajectories/cost_function.py b/Trajectories/cost_function.py
index b6681c927f75606fad29e43a6f30beba2b4649b1..0c09635cb89170e81975a4e3b939e1f9abc7685b 100755
--- a/Trajectories/cost_function.py
+++ b/Trajectories/cost_function.py
@@ -174,41 +174,44 @@ print("Iterating on vertices...")
 for n in g.vertices():
     if n.in_degree() != 0:
         # Define properties for beginning of the segment:
+        ip_beg = min(n_days_avg, len(g.vp.inst_eddies[n]))
         properties = node_to_prop(
-            g.vp.inst_eddies[n][:n_days_avg],
+            g.vp.inst_eddies[n][:ip_beg],
             g.gp.e_overestim,
             SHPC,
             args.orientation,
         )
         first_av_rad[n], first_av_ros[n] = calculate_radii_rossby(properties)
         pos_first[n] = properties[0]["pos"]  # in rad
+    else:
+        ip_beg = 0
 
     if n.out_degree() != 0:
         # Define properties for end of the segment:
 
         len_seg = len(g.vp.inst_eddies[n])
 
-        if n.in_degree() == 0 or len_seg > n_days_avg:
+        if ip_beg < len_seg:
             # We have to read more from the shapefiles and redefine
             # properties.
+            ip_end = max(len_seg - n_days_avg, 0)
 
-            if n.in_degree() == 0 or len_seg >= 2 * n_days_avg:
+            if ip_beg <= ip_end:
                 # We cannot use part of properties from the beginning
                 # of the segment.
                 properties = node_to_prop(
-                    g.vp.inst_eddies[n][-n_days_avg:],
+                    g.vp.inst_eddies[n][ip_end:],
                     g.gp.e_overestim,
                     SHPC,
                     args.orientation,
                 )
             else:
-                # assertion: n.in_degree() != 0 and n_days_avg <
-                # len_seg < 2 * n_days_avg
+                # assertion: ip_end < ip_beg < len_seg
 
                 # We can use part of the properties from the beginning
                 # of the segment.
-                properties = properties[len_seg - n_days_avg:] + node_to_prop(
-                    g.vp.inst_eddies[n][n_days_avg:],
+                properties = properties[ip_end:] + node_to_prop(
+                    g.vp.inst_eddies[n][ip_beg:],
                     g.gp.e_overestim,
                     SHPC,
                     args.orientation,