diff --git a/code/GP_Code.py b/code/GP_Code.py
index 681de80..13d87b3 100644
--- a/code/GP_Code.py
+++ b/code/GP_Code.py
@@ -1,11 +1,10 @@
 # -*- coding: utf-8 -*-
 """
-Created on Tue Oct 22 11:43:41 2024
+Created on Tue Oct 22 15:21:41 2024
 
 @author: diana
 """
 
-
 import glob
 import os
 import rlxnix as rlx
@@ -17,8 +16,8 @@ from scipy.integrate import quad
 
 ### FUNCTIONS ###
 def binary_spikes(spike_times, duration, dt): 
-    """ Converts the spike times to a binary representation. 
-    Zeros when there is no spike, One when there is.     
+    """Converts the spike times to a binary representation. 
+    Zeros when there is no spike, one when there is.
 
     Parameters
     ----------
@@ -27,15 +26,14 @@ def binary_spikes(spike_times, duration, dt):
     duration : float
         The trial duration.
     dt : float
-        the temporal resolution.
+        The temporal resolution.
 
     Returns
     -------
     binary : np.array
         The binary representation of the spike times.
-
     """
-    binary = np.zeros(int(np.round(duration / dt))) #Vektor, der genauso lang ist wie die stim time
+    binary = np.zeros(int(np.round(duration / dt)))  #Vektor, der genauso lang ist wie die stim time
     spike_indices = np.asarray(np.round(spike_times / dt), dtype=int)
     binary[spike_indices] = 1
     return binary
@@ -44,9 +42,6 @@ def binary_spikes(spike_times, duration, dt):
 def firing_rate(binary_spikes, box_width, dt=0.000025):
     """Calculate the firing rate from binary spike data.
 
-    This function computes the firing rate using a boxcar (moving average) 
-    filter of a specified width.
-
     Parameters
     ----------
     binary_spikes : np.array
@@ -62,7 +57,7 @@ def firing_rate(binary_spikes, box_width, dt=0.000025):
         An array representing the firing rate at each time step.
     """
     box = np.ones(int(box_width // dt))
-    box /= np.sum(box) * dt #Normalization of box kernel to an integral of 1
+    box /= np.sum(box) * dt  # Normalization of box kernel to an integral of 1
     rate = np.convolve(binary_spikes, box, mode="same")
     return rate
 
@@ -90,7 +85,96 @@ def powerspectrum(rate, dt):
     return frequency, power
 
 
+def calculate_integral(frequency, power, point, delta):   
+    """
+    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
+        Half-width of the range for integration around the point.
+
+    Returns
+    -------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    """
+    indices = (frequency >= point - delta) & (frequency <= point + delta)
+    integral = np.trapz(power[indices], frequency[indices])
+            
+    left_indices = (frequency >= point - 5 * delta) & (frequency < point - delta)
+    right_indices = (frequency > point + delta) & (frequency <= point + 5 * delta)
+           
+    l_integral = np.trapz(power[left_indices], frequency[left_indices])
+    r_integral = np.trapz(power[right_indices], frequency[right_indices])
+           
+    local_mean = np.mean([l_integral, r_integral])
+    return integral, local_mean
+
+
+def valid_integrals(integral, local_mean, threshold, point):
+    """
+    Check if the integral exceeds the threshold compared to the local mean and 
+    provide feedback on whether the given point is valid or not.
+
+    Parameters
+    ----------
+    integral : float
+        The calculated integral around the point.
+    local_mean : float
+        The local mean value (adjacent integrals).
+    threshold : float
+        Threshold value to compare integrals with local mean.
+    point : float
+        The harmonic frequency point being evaluated.
+
+    Returns
+    -------
+    valid : bool
+        True if the integral exceeds the local mean by the threshold, otherwise False.
+    message : str
+        A message stating whether the point is valid or not.
+    """
+    valid = integral > (local_mean * threshold)
+    if valid:
+        message = f"The point {point} is valid, as its integral exceeds the threshold."
+    else:
+        message = f"The point {point} is not valid, as its integral does not exceed the threshold."
+    return valid, message
+
+
 def prepare_harmonics(frequencies, categories, num_harmonics, colors):
+    """
+    Prepare harmonic frequencies and assign colors based on categories.
+    
+    Parameters
+    ----------
+    frequencies : list
+        Base frequencies to generate harmonics.
+    categories : list
+        Corresponding categories for the base frequencies.
+    num_harmonics : list
+        Number of harmonics for each base frequency.
+    colors : list
+        List of colors corresponding to the categories.
+
+    Returns
+    -------
+    points : list
+        A flat list of harmonic frequencies.
+    color_mapping : dict
+        A dictionary mapping each category to its corresponding color.
+    points_categories : dict
+        A mapping of categories to their harmonic frequencies.
+    """
     points_categories = {}
     for idx, (freq, category) in enumerate(zip(frequencies, categories)):
         points_categories[category] = [freq * (i + 1) for i in range(num_harmonics[idx])]
@@ -101,12 +185,9 @@ def prepare_harmonics(frequencies, categories, num_harmonics, colors):
     return points, color_mapping, points_categories
 
 
-def plot_power_spectrum_with_integrals(frequency, power, points, delta):
+def find_exceeding_points(frequency, power, points, delta, threshold):
     """
