diff --git a/Source/main_yael.py b/Source/TODEL_main_yael.py
similarity index 98%
rename from Source/main_yael.py
rename to Source/TODEL_main_yael.py
index d233716..82f551d 100644
--- a/Source/main_yael.py
+++ b/Source/TODEL_main_yael.py
@@ -98,7 +98,7 @@ ax.set_xlabel("energy E")
 ax.set_ylabel("quantity mu")
 
 fig.tight_layout()
-fig.savefig("mu")
+# fig.savefig("mu")
 
 fig, ax = plt.subplots(1)
 ax.scatter(all_E, A, color="C0", s=5, marker="o", label="Data")
@@ -106,7 +106,7 @@ ax.set_xlabel("energy E")
 ax.set_ylabel("area A")
 ax.legend()
 fig.tight_layout()
-fig.savefig("area")
+# fig.savefig("area")
 plt.show()
 
 
diff --git a/Source/poincare_test.py b/Source/TODEL_poincare_test.py
similarity index 100%
rename from Source/poincare_test.py
rename to Source/TODEL_poincare_test.py
diff --git a/Source/energies.py b/Source/energies.py
index d3cffee..354d923 100644
--- a/Source/energies.py
+++ b/Source/energies.py
@@ -1,9 +1,23 @@
+#!/usr/bin/env python
+"""
+Energies
+
+Compute energies (kinetic or total)
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-12-13
+"""
+
 import numpy as np
 
 def kinetic(W: np.ndarray) -> np.ndarray:
     """Computes the kinetic energy.
-    :param W: Phase-space vectors
-    :returns T: Kinetic energy
+    @param 
+        - W: Phase-space vectors
+    @returns 
+        - T: Kinetic energy
     """
     U = W[1,0]
     V = W[1,1]
diff --git a/Source/initial_conditions.py b/Source/initial_conditions.py
index 2ae0cd8..4f48950 100644
--- a/Source/initial_conditions.py
+++ b/Source/initial_conditions.py
@@ -1,3 +1,15 @@
+#!/usr/bin/env python
+"""
+Initial Conditions
+
+Generate initial conditions depending on different criteria
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-11-29
+"""
+
 import numpy as np
 import matplotlib.pyplot as plt
 POS_MIN = -1
@@ -12,12 +24,14 @@ def mesh_grid(N: int = N_PART,
               ymin: float = POS_MIN, 
               ymax: float = POS_MAX) -> np.ndarray:
     """Generates a set of regularly sampled particles with no velocity
-    :param N: number of particles
-    :parma xmin: Minimum value for position x
-    :param xmax: Maximum value for position x
-    :param ymin: Minimum value for position y
-    :param ymax: Maximum value for position y
-    :returns W: Phase-space vector
+    @ params:
+        - N: number of particles
+        - xmin: minimum value for position x
+        - xmax: maximum value for position x
+        - ymin: minimum value for position y
+        - ymax: maximum value for position y
+    @ returns:
+        - W: phase-space vector
     """
     X = np.linspace(xmin, xmax, N)
     Y = np.linspace(ymin, ymax, N)
@@ -30,10 +44,14 @@ def one_part(x0: float = 0,
              v0: float = 0) -> np.ndarray:
     """Generates a particle at position (x0, y0) 
     with velocities (u0,v0)
-    :param x0: Initial position x
-    :param y0: Initial position y
-    :param u0: Initial velocity u
-    :param v0: Initial velocity v
+    @ params:
+        - N: number of particles
+        - x0: initial position x
+        - y0: initial position y
+        - u0: initial velocity u
+        - v0: initial velocity v
+    @ returns:
+        - W: phase-space vector
     """
     X = x0
     Y = y0
@@ -48,6 +66,18 @@ def n_energy_part(potential,
                   xmax: float = +1,
                   ymin: float = -0.5,
                   ymax: float = +1):
