diff --git a/LoopStructural/modelling/features/builders/_fault_builder.py b/LoopStructural/modelling/features/builders/_fault_builder.py
index 8aa5836e4..eae038904 100644
--- a/LoopStructural/modelling/features/builders/_fault_builder.py
+++ b/LoopStructural/modelling/features/builders/_fault_builder.py
@@ -29,7 +29,9 @@ def __init__(
             the maximum area around the model domain that a fault is modelled. For high displacement faults this
             may need to be large, smaller values will be result in fewer degrees of freedom = quicker interpolation
         """
-        from LoopStructural.modelling.features.fault import FaultSegment  # defer import until needed
+        from LoopStructural.modelling.features.fault import (
+            FaultSegment,
+        )  # defer import until needed
 
         StructuralFrameBuilder.__init__(
             self, interpolator, interpolators, frame=FaultSegment, **kwargs
diff --git a/LoopStructural/modelling/intrusions/intrusion_frame_builder.py b/LoopStructural/modelling/intrusions/intrusion_frame_builder.py
index 637bc8fb8..04c4c2ef3 100644
--- a/LoopStructural/modelling/intrusions/intrusion_frame_builder.py
+++ b/LoopStructural/modelling/intrusions/intrusion_frame_builder.py
@@ -7,6 +7,7 @@
 )
 from LoopStructural.utils import getLogger
 
+
 logger = getLogger(__name__)
 
 import numpy as np
@@ -154,10 +155,14 @@ def create_grid_for_indicator_fxs(self, spacing=None):
 
     def add_contact_anisotropies(self, series_list=None, **kwargs):
         """
-        Add to the intrusion network the anisotropies likely exploited by the intrusion (series-type geological features)
-        Given a list of series-type features, evaluates contact points on each series and
-        compute mean value and standard deviation. Different contacts of the same series are indentify using clustering algortihm.
-        Mean and std deviation values will be used to identify each contact thoughout the model using the indicator functions.
+        Add to the intrusion network the anisotropies likely
+        exploited by the intrusion (series-type geological features)
+        Given a list of series-type features, evaluates contact points
+        on each series and compute mean value and standard deviation.
+        Different contacts of the same series are indentify using
+        clustering algortihm. Mean and std deviation values will
+        be used to identify each contact thoughout the model using
+        the indicator functions.
 
         Parameters
         ----------
@@ -165,6 +170,9 @@ def add_contact_anisotropies(self, series_list=None, **kwargs):
 
         Returns
         -------
+
+        Note
+        -----
         assigns to self.anisotropies_series_parameters a list like (for each strtaigraphic contact) =
         [series_name, mean of scalar field vals, standar dev. of scalar field val]
 
@@ -172,57 +180,53 @@ def add_contact_anisotropies(self, series_list=None, **kwargs):
         if self.intrusion_network_type == "shortest path":
             n_clusters = self.number_of_contacts
 
-        if series_list == None:
-            self.anisotropies_series_list = None
-
-        else:
-            self.anisotropies_series_list = series_list
-            series_parameters = {}
-            for i in range(len(series_list)):
-
-                series_list[
-                    i
-                ].faults_enabled = True  # data points evaluated in faulted hostrock
-                data_temp = self.intrusion_network_data[
-                    self.intrusion_network_data["intrusion_anisotropy"]
-                    == series_list[i].name
-                ].copy()
-                data_array_temp = data_temp.loc[:, ["X", "Y", "Z"]].to_numpy()
-                series_i_vals = series_list[i].evaluate_value(data_array_temp)
-                series_array = np.zeros([len(data_array_temp), 4])
-                series_array[:, :3] = data_array_temp
-                series_array[:, 3] = series_i_vals
-
-                n_contacts = n_clusters[i]
-
-                # use scalar field values to find different contacts
-                series_i_vals_mod = series_i_vals.reshape(len(series_i_vals), 1)
-                contact_clustering = KMeans(n_clusters=n_contacts, random_state=0).fit(
-                    series_i_vals_mod
-                )
+        self.anisotropies_series_list = series_list
+        series_parameters = {}
+        for i, series in enumerate(series_list):
+
+            series.faults_enabled = True  # data points evaluated in faulted hostrock
+            data_temp = self.intrusion_network_data[
+                self.intrusion_network_data["intrusion_anisotropy"] == series.name
+            ].copy()
+            data_array_temp = data_temp.loc[:, ["X", "Y", "Z"]].to_numpy()
+            series_i_vals = series.evaluate_value(data_array_temp)
+            series_array = np.zeros((len(data_array_temp), 4))
+            series_array[:, :3] = data_array_temp
+            series_array[:, 3] = series_i_vals
+
+            n_contacts = n_clusters[i]
+
+            # use scalar field values to find different contacts
+            series_i_vals_mod = series_i_vals.reshape(len(series_i_vals), 1)
+            # TODO create global loopstructural random state variable
+            contact_clustering = KMeans(n_clusters=n_contacts, random_state=0).fit(
+                series_i_vals_mod
+            )
 
-                for j in range(n_contacts):
-                    z = np.ma.masked_not_equal(contact_clustering.labels_, j)
-                    y = np.ma.masked_array(series_i_vals, z.mask)
-                    series_ij_vals = np.ma.compressed(y)
-                    series_ij_mean = np.mean(series_ij_vals)
-                    series_ij_std = np.std(series_ij_vals)
-                    series_ij_name = series_list[i].name + "_" + str(series_ij_mean)
-
-                    series_parameters[series_ij_name] = [
-                        series_list[i],
-                        series_ij_mean,
-                        series_ij_std,
-                    ]
+            for j in range(n_contacts):
+                z = np.ma.masked_not_equal(contact_clustering.labels_, j)
+                y = np.ma.masked_array(series_i_vals, z.mask)
+                series_ij_vals = np.ma.compressed(y)
+                series_ij_mean = np.mean(series_ij_vals)
+                series_ij_std = np.std(series_ij_vals)
+                series_ij_name = f"{series.name}_{str(series_ij_mean)}"
+
+                series_parameters[series_ij_name] = [
+                    series,
+                    series_ij_mean,
+                    series_ij_std,
+                ]
 
             self.anisotropies_series_parameters = series_parameters
 
     def add_faults_anisotropies(self, fault_list=None):
         """
-        Add to the intrusion network the anisotropies likely exploited by the intrusion (fault-type geological features)
-        Given a list of fault features, evaluates the contact points on each fault and
-        compute mean value and standard deviation.
-        These values will be used to identify each fault with the indicator function.
+        Add to the intrusion network the anisotropies likely
+        exploited by the intrusion (fault-type geological features)
+        Given a list of fault features, evaluates the contact
+        points on each fault and compute mean value and standard
+        deviation. These values will be used to identify each
+        fault with the indicator function.
 
         Parameters
         ----------
@@ -232,35 +236,34 @@ def add_faults_anisotropies(self, fault_list=None):
         -------
 
         """
-        if fault_list == None:
-            self.anisotropies_fault_list = []
-        else:
-            self.anisotropies_fault_list = fault_list
-            faults_parameters = {}
-            for i in range(len(fault_list)):
+        self.anisotropis_fault_list = fault_list
 
-                data_temp = self.intrusion_network_data[
-                    self.intrusion_network_data["intrusion_anisotropy"]
-                    == fault_list[i].name
-                ].copy()
-                data_array_temp = data_temp.loc[:, ["X", "Y", "Z"]].to_numpy()
+        faults_parameters = {}
+        for fault in fault_list:
 
-                if data_temp.empty:
-                    fault_i_mean = 0
-                    fault_i_std = 0.1
+            data_temp = self.intrusion_network_data[
+                self.intrusion_network_data["intrusion_anisotropy"] == fault.name
+            ].copy()
+            data_array_temp = data_temp.loc[:, ["X", "Y", "Z"]].to_numpy()
 
-                else:
-                    fault_i_vals = fault_list[i][0].evaluate_value(data_array_temp)
-                    fault_i_mean = np.mean(fault_i_vals)
-                    fault_i_std = np.std(fault_i_vals)
-
-                faults_parameters[fault_list[i].name] = [
-                    fault_list[i],
-                    fault_i_mean,
-                    fault_i_std,
-                ]
+            if data_temp.empty:
+                fault_i_mean = 0
+                fault_i_std = 0.1
 
-            self.anisotropies_fault_parameters = faults_parameters
+            else:
+                # evaluate the value of the fault, evaluate_value called on a fault
+                # will return the scalar field value of the main structural feature
+                fault_i_vals = fault.evaluate_value(data_array_temp)
+                fault_i_mean = np.mean(fault_i_vals)
+                fault_i_std = np.std(fault_i_vals)
+
+            faults_parameters[fault.name] = [
+                fault,
+                fault_i_mean,
+                fault_i_std,
+            ]
+
+        self.anisotropies_fault_parameters = faults_parameters
 
     def add_steps_stratigraphic_units(self):
         """
@@ -277,9 +280,11 @@ def add_steps_stratigraphic_units(self):
 
         """
 
-        if self.model.stratigraphic_column == None:
+        if self.model.stratigraphic_column is None:
             logger.error(
-                "Set stratigraphic column to model using model.set_stratigraphic_column(name_of_stratigraphic_column)"
+                "No stratigraphic column. \
+                Set stratigraphic column using \
+                model.set_stratigraphic_column(name_of_stratigraphic_column)"
             )
 
         # set data
@@ -289,16 +294,16 @@ def add_steps_stratigraphic_units(self):
         )
         std_backup = 25
 