-    Create a figure of the power spectrum and calculate integrals around specified points.
-
-    This function generates the plot of the power spectrum and calculates integrals
-    around specified harmonic points, but it does not color the regions or add vertical lines.
+    Find the points where the integral exceeds the local mean by a given threshold.
 
     Parameters
     ----------
@@ -115,60 +196,35 @@ def plot_power_spectrum_with_integrals(frequency, power, points, delta):
     power : np.array
         An array of power spectral density values.
     points : list
-        A list of harmonic frequencies to highlight.
+        A list of harmonic frequencies to evaluate.
     delta : float
-        Half-width of the range for integration around each point.
+        Half-width of the range for integration around the point.
+    threshold : float
+        Threshold value to compare integrals with local mean.
 
     Returns
     -------
-    integrals : list
-        List of calculated integrals for each point.
-    local_means : list
-        List of local mean values (adjacent integrals).
-    fig : matplotlib.figure.Figure
-        The created figure object with the power plot.
-    ax : matplotlib.axes.Axes
-        The axes object for further modifications.
+    exceeding_points : list
+        A list of points where the integral exceeds the local mean by the threshold.
     """
+    exceeding_points = []
     
-    fig, ax = plt.subplots()
-    ax.plot(frequency, power)  # Plot power spectrum
-    
-    integrals = []
-    local_means = []
-
     for point in points:
-        # Define indices for the integration window
-        indices = (frequency >= point - delta) & (frequency <= point + delta)
-        # Calculate integral around the point
-        integral = np.trapz(power[indices], frequency[indices])
-        integrals.append(integral)
-        
-        # Calculate adjacent region integrals for local mean
-        left_indices = (frequency >= point - 5 * delta) & (frequency < point - delta)
-        right_indices = (frequency > point + delta) & (frequency <= point + 5 * delta)
+        # Calculate the integral and local mean for the current point
+        integral, local_mean = calculate_integral(frequency, power, point, delta)
         
-        l_integral = np.trapz(power[left_indices], frequency[left_indices])
-        r_integral = np.trapz(power[right_indices], frequency[right_indices])
+        # Check if the integral exceeds the threshold
+        valid, message = valid_integrals(integral, local_mean, threshold, point)
         
-        local_mean = np.mean([l_integral, r_integral])
-        local_means.append(local_mean)
-
-    ax.set_xlim([0, 1200])  # Set x-axis limit
-    ax.set_xlabel('Frequency (Hz)')
-    ax.set_ylabel('Power')
-    ax.set_title('Power Spectrum with Integrals (Uncolored)')
-
-    return integrals, local_means, fig, ax
+        if valid:
+            exceeding_points.append(point)
+    
+    return exceeding_points
 
 
-def highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig_orig, ax_orig):
+def plot_highlighted_integrals(frequency, power, exceeding_points, delta, threshold, color_mapping, points_categories):
     """
-    Create a new figure by highlighting integrals that exceed the threshold.
-
-    This function generates a new figure with colored shading around points where the integrals exceed
-    the local mean by a given threshold and adds vertical lines at the boundaries of adjacent regions.
-    It leaves the original figure unchanged.
+    Plot the power spectrum and highlight integrals that exceed the threshold.
 
     Parameters
     ----------