+    """Generates N particles with an energy E in a potential.
+    @ params:
+        - potential: gravitational potential
+        - N: number of particles
+        - E: total energy
+        - xmin: minimum value for position x
+        - xmax: maximum value for position x
+        - ymin: minimum value for position y
+        - ymax: maximum value for position y
+    @ returns:
+        - W: an array of all the positions and velocities.
+    """
     X = []
     Y = []
     POT = []
@@ -81,6 +111,22 @@ def n_energy_2part(potential,
                    xmax: float = +1,
                    ymin: float = -0.5,
                    ymax: float = +1):
+    """Generate a sample of 2N particles with the energy E in a potential in 
+    two sets: one "normal" set (see n_energy_part), and a slightly shifted set
+    with a separation sep.
+    @ params:
+        - potential: gravitational potential
+        - N: number of particles
+        - E: total energy
+        - sep: the separation between the two sets
+        - xmin: minimum value for position x
+        - xmax: maximum value for position x
+        - ymin: minimum value for position y
+        - ymax: maximum value for position y
+    @ returns:
+        - (W1, W2): the two arrays of all the positions and velocities for 
+        each set.
+    """
     W_1 = n_energy_part(potential, N, E)
     W_2 = np.zeros_like(W_1)
     alpha = np.random.uniform(0, 2*np.pi, N)
diff --git a/Source/integrator.py b/Source/integrator.py
index 527d6da..ceb765a 100644
--- a/Source/integrator.py
+++ b/Source/integrator.py
@@ -1,7 +1,18 @@
+#!/usr/bin/env python
+"""
+Integrator
+
+Integrate differential equations.
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-11-29
+"""
 import numpy as np
 
 def euler(x0, y0, h, n, func): # FIXME cannot be used with vectors
-    """
+    """DEPRECIATED - DO NOT USE
     Euler method adapted for state vector [[x, y], [u, v]]
     :param x0: initial time value
     :param y0: initial state vector [[x, y], [u, v]]
@@ -25,7 +36,7 @@ def euler(x0, y0, h, n, func): # FIXME cannot be used with vectors
 
 # Updated RK2 integrator
 def rk2(x0, y0, h, n, func): # FIXME cannot be used with vectors
-    """
+    """ DEPRECIATED - DO NOT USE
     RK2 method adapted for state vector [[x, y], [u, v]]
     :param x0: initial time value
     :param y0: initial state vector [[x, y], [u, v]]
@@ -54,14 +65,15 @@ def rk4(t0: float,
         h: float, 
         n: int, 
         func):