-        for step_i in self.intrusion_steps.keys():
-            step_structure = self.intrusion_steps[step_i].get("structure")
-            unit_from_name = self.intrusion_steps[step_i].get("unit_from")
-            series_from_name = self.intrusion_steps[step_i].get("series_from")
+        for step_i, step in self.intrusion_steps.items():
+            step_structure = step.get("structure")
+            unit_from_name = step.get("unit_from")
+            series_from_name = step.get("series_from")
             unit_from_id = self.model.stratigraphic_column[series_from_name.name][
                 unit_from_name
             ].get("id")
 
-            unit_to_name = self.intrusion_steps[step_i].get("unit_to")
-            series_to_name = self.intrusion_steps[step_i].get("series_to")
+            unit_to_name = step.get("unit_to")
+            series_to_name = step.get("series_to")
             unit_to_id = self.model.stratigraphic_column[series_to_name.name][
                 unit_to_name
             ].get("id")
@@ -346,16 +351,16 @@ def add_steps_stratigraphic_units(self):
                 contact_1_std = np.std(contact_1_vals)
 
                 if contact_0_mean <= contact_1_mean:
-                    self.intrusion_steps[step_i]["unit_from_mean"] = contact_0_mean
-                    self.intrusion_steps[step_i]["unit_from_std"] = contact_0_std
-                    self.intrusion_steps[step_i]["unit_to_mean"] = contact_1_mean
-                    self.intrusion_steps[step_i]["unit_to_std"] = contact_1_std
+                    step["unit_from_mean"] = contact_0_mean
+                    step["unit_from_std"] = contact_0_std
+                    step["unit_to_mean"] = contact_1_mean
+                    step["unit_to_std"] = contact_1_std
 
                 else:
-                    self.intrusion_steps[step_i]["unit_from_mean"] = contact_1_mean
-                    self.intrusion_steps[step_i]["unit_from_std"] = contact_1_std
-                    self.intrusion_steps[step_i]["unit_to_mean"] = contact_0_mean
-                    self.intrusion_steps[step_i]["unit_to_std"] = contact_0_std
+                    step["unit_from_mean"] = contact_1_mean
+                    step["unit_from_std"] = contact_1_std
+                    step["unit_to_mean"] = contact_0_mean
+                    step["unit_to_std"] = contact_0_std
 
             else:  # step between different strat units
                 data_points_from_xyz = (
@@ -377,10 +382,10 @@ def add_steps_stratigraphic_units(self):
                     unit_from_max = self.model.stratigraphic_column[
                         series_from_name.name
                     ][unit_from_name].get("max")
-                    self.intrusion_steps[step_i]["unit_from_mean"] = (
+                    step["unit_from_mean"] = (
                         unit_from_min + (unit_from_max - unit_from_min) / 2
                     )
-                    self.intrusion_steps[step_i]["unit_from_std"] = std_backup
+                    step["unit_from_std"] = std_backup
                 else:
                     series_values = series_from_name.evaluate_value(
                         data_points_from_xyz
@@ -395,15 +400,11 @@ def add_steps_stratigraphic_units(self):
                         )
                     else:
                         series_values_mod = series_values
-                    self.intrusion_steps[step_i]["unit_from_mean"] = np.nanmean(
-                        series_values_mod
-                    )
-                    self.intrusion_steps[step_i]["unit_from_std"] = np.nanstd(
-                        series_values_mod
-                    )
+                    step["unit_from_mean"] = np.nanmean(series_values_mod)
+                    step["unit_from_std"] = np.nanstd(series_values_mod)
 
-                if self.intrusion_steps[step_i]["unit_from_std"] == 0:
-                    self.intrusion_steps[step_i]["unit_from_std"] = std_backup
+                if step["unit_from_std"] == 0:
+                    step["unit_from_std"] = std_backup
 
                 data_points_to_xyz = (
                     intrusion_network_data[
@@ -423,10 +424,8 @@ def add_steps_stratigraphic_units(self):
                     unit_to_max = self.model.stratigraphic_column[series_to_name.name][
                         unit_to_name
                     ].get("max")
-                    self.intrusion_steps[step_i]["unit_to_mean"] = (
-                        unit_to_min + (unit_to_max - unit_to_min) / 2
-                    )
-                    self.intrusion_steps[step_i]["unit_to_std"] = std_backup
+                    step["unit_to_mean"] = unit_to_min + (unit_to_max - unit_to_min) / 2
+                    step["unit_to_std"] = std_backup
                 else:
                     series_values = series_to_name.evaluate_value(data_points_to_xyz)
                     # print(len(step_structure_points_vals), len(series_values))
@@ -439,15 +438,11 @@ def add_steps_stratigraphic_units(self):
                         )
                     else:
                         series_values_mod = series_values
-                    self.intrusion_steps[step_i]["unit_to_mean"] = np.nanmean(
-                        series_values_mod
-                    )
-                    self.intrusion_steps[step_i]["unit_to_std"] = np.nanstd(
-                        series_values_mod
-                    )
+                    step["unit_to_mean"] = np.nanmean(series_values_mod)
+                    step["unit_to_std"] = np.nanstd(series_values_mod)
 
-                if self.intrusion_steps[step_i]["unit_to_std"] == 0:
-                    self.intrusion_steps[step_i]["unit_to_std"] = std_backup
+                if step["unit_to_std"] == 0:
+                    step["unit_to_std"] = std_backup
 
     def set_intrusion_network_parameters(
         self, intrusion_data, intrusion_network_input, **kwargs
@@ -565,13 +560,16 @@ def indicator_function_contacts(self, delta=None):
         Function to compute indicator function for list of contacts anisotropies
         For each point of the grid, this function assignes a 1 if contact i is present, 0 otherwise.
 
-        A contact is defined as an isovalue of an scalar field defining the geological feature  of which the contact is part of.
-        Each point is evaluated in the feature scalar field and is identified as part of the contact if its value is around the contact isovalue.
+        A contact is defined as an isovalue of an scalar field defining the
+        geological feature  of which the contact is part of.
+        Each point is evaluated in the feature scalar field and is
+        identified as part of the contact if its value is around the contact isovalue.
 
         Parameters
         ----------
         delta : list of numbers, same lenght as number of anisotropies (series).
-            delta multiplies the standard deviation to increase probability of finding the contact on a grid point.
+            delta multiplies the standard deviation to increase
+            probability of finding the contact on a grid point.
 
         Returns
         ----------
@@ -613,8 +611,10 @@ def indicator_function_faults(self, delta=None):
         Function to compute indicator function for list of faults anisotropies
         For each point of the grid, this function assignes a 1 if fault i is present, 0 otherwise.
 
-        A fault surface is defined as an isovalue 0 of the scalar field representing the fault surface (coordinate 0 of its structural frame)
-        Each point of the grid is evaluated in this scalar field and is identified as part of the fault if its value is around 0.
+        A fault surface is defined as an isovalue 0 of the scalar field representing
+         the fault surface (coordinate 0 of its structural frame)
+        Each point of the grid is evaluated in this scalar field
+        and is identified as part of the fault if its value is around 0.
 
         Parameters
         ----------
@@ -653,7 +653,8 @@ def indicator_function_faults(self, delta=None):
     def compute_velocity_field(self, indicator_fx_contacts, indicator_fx_faults):
         """
         Function to compute velocity field to be used in the shortest path search.
-        A velocity parameter is assign to each point of the grid, depending on which anisotropy is found at each grid point
+        A velocity parameter is assign to each point of the grid, depending on
+        which anisotropy is found at each grid point
 
         Velocity field
             K(x) = sum(Ic+Kc) + (1-sum(Ic))*sum(If*Kf)
@@ -723,7 +724,7 @@ def create_intrusion_network(self, **kwargs):
         """
 
         # --- check type of intrusion network
-        if self.intrusion_network_type == None:
+        if not self.intrusion_network_type:
             logger.error(
                 "Specify type of intrusion network: 'interpolated' or 'shortest path'"
             )
@@ -754,8 +755,6 @@ def create_intrusion_network(self, **kwargs):
 
                 # evaluate grid points in series
 
-                faults = []
-
                 for i, step_i in enumerate(self.intrusion_steps.keys()):
                     delta_contact = self.delta_contacts[i]
                     if type(delta_contact) is list:
@@ -769,20 +768,25 @@ def create_intrusion_network(self, **kwargs):
                     series_from_name = self.intrusion_steps[step_i].get("series_from")
                     series_to_name = self.intrusion_steps[step_i].get("series_to")
 
+                    # TODO LG to as FA why this is done? aren't faults
+                    # enabled by default? If they ar disabled, should
+                    # they get turned back off after this?
                     series_from_name.faults_enabled = True
                     series_to_name.faults_enabled = True
 
                     fault_gridpoints_vals = step_fault[0].evaluate_value(grid_points)
-                    fault_gridpoints_grads = step_fault[0].evaluate_gradient(
+
+                    series_from_gridpoints_vals = series_from_name.evaluate_value(
                         grid_points
                     )
-
-                    series_from_gridpoints_vals = series_from_name.evaluate_value(grid_points)
+                    # if
                     if series_from_name == series_to_name:
                         series_to_gridpoints_vals = series_from_gridpoints_vals
 
                     else:
-                        series_to_gridpoints_vals = series_to_name.evaluate_value(grid_points)
+                        series_to_gridpoints_vals = series_to_name.evaluate_value(
+                            grid_points
+                        )
 
                     contacts0_val_min = self.intrusion_steps[step_i].get(
                         "unit_from_mean"
@@ -951,7 +955,7 @@ def create_intrusion_network(self, **kwargs):
 
             # create dataframe containing grid points, velocity values
             nodes = np.linspace(0, len(grid_points), len(grid_points))
-            medium = np.zeros([len(grid_points), 7])
+            medium = np.zeros((len(grid_points), 7))
             medium[:, 0] = nodes.T
             medium[:, 1:4] = grid_points
             medium[:, 4] = velocity_field
@@ -1121,7 +1125,7 @@ def set_intrusion_frame_data(self, intrusion_frame_data, intrusion_network_point
         scaled_inet_points = intrusion_network_points[:, :3]
         coord_0_values = pd.DataFrame(scaled_inet_points, columns=["X", "Y", "Z"])
         coord_0_values["val"] = 0
-        coord_0_values["coord"] = 0 
+        coord_0_values["coord"] = 0
         coord_0_values["feature_name"] = self.name
         coord_0_values["w"] = 1
 
diff --git a/LoopStructural/modelling/intrusions/intrusion_support_functions.py b/LoopStructural/modelling/intrusions/intrusion_support_functions.py
index 805a2ea9d..c4b8572fe 100644
--- a/LoopStructural/modelling/intrusions/intrusion_support_functions.py
+++ b/LoopStructural/modelling/intrusions/intrusion_support_functions.py
@@ -7,6 +7,53 @@
 logger = getLogger(__name__)
 
 
+def sort_2_arrays(main_array, array):
+    """
+    Sort two arrays, considering values of only the main array
+
+    Parameters
+    ----------
+    main array: numpy array, considered to sort secondary array
+    array: numpy aray, array to be sorted
+
+    Returns
+    -------
+    sorted arrays
+    """
+    # function to sort 2 arrays, considering values of only the main array
+
+    for i in range(len(main_array)):
+        swap = i + np.argmin(main_array[i:])
+        (main_array[i], main_array[swap]) = (main_array[swap], main_array[i])
+        (array[i], array[swap]) = (array[swap], array[i])
+
+    return main_array, array
+
+
+def findMinDiff(arr, n):
+    """
+    Find the min diff by comparing difference of all possible pairs in given array
+
+    Parameters
+    ----------
+    arr: numpy array with values to compare and fin the minimum difference
+
+    Returns
+    -------
+    minimum difference between values in arr
+
+    """
+    # Initialize difference as infinite
+    diff = 10 ** 20
+
+    for i in range(n - 1):
+        for j in range(i + 1, n):
+            if abs(arr[i] - arr[j]) < diff:
+                diff = abs(arr[i] - arr[j])
+
+    return diff
+
+
 def array_from_coords(df, section_axis, df_axis):
     """
     Create numpy array representing a section of the model
@@ -103,6 +150,250 @@ def find_inout_points(velocity_field_array, velocity_parameters):
     return inlet_point, outlet_point
 
 
+def shortest_path(inlet, outlet, time_map):
+    """
+    Look for the shortest path between inlet and outlet, using a time map.
+    In practice, for a given point looks for the neighbour elements which is closer in time.
+    The search starts in the inlet, until it reaches the outlet.
+
+    Parameters
+    ----------
+    inlet: list of indexes (row and column) of inlet point.
+    outlet: list of indexes (row and column) of outlet point.
+    time_map: numpy array, contains values of time, computed using the fast-marching method.
+
+    Returns
+    -------
+    inet: (intrusion network) numpy array of same shape of time_map array.
+        0 where the intrusion network is found, 1 above it, and -1 below it
+
+    """
+    inet = np.ones_like(time_map)  # array to save shortest path with zeros
+    temp_inlet = inlet  # temporary inlet
+    inet[temp_inlet[0], temp_inlet[1]] = 0
+    i = 0
+
+    n_rows = len(time_map)
+    n_cols = len(time_map[0])
+
+    while True:
+        i = i + 1
+        time_temp_inlet = time_map[
+            temp_inlet[0], temp_inlet[1]
+        ]  # obtain time value of temporary outlet
+        neighbors = element_neighbour(
+            temp_inlet, time_map, inet
+        )  # identify neighbours elements of temporary outlet
+        direction = index_min(
+            neighbors
+        )  # obtain the location (index min) of minimun difference
+
+        if direction == 10:
+            break
+
+        temp_inlet = new_inlet(temp_inlet, direction)
+        row = temp_inlet[0]
+        col = temp_inlet[1]
+
+        # if row >= n_rows or col >= n_cols:
+        #     break
+        # else:
+        inet[row, col] = 0
+
+        if temp_inlet[0] == outlet[0] and temp_inlet[1] == outlet[1]:
+            break
+        else:
+            continue
+
+    # Assing -1 to points below intrusion network
+    for j in range(len(inet[0])):  # columns
+        for h in range(len(inet)):  # rows
+            if inet[h, j] == 0:
+                break
+
+        inet[(h + 1) :, j] = -1
+
+    return inet
+
+
+def element_neighbour(index, array, inet):
+    """
+    Identify the value of neighbours elements for a given element.
+
+    Parameters
+    ----------
+    index: list of indexes (row and column) of elements.
+    array: numpy array, containing values
+    inet: numpy array, temporal intrusion network. If one of the elements is already inet=0, the assign -1.
+
+    Returns
+    -------
+    values: numpy array, 1x5 with the values of the neighbours of a particular element
+
+    """
+
+    rows = len(array) - 1  # max index of rows of time_map array
+    cols = len(array[0]) - 1  # max index of columns of time_map arrays
+    values = np.zeros(
+        8
+    )  # 8 - array to save values (element above, element to the left, element to the right)
+    # values[8] = 10
+    index_row = index[0]
+    index_col = index[1]
+
+    if index_row == 0:
+        values[0] = -1
+        values[1] = -1
+        values[2] = -1
+
+    if index_row == rows:
+        values[5] = -1
+        values[6] = -1
+        values[7] = -1
+
+    if index_col == 0:
+        values[0] = -1
+        values[3] = -1
+        values[5] = -1
+
+    if index_col == cols:
+        values[2] = -1
+        values[4] = -1
+        values[7] = -1
+
+    for k in range(8):
+        if values[k] > -1:
+            if k == 0:
+                values[0] = array[index[0] - 1, index[1] - 1]
+
+            if k == 1:
+                values[1] = array[index[0] - 1, index[1]]
+
+            if k == 2:
+                values[2] = array[index[0] - 1, index[1] + 1]
+
+            if k == 3:
+                values[3] = array[index[0], index[1] - 1]
+
+            if k == 4:
+                values[4] = array[index[0], index[1] + 1]
+
+            if k == 5:
+                values[5] = array[index[0] + 1, index[1] - 1]
+
+            if k == 6:
+                values[6] = array[index[0] + 1, index[1]]
+
+            if k == 7:
+                values[7] = array[index[0] + 1, index[1] + 1]
+
+        else:
+            continue
+
+    # check if some of the neighbours is already part of the intrusion network
+    for h in range(8):
+        if values[h] > -1:
+            if h == 0:
+                if inet[index[0] - 1, index[1] - 1] == 0:
+                    values[0] = -2
+            if h == 1:
+                if inet[index[0] - 1, index[1]] == 0:
+                    values[1] = -2
+            if h == 2:
+                if inet[index[0] - 1, index[1] + 1] == 0:
+                    values[2] = -2
+            if h == 3:
+                if inet[index[0], index[1] - 1] == 0:
+                    values[3] = -2
+            if h == 4:
+                if inet[index[0], index[1] + 1] == 0:
+                    values[4] = -2
+            if h == 5:
+                if inet[index[0] + 1, index[1] - 1] == 0:
+                    values[5] = -2
+            if h == 6:
+                if inet[index[0] + 1, index[1]] == 0:
+                    values[6] = -2
+            if h == 7:
+                if inet[index[0] + 1, index[1] + 1] == 0:
+                    values[7] = -2
+        else:
+            continue
+
+    return values
+
+
+def index_min(array):
+    """
+    Given an array of 1x8, dentify the index of the minimum value within the array.
+
+    Parameters
+    ----------
+    array: numpy array, 1x8
+
+    Returns
+    -------
+    index_min: integer, index of minimum value in array
+
+    """
+
+    # return the index value of the minimum value in an array of 1x8
+    # print(array)
+    index_array = {}
+
+    for i in range(
+        8
+    ):  # create a dictionary assining positions from 0 to 7 to the values in the array
+        a = i
+        if array[i] >= 0:
+            index_array.update({i: array[i]})
+
+    if len(index_array.values()) > 0:
+
+        minimum_val = min(index_array.values())
+
+        for key, value in index_array.items():
+            if value == minimum_val:
+                index_min = key
+
+    else:
+        index_min = 10
+
+    return index_min
+
+
+def new_inlet(inlet, direction):
+    """
+    Determine new inlet indexes, given current inlet and direction of minimum difference in time map.
+
+    Parameters
+    ----------
+    inlet: list of indexes [row, column]
+    direction: integers e[0,7] (0: above-left, 1: above, 2: above right, 3: left, 4: right, 5: below left, 6: below, 7:below right)
+
+    Returns
+    -------
+    new_inlet: list of indexes [row, column]
+
+    """
+    pot_new_inlets = {}
+
+    pot_new_inlets.update({"0": np.array([inlet[0] - 1, inlet[1] - 1])})
+    pot_new_inlets.update({"1": np.array([inlet[0] - 1, inlet[1]])})
+    pot_new_inlets.update({"2": np.array([inlet[0] - 1, inlet[1] + 1])})
+    pot_new_inlets.update({"3": np.array([inlet[0], inlet[1] - 1])})
+    pot_new_inlets.update({"4": np.array([inlet[0], inlet[1] + 1])})
+    pot_new_inlets.update({"5": np.array([inlet[0] + 1, inlet[1] - 1])})
+    pot_new_inlets.update({"6": np.array([inlet[0] + 1, inlet[1]])})
+    pot_new_inlets.update({"7": np.array([inlet[0] + 1, inlet[1] + 1])})
+    new_inlet = np.zeros(2)
+
+    for key, value in pot_new_inlets.items():
+        if key == str(direction):
+            new_inlet = value
+    return new_inlet
+
+
 def grid_from_array(array, fixed_coord, lower_extent, upper_extent):
 
     """