@@ -176,75 +232,62 @@ def highlight_integrals_with_threshold(frequency, power, points, delta, threshol
         An array of frequencies corresponding to the power values.
     power : np.array
         An array of power spectral density values.
-    points : list
-        A list of harmonic frequencies to highlight.
+    exceeding_points : list
+        A list of harmonic frequencies that exceed the threshold.
     delta : float
         Half-width of the range for integration around each point.
     threshold : float
         Threshold value to compare integrals with local mean.
-    integrals : list
-        List of calculated integrals for each point.
-    local_means : list
-        List of local mean values (adjacent integrals).
     color_mapping : dict
-        A mapping of point categories to colors.
+        A dictionary mapping each category to its color.
     points_categories : dict
         A mapping of categories to lists of points.
-    fig_orig : matplotlib.figure.Figure
-        The original figure object (remains unchanged).
-    ax_orig : matplotlib.axes.Axes
-        The original axes object (remains unchanged).
 
     Returns
     -------
-    fig_new : matplotlib.figure.Figure
-        The new figure object with color highlights and vertical lines.
+    fig : matplotlib.figure.Figure
+        The created figure object with highlighted integrals.
     """
+    fig, ax = plt.subplots()
+    ax.plot(frequency, power)  # Plot power spectrum
     
-    # Create a new figure based on the original power spectrum
-    fig_new, ax_new = plt.subplots()
-    ax_new.plot(frequency, power)  # Plot the same power spectrum
-    
-    # Loop through each point and check if the integral exceeds the threshold
-    for i, point in enumerate(points):
-        exceeds = integrals[i] > (local_means[i] * threshold)
-
-        if exceeds:
+    for point in exceeding_points:
+        integral, local_mean = calculate_integral(frequency, power, point, delta)
+        valid, _ = valid_integrals(integral, local_mean, threshold, point)
+        if valid:  
             # Define color based on the category of the point
             color = next((c for cat, c in color_mapping.items() if point in points_categories[cat]), 'gray')
             # Shade the region around the point where the integral was calculated
-            ax_new.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
-            print(f"Integral around {point:.2f} Hz: {integrals[i]:.5e}")
-            
+            ax.axvspan(point - delta, point + delta, color=color, alpha=0.3, label=f'{point:.2f} Hz')
+            print(f"Integral around {point:.2f} Hz: {integral:.5e}")
+
+                        
             # Define left and right boundaries of adjacent regions
             left_boundary = frequency[np.where((frequency >= point - 5 * delta) & (frequency < point - delta))[0][0]]
             right_boundary = frequency[np.where((frequency > point + delta) & (frequency <= point + 5 * delta))[0][-1]]
 
             # Add vertical dashed lines at the boundaries of the adjacent regions
-            ax_new.axvline(x=left_boundary, color="k", linestyle="--")
-            ax_new.axvline(x=right_boundary, color="k", linestyle="--")
-
-    # Update plot legend and return the new figure
-    ax_new.set_xlim([0, 1200])
-    ax_new.set_xlabel('Frequency (Hz)')
-    ax_new.set_ylabel('Power')
-    ax_new.set_title('Power Spectrum with Highlighted Integrals')
-    ax_new.legend()
+            ax.axvline(x=left_boundary, color="k", linestyle="--")
+            ax.axvline(x=right_boundary, color="k", linestyle="--")
+            
+            
+    ax.set_xlim([0, 1200])
+    ax.set_xlabel('Frequency (Hz)')
+    ax.set_ylabel('Power')
+    ax.set_title('Power Spectrum with Highlighted Integrals')
+    ax.legend()
     
-    return fig_new
-
-
-
+    return fig
 
 
 ### Data retrieval ###
-datafolder = "../data" # Geht in der Hierarchie einen Ordern nach oben (..) und dann in den Ordner 'data'
+datafolder = "../data"
 example_file = os.path.join("..", "data", "2024-10-16-ad-invivo-1.nix")
 dataset = rlx.Dataset(example_file)
 sams = dataset.repro_runs("SAM") 
 sam = sams[2]
 
-## Daten für Funktionen
+## Data for functions
 df = sam.metadata["RePro-Info"]["settings"]["deltaf"][0][0]
 stim = sam.stimuli[1]
 potential, time = stim.trace_data("V-1")
@@ -253,29 +296,34 @@ duration = stim.duration
 dt = stim.trace_info("V-1").sampling_interval
 
 
-### Anwendung Functionen ### 
+### Apply Functions to calculate data ###
 b = binary_spikes(spikes, duration, dt)
 rate = firing_rate(b, box_width=0.05, dt=dt)
 frequency, power = powerspectrum(b, dt)
 
-## Important stuff 
+
+### Important stuff ###
+## Frequencies
 eodf = stim.metadata[stim.name]["EODf"][0][0]
 stimulus_frequency = eodf + df
-AM = 50 # Hz
-#print(f"EODf: {eodf}, Stimulus Frequency: {stimulus_frequency}, AM: {AM}")
-
+AM = 50  # Hz
 frequencies = [AM, eodf, stimulus_frequency]
+
 categories = ["AM", "EODf", "Stimulus frequency"]
 num_harmonics = [4, 2, 2]
 colors = ["green", "orange", "red"]
+
 delta = 2.5
 threshold = 10
 
-### 
-points, color_mapping, points_categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
-
-# First, create the power spectrum plot with integrals (without coloring)
-integrals, local_means, fig1, ax1 = plot_power_spectrum_with_integrals(frequency, power, points, delta)
-
-# Then, create a new separate figure where integrals exceeding the threshold are highlighted
-fig2 = highlight_integrals_with_threshold(frequency, power, points, delta, threshold, integrals, local_means, color_mapping, points_categories, fig1, ax1)
\ No newline at end of file
+### Apply functions to make powerspectrum ###
+integral, local = calculate_integral(frequency, power, eodf, delta)
+valid = valid_integrals(integral, local, threshold, eodf)
+points, color, categories = prepare_harmonics(frequencies, categories, num_harmonics, colors)
+print(len(points))
+exceeding = find_exceeding_points(frequency, power, points, delta, threshold)
+print(len(exceeding))
+
+## Plot power spectrum and highlight integrals
+fig = plot_highlighted_integrals(frequency, power, points, delta, threshold, color, categories)
+plt.show()