Store position in rad

We avoid multiple conversions to rad of the same value, in the loop on
edges, since a segment can be head or tail of several edges.

Trajectories/cost_function.py

@@ -43,7 +43,7 @@ def calculate_radii_rossby(properties):
     n_valid_Rossby = 0
 
     for prop in properties:
-        f = 2 * Omega * math.sin(math.radians(prop["pos"][1])) # in s-1
+        f = 2 * Omega * math.sin(prop["pos"][1]) # in s-1
         radius = prop["radius"] * 1000 # in m
 
         if abs(prop["speed"]) < 100:
@@ -75,7 +75,7 @@ def node_to_prop(node_list, e_overestim, SHPC, orientation):
         ishape = SHPC.comp_ishape(date_index, eddy_index, i_slice, orientation)
         shapeRec = SHPC.get_reader(i_slice, orientation, "extremum")\
                        .shapeRecord(ishape)
-        prop = {"pos": shapeRec.shape.points[0], "speed": shapeRec.record.speed}
+        prop = {"pos": [math.radians(x) for x in shapeRec.shape.points[0]], "speed": shapeRec.record.speed}
         prop["radius"] \
             = SHPC.get_reader(i_slice, orientation, "max_speed_contour")\
                   .record(ishape).r_eq_area
@@ -156,7 +156,7 @@ for n in g.vertices():
         properties = node_to_prop(g.vp.inst_eddies[n][:n_days_avg],
                                   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 degrees
+        pos_first[n] = properties[0]["pos"] # in rad
 
     if n.out_degree() != 0:
         # Define properties for end of the segment:
@@ -192,7 +192,7 @@ for n in g.vertices():
             last_av_rad[n] = first_av_rad[n]
             last_av_ros[n] = first_av_ros[n]
 
-        pos_last[n] = properties[- 1]["pos"] # in degrees
+        pos_last[n] = properties[- 1]["pos"] # in rad
 
 t1 = time.perf_counter()
 timings.write(f"iterating on vertices: {t1 - t0:.0f} s\n")
@@ -203,18 +203,17 @@ for edge in g.edges():
     source_node = edge.source()
     target_node = edge.target()
 
-    latitude = math.radians((pos_last[source_node][1] +
-                             pos_first[target_node][1]) / 2)
-    # (latitude needed for conversion of degrees to kilometers)
+    latitude = (pos_last[source_node][1] +
+                             pos_first[target_node][1]) / 2
 
     # Because of the wrapping issue (360° wrapping incorrectly to 0°),
     # we check for that here:
     lon_diff = abs(pos_last[source_node][0] - pos_first[target_node][0])
-    if lon_diff > 300: lon_diff = 360 - lon_diff
+    if lon_diff > math.radians(300): lon_diff = 2 * math.pi - lon_diff
 
-    Delta_Cent = math.sqrt((lon_diff * 111.32 * math.cos(latitude))**2
+    Delta_Cent = math.sqrt((lon_diff * 6378.166175 * math.cos(latitude))**2
                            + ((pos_last[source_node][1]
-                               - pos_first[target_node][1]) * 110.574)**2)
+                               - pos_first[target_node][1]) * 6335.423523)**2)
 
     # Rossbies:
     if (first_av_ros[target_node] and last_av_ros[source_node]):