-    """
-    RK4 method adapted for state vector [[x, y], [u, v]]
-    :param x0: initial time value
-    :param y0: initial state vector [[x, y], [u, v]]
-    :param h: step size (time step)
-    :param n: number of steps
-    :param func: RHS of differential equation
-    :returns: x array, solution array
+    """RK4 method adapted for state vector [[x, y], [u, v]]
+    @ params
+        - x0: initial time value
+        - y0: initial state vector [[x, y], [u, v]]
+        - h: step size (time step)
+        - n: number of steps
+        - func: RHS of differential equation
+    @returns: 
+        - t, W: time and state (solution) arrays, 
     """
     time = np.zeros(n)
     W = np.zeros((n,) + np.shape(W0))
@@ -84,4 +96,5 @@ def rk4(t0: float,
 
 
 def integrator_type(x0, y0, h, n, func, int_type):
+    """DEPRECIATED - DO NOT USE"""
     return int_type(x0, y0, h, n, func)
diff --git a/Source/plot_area.py b/Source/plot_area.py
new file mode 100644
index 0000000..e942745
--- /dev/null
+++ b/Source/plot_area.py
@@ -0,0 +1,78 @@
+#!/usr/bin/env python
+"""
+Plot: Area
+
+Plots areas
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-11-29
+"""
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+if "YII_1" in plt.style.available: plt.style.use("YII_1")
+
+OUT_DIR = "./Output/"
+FILENAME_PREFIX = "phase_separation_"
+EXTENSION = ".csv"
+
+def plot_area(filelist: list, mu_c = 1e-4) -> int:
+    """
+    Plot all the Poincaré sections in the file list.
+    @params:
+        - filelist: the list of files in the output directory, with the format 
+        "poincare_sections_{linear, parallel}_[1/E].csv"
+        - title: title of the figure
+    @returns: 
+        - 0.
+    """
+    orderlist = np.argsort([(int(file
+                                 .replace(FILENAME_PREFIX, "")
+                                 .replace(EXTENSION, "")))
+                            for file in filelist])
+    filelist = np.array(filelist)[orderlist]
+    N = len(filelist)
+    E =  np.linspace(1/100, 1/6, N)
+    mu = []
+    for filename in filelist:
+        with open(OUT_DIR + filename) as file:
+            data = file.readlines()
+            data = [np.float64(d.replace("\n", "")) for d in data]
+            mu.append(data)
+    mu = np.array(mu)
+
+    fig, ax = plt.subplots(1)
+    ax.scatter([], [], s=1, color="k", label="Data")
+    for i in range(len(mu)):
+        Y = mu[i]
+        ax.scatter([E[i]]*len(Y), Y, s=1, color="k", alpha=0.1)
+    ax.scatter(E, np.mean(mu, axis=1), s=5, 
+               color="C1", marker="o", label="Mean")
+    ax.scatter(E, np.median(mu, axis=1), s=5, 
+               color="C3", marker="s", label="Median")
+    ax.plot(E, [mu_c]*len(E), 
+            color="C5", label="Critical value $\\mu_\\mathrm{{c}}$")
+    ax.text(0.01, 1e-5, "Regular", va="bottom", ha="left", color="C5")
+    ax.text(0.01, 1e-3, "Chaotic", va="top", ha="left", color="C5")
+    ax.set_xlabel("Energy $E$")
+    ax.set_ylabel("Phase-space squared distance $\\mu$")
+    ax.set_yscale("log")
+    ax.legend()
+    fig.savefig("mu.pdf")
+
+    fig, ax = plt.subplots(1)
+    N_reg = np.count_nonzero(mu < mu_c, axis=1)
+    N = np.shape(mu)[1]
+    Area = N_reg / N
+    ax.scatter(E, Area, s=5, color="C0")
+    ax.set_xlabel("Energy $E$")
+    ax.set_ylabel("Area $N_\\mathrm{{reg}}/N_\\mathrm{{part}}$")
+    fig.savefig("area.pdf")
+    return 0
+
+filelist = [f for f in os.listdir(OUT_DIR) if FILENAME_PREFIX in f]
+plot_area(filelist)
+plt.show()
diff --git a/Source/plot_poincare_sections.py b/Source/plot_poincare_sections.py
index fb60c09..f74b9cc 100644
--- a/Source/plot_poincare_sections.py
+++ b/Source/plot_poincare_sections.py
@@ -61,6 +61,10 @@ def plot_poincare_sections(filelist: list, title:str = "") -> int:
     while i < N:
         i += 1
         axs[i].axis('off')
+    if "linear" in title: kind = "linear"
+    elif "parallel" in title: kind = "parallel"
+    else: kind = "error"
+    fig.savefig("pcs_{}.pdf".format(kind))
     return 0
 print("\033[32m" 
       + "[P]arallel or [L]inear algorithm result, or [B]oth?" 
diff --git a/Source/plot_poincare_sections_zoom.py b/Source/plot_poincare_sections_zoom.py
new file mode 100644
index 0000000..a338fd4
--- /dev/null
+++ b/Source/plot_poincare_sections_zoom.py
@@ -0,0 +1,101 @@
+#!/usr/bin/env python
+"""
+Plot: Poincaré Sections (Linear and Parallel)
+
+Plots the Poincaré sections for different energies, computed either with linear
+or parallel algorithms.
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-11-29
+"""
+import os
+import numpy as np
+import matplotlib.pyplot as plt
+
+if "YII_1" in plt.style.available: plt.style.use("YII_1")
+
+OUT_DIR = "./Output/"
+FILENAME_PREFIX = "poincare_sections_"
+EXTENSION = ".csv"
+
+def plot_poincare_sections(filelist: list, title:str = "") -> int:
+    """
+    Plot all the Poincaré sections in the file list.
+    @params:
+        - filelist: the list of files in the output directory, with the format 
+        "poincare_sections_{linear, parallel}_[1/E].csv"
+        - title: title of the figure
+    @returns: 
+        - 0.
+    """
+    orderlist = np.argsort([(int(file
+                                 .replace(FILENAME_PREFIX, "")
+                                 .replace(EXTENSION, "")
+                                 .replace("linear_", "")
+                                 .replace("parallel_", ""))) 
+                            for file in filelist])
+    filelist = np.array(filelist)[orderlist]
+    N = len(filelist)
+    #fig, axs = plt.subplot_mosaic("ABC\nDEF")
+    fig1, axs1 = plt.subplot_mosaic("AB")
+    fig2, axs2 = plt.subplot_mosaic("CD")
+    fig3, axs3 = plt.subplot_mosaic("EF")
+
+    axs = list(axs1.values()) + list(axs2.values()) + list(axs3.values())
+    fig1.suptitle(title)
+    fig2.suptitle(title)
+    fig3.suptitle(title)
+    for i in range(N):
+        ax = axs[i]
+        filename = filelist[i]
+        inv_E = (filename
+                 .replace(FILENAME_PREFIX, "")
+                 .replace(EXTENSION, "")
+                 .replace("linear_", "")
+                 .replace("parallel_", ""))
+        data = np.loadtxt(OUT_DIR + filename)
+        y_section = data[0]
+        v_section = data[1]
+        ax.scatter(y_section, v_section, 
+                   s=.1, color="C3", marker=",", alpha=0.5,
+                   label="$E = 1/{}$".format(inv_E))
+        ax.set_xlabel("$y$")
+        ax.set_ylabel("$v$")
+        ax.legend(loc="upper right")
+    while i < N:
+        i += 1
+        axs[i].axis('off')
+    if "linear" in title: kind = "linear"
+    elif "parallel" in title: kind = "parallel"
+    else: kind = "error"
+    fig1.savefig("pcs_zoom_1_{}.pdf".format(kind))
+    fig2.savefig("pcs_zoom_2_{}.pdf".format(kind))
+    fig3.savefig("pcs_zoom_3_{}.pdf".format(kind))
+    return 0
+
+print("\033[32m" 
+      + "[P]arallel or [L]inear algorithm result, or [B]oth?" 
+      + "\033[0m")
+answer = input("\033[32m" + "> " + "\033[0m").upper()
+
+if answer == "P":
+    FILENAME_PREFIX += "parallel_"
+elif answer == "L":
+    FILENAME_PREFIX += "linear_"
+
+filelist = [fname for fname in os.listdir(OUT_DIR) if FILENAME_PREFIX in fname]
+
+if answer in ["L", "B"]:
+    filelist_linear = [fname for fname in filelist if "linear_" in fname]
+    plot_poincare_sections(filelist_linear, 
+                           title=("Poincaré Sections "
+                                  "(results from the linear algorithm)"))
+if answer in ["P", "B"]:
+    filelist_parallel = [fname for fname in filelist if "parallel_" in fname]
+    plot_poincare_sections(filelist_parallel, 
+                           title=("Poincaré Sections "
+                                  "(results from the parallel algorithm)"))
+
+plt.show()
diff --git a/Source/poincare_sections.py b/Source/poincare_sections.py
index 85828f0..1ec7d07 100644
--- a/Source/poincare_sections.py
+++ b/Source/poincare_sections.py
@@ -1,5 +1,27 @@
+#!/usr/bin/env python
+"""
+Poincaré Sections
+
+Computes the Poincaré Sections
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-11-29
+"""
+
 import numpy as np
 def pcs_find(pos_x, pos_y, vel_x, vel_y):
+    """Find Poincaré sections (PCS; x = 0)
+    @ params:
+        - pos_x: position along the x axis
+        - pos_y: position along the y axis
+        - pos_x: velocity along the x axis
+        - pos_y: velocity along the y axis
+    @ returns: (tuple)
+        - pcs_pos_y: position of the points in the PCS along the y axis
+        - pcs_vel_y: velocity of the points in the PCS along the y axis
+    """
     if np.ndim(pos_x) == 1: 
         pos_x = np.array([pos_x])
         pos_y = np.array([pos_y])
@@ -25,8 +47,10 @@ def pcs_find(pos_x, pos_y, vel_x, vel_y):
                 pcs_vel_y.append(v0)
             i += 1
     return pcs_pos_y, pcs_vel_y
+
 def pcs_find_legacy(pos_x, pos_y, vel_x, vel_y):
-    """Depreciated legacy function that should not be used
+    """DEPRECIATED - DO NOT USE
+    Depreciated legacy function that should not be used
     """
     if np.ndim(pos_x) == 1: 
         pos_x = np.array([pos_x])
