diff --git a/cost_function.py b/cost_function.py
index 26e616fb780ec871a1b3f90de108224eabd7529e..925947d693f94dd06437b52d7f4d073c7792bf24 100755
--- a/cost_function.py
+++ b/cost_function.py
@@ -28,9 +28,11 @@ import argparse
 
 Omega = 2 * math.pi / 86164.
 
-def calculate_radii_rossby(list_eddies, e_overestim, handlers, array_d_init):
-    """Compute average on list_eddies of Rossby number and radius of
-    maximum speed contour.
+def calculate_radii_rossby(properties):
+    """Compute average on some instantaneous eddies of Rossby number and
+    radius of maximum speed contour. The required properties for each
+    eddy are position, radius and speed. "properties" is a list of
+    dictionaries. Each dictionary in the list contains the three properties.
 
     """
 
@@ -38,26 +40,17 @@ def calculate_radii_rossby(list_eddies, e_overestim, handlers, array_d_init):
     avg_Rossby = 0
     n_valid_Rossby = 0
 
-    for n in list_eddies:
-        date_index, eddy_index = report_graph.node_to_date_eddy(n, e_overestim)
-        i_SHPC = get_SHPC(array_d_init, date_index)
-        ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
-                                         eddy_index)
+    for prop in properties:
+        f = 2 * Omega * math.sin(math.radians(prop["pos"][1])) # in s-1
+        radius = prop["radius"] * 1000 # in m
 
-        # Now that we have the location in the shapefiles, we need to
-        # get the radius and the rossby number:
-        shapeRec = handlers[i_SHPC]["readers"]["extremum"].shapeRecord(ishape)
-        f = 2 * Omega * math.sin(math.radians(shapeRec.shape.points[0][1]))
-        radius = handlers[i_SHPC]["readers"]["max_speed_contour"]\
-            .record(ishape).r_eq_area * 1000 # in m
-
-        if abs(shapeRec.record.speed) < 100:
-            avg_Rossby += shapeRec.record.speed / (f * radius)
+        if abs(prop["speed"]) < 100:
+            avg_Rossby += prop["speed"] / (f * radius)
             n_valid_Rossby += 1
 
         avg_rad += radius # in m
 
-    avg_rad /= len(list_eddies)
+    avg_rad /= len(properties)
     if n_valid_Rossby != 0: avg_Rossby /= n_valid_Rossby
     return avg_rad, avg_Rossby
 
@@ -66,6 +59,30 @@ def get_SHPC(array_d_ini, date_index):
     assert i_SHPC >= 0
     return i_SHPC
 
+def node_to_prop(node_list, e_overestim, array_d_init, handlers):
+    """node_list is a list of node identification numbers for
+    instantaneous eddies. This function returns some properties of the
+    eddies, read from shapefiles: position of extremum, radius of
+    outermost contour or maximum speed contour, and speed. The three
+    properties are in a dictionary, for each eddy.
+
+    """
+
+    properties = []
+
+    for n in node_list:
+        date_index, eddy_index = report_graph.node_to_date_eddy(n, e_overestim)
+        i_SHPC = get_SHPC(array_d_init, date_index)
+        ishape = util_eddies.comp_ishape(handlers[i_SHPC], date_index,
+                                         eddy_index)
+        shapeRec = handlers[i_SHPC]["readers"]["extremum"].shapeRecord(ishape)
+        prop = {"pos": shapeRec.shape.points[0], "speed": shapeRec.record.speed}
+        prop["radius"] = handlers[i_SHPC]["readers"]["max_speed_contour"]\
+            .record(ishape).r_eq_area
+        properties.append(prop)
+
+    return properties
+
 t0 = time.perf_counter()
 timings = open("timings.txt", "w")
 parser = argparse.ArgumentParser()
@@ -156,20 +173,23 @@ for n in g.vertices():
             # to average
 
             # First 7 days calculation
+            properties = node_to_prop(segment[:n_days_avg], e_overestim,
+                                      array_d_init, handlers)
             g.vp.first_av_rad[n], g.vp.first_av_ros[n] \
-                = calculate_radii_rossby(segment[:n_days_avg], e_overestim,
-                                         handlers, array_d_init)
+                = calculate_radii_rossby(properties)
 
             # Last 7 days calculation:
+            properties = node_to_prop(segment[- n_days_avg:], e_overestim,
+                                      array_d_init, handlers)
             g.vp.last_av_rad[n], g.vp.last_av_ros[n] \
-                = calculate_radii_rossby(segment[- n_days_avg:], e_overestim,
-                                         handlers, array_d_init)
+                = calculate_radii_rossby(properties)
         else:
             # The number of eddies in a segment is lower than the number
             # of days over which to average. The values will be the same
             # except for the positions.
-            avg_rad, avg_Rossby = calculate_radii_rossby(segment, e_overestim,
-                                                         handlers, array_d_init)
+            properties = node_to_prop(segment, e_overestim, array_d_init,
+                                      handlers)
+            avg_rad, avg_Rossby = calculate_radii_rossby(properties)
 
             # Average and assign the rossbies:
             g.vp.first_av_ros[n] = avg_Rossby