diff --git a/code/useful_functions.py b/code/useful_functions.py
index ed94b8e..30dbd9c 100644
--- a/code/useful_functions.py
+++ b/code/useful_functions.py
@@ -47,7 +47,7 @@ def all_coming_together(freq_array, power_array, points_list, categories, num_ha
         color = colors[i]
         
         # Step 1: Calculate the integral for the point
-        integral, local_mean, _ = calculate_integral(freq_array, power_array, point, delta)
+        integral, local_mean = calculate_integral_2(freq_array, power_array, point, delta)
         
         # Step 2: Check if the point is valid
         valid = valid_integrals(integral, local_mean, point, threshold)
@@ -150,6 +150,42 @@ def calculate_integral(freq, power, point, delta = 2.5):
     local_mean = np.mean([l_integral, r_integral])
     return integral, local_mean, p_power
 
+def calculate_integral_2(freq, power, point, delta = 2.5):   
+    """
+    Calculate the integral around a single specified point.
+
+    Parameters
+    ----------
+    frequency : np.array
+        An array of frequencies corresponding to the power values.
+    power : np.array
+        An array of power spectral density values.
+    point : float
+        The harmonic frequency at which to calculate the integral.
+    delta : float, optional
+        Radius of the range for integration around the point. The default is 2.5.
+
+    Returns
+    -------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    p_power : float
+        The local maxiumum power.
+    """
+    indices = (freq >= point - delta) & (freq <= point + delta)
+    integral = np.trapz(power[indices], freq[indices])
+            
+    left_indices = (freq >= point - 5 * delta) & (freq < point - delta)
+    right_indices = (freq > point + delta) & (freq <= point + 5 * delta)
+           
+    l_integral = np.trapz(power[left_indices], freq[left_indices])
+    r_integral = np.trapz(power[right_indices], freq[right_indices])
+           
+    local_mean = np.mean([l_integral, r_integral])
+    return integral, local_mean
+
 def contrast_sorting(sams, con_1 = 20, con_2 = 10, con_3 = 5, stim_count = 3, stim_dur = 2):
     '''
     sorts the sams into three contrasts