diff --git a/Source/potentials.py b/Source/potentials.py
index 282a9e4..d81d996 100644
--- a/Source/potentials.py
+++ b/Source/potentials.py
@@ -1,3 +1,15 @@
+#!/usr/bin/env python
+"""
+Potentials
+
+Functions of the different potentials (and their derivatives for the evolution)
+
+@ Author: Moussouni, Yaël (MSc student) & Bhat, Junaid Ramzan (MSc student)
+@ Institution:  Université de Strasbourg, CNRS, Observatoire astronomique
+                de Strasbourg, UMR 7550, F-67000 Strasbourg, France
+@ Date: 2024-11-29
+"""
+
 import numpy as np
 
 MAX_VAL = 1e3
@@ -11,7 +23,7 @@ def kepler_potential(W_grid: np.ndarray,
         - W: Phase-space vector
         - position_only: True if W is np.array([X, Y])
     @returns:
-        - V: computed potential
+        - computed potential
     """
     if position_only: 
         X = W_grid[0]
@@ -29,8 +41,11 @@ def kepler_potential(W_grid: np.ndarray,
 def hh_potential(W_grid: np.ndarray, 
                  position_only=False) -> np.ndarray:
     """Computes the Hénon-Heiles potential.
-    :param W: Phase-space vector
-    :output V: Potential
+    @params:
+        - W: Phase-space vector
+        - position_only: True if W is np.array([X, Y])
+    @returns:
+        - POT: Potential
     """
     if position_only:
         X = W_grid[0]
@@ -48,10 +63,12 @@ def hh_potential(W_grid: np.ndarray,
     return POT        
 
 def hh_evolution(t: np.ndarray, W: np.ndarray):
-    """Computes the evolution from the HH potential
-    :param t: Time (not used)
-    :param W: Phase space vector
-    :returns dot W: Time derivative of the phase space vector
+    """Computes the evolution from the HH potential derivative
+    @params
+        - t: Time (not used)
+        - W: Phase space vector
+    &returns 
+        - dot W: Time derivative of the phase space vector
     """
     X = W[0 ,0]
     Y = W[0, 1]
diff --git a/Source/test_evolution.py b/Source/test_evolution.py
index 33b88a3..08e9ab5 100644
--- a/Source/test_evolution.py
+++ b/Source/test_evolution.py
@@ -64,6 +64,8 @@ if __name__ == "__main__":
     ax.set_ylabel("$y$")
     ax.set_aspect("equal")
 
+    plt.savefig("evolution.png")
+
     plt.show(block=True)
 
 
diff --git a/Source/test_initial_E.py b/Source/test_initial_E.py
index 3c27ed5..9098e51 100644
--- a/Source/test_initial_E.py
+++ b/Source/test_initial_E.py
@@ -48,6 +48,8 @@ if __name__ == "__main__":
     ax.set_ylabel("$y$")
     ax.set_aspect("equal")
 
+    fig.savefig("initial_E.png")
+
     plt.show(block=True)
 
 
diff --git a/Source/test_potentials.py b/Source/test_potentials.py
index 84f923c..b205dbc 100644
--- a/Source/test_potentials.py
+++ b/Source/test_potentials.py
@@ -32,6 +32,8 @@ def kepler(W):
     ax.set_aspect("equal")
     ax.set_xlabel("$x$")
     ax.set_ylabel("$y$")
+
+    fig.savefig("pot_kepler.png")
     return 0
 
 def hh(W):
@@ -46,6 +48,8 @@ def hh(W):
     ax.set_aspect("equal")
     ax.set_xlabel("$x$")
     ax.set_ylabel("$y$")
+
+    fig.savefig("pot_hh.png")
     return 0
 
 if __name__ == "__main__":
diff --git a/area.pdf b/area.pdf
index 6cf7ee3..bc5f5e3 100644
Binary files a/area.pdf and b/area.pdf differ
diff --git a/evolution.png b/evolution.png
new file mode 100644
index 0000000..4656f87
Binary files /dev/null and b/evolution.png differ
diff --git a/initial_E.png b/initial_E.png
new file mode 100644
index 0000000..6f051ae
Binary files /dev/null and b/initial_E.png differ
diff --git a/mu.pdf b/mu.pdf
index f21bfab..2954596 100644
Binary files a/mu.pdf and b/mu.pdf differ
diff --git a/pcs_linear.pdf b/pcs_linear.pdf
new file mode 100644
index 0000000..3f97c33
Binary files /dev/null and b/pcs_linear.pdf differ
diff --git a/pcs_parallel.pdf b/pcs_parallel.pdf
new file mode 100644
index 0000000..b339961
Binary files /dev/null and b/pcs_parallel.pdf differ
diff --git a/pcs_zoom_1_linear.pdf b/pcs_zoom_1_linear.pdf
new file mode 100644
index 0000000..441393f
Binary files /dev/null and b/pcs_zoom_1_linear.pdf differ
diff --git a/pcs_zoom_1_parallel.pdf b/pcs_zoom_1_parallel.pdf
new file mode 100644
index 0000000..26c4068
Binary files /dev/null and b/pcs_zoom_1_parallel.pdf differ
diff --git a/pcs_zoom_2_linear.pdf b/pcs_zoom_2_linear.pdf
new file mode 100644
index 0000000..3fb9257
Binary files /dev/null and b/pcs_zoom_2_linear.pdf differ
diff --git a/pcs_zoom_2_parallel.pdf b/pcs_zoom_2_parallel.pdf
new file mode 100644
index 0000000..49cb064
Binary files /dev/null and b/pcs_zoom_2_parallel.pdf differ
diff --git a/pcs_zoom_3_linear.pdf b/pcs_zoom_3_linear.pdf
new file mode 100644
index 0000000..77e92dc
Binary files /dev/null and b/pcs_zoom_3_linear.pdf differ
diff --git a/pcs_zoom_3_parallel.pdf b/pcs_zoom_3_parallel.pdf
new file mode 100644
index 0000000..abce445
Binary files /dev/null and b/pcs_zoom_3_parallel.pdf differ
diff --git a/pot_hh.png b/pot_hh.png
new file mode 100644
index 0000000..59b7f73
Binary files /dev/null and b/pot_hh.png differ
diff --git a/pot_kepler.png b/pot_kepler.png
new file mode 100644
index 0000000..4dfbbfe
Binary files /dev/null and b/pot_kepler.png differ