update

Source/TODEL_area.py

@@ -0,0 +1,67 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import potentials as pot
+import integrator as itg
+import initial_conditions as init
+
+# Parameters
+h = 0.01  # Step size for integrator
+N_inter = 500 # Number of iterations for distance calculation
+eps = 1e-3  # Perturbation distance
+mu_c = 0.6  # Threshold for classifying ergodic vs regular
+energy_values = np.array([0.01, 0.02, 0.04, 0.08, 0.10, 0.12, 0.14, 0.16, 0.18])
+N_PART = 100
+
+def classify_point(P1, P1_prime):
+    pos_t1, positions1 = itg.rk4(0, P1, h, N_inter, pot.hh_evolution)
+    pos_t2, positions2 = itg.rk4(0, P1_prime, h, N_inter, pot.hh_evolution)
+
+    # Compute the sum of squared distances (μ)
+    mu = 0
+    for i in range(N_inter):
+        y1, y_dot1 = positions1[i, 0, 1], positions1[i, 1, 1]
+        y2, y_dot2 = positions2[i, 0, 1], positions2[i, 1, 1]
+        x1, x_dot1 = positions1[i, 0, 0], positions1[i, 1, 0]
+        x2, x_dot2 = positions2[i, 0, 0], positions2[i, 1, 0]
+
+        squared_distance = (np.sqrt((x2 - x1)**2 + (y2 - y1) ** 2 + (x_dot2 - x_dot1) ** 2 +(y_dot2 - y_dot1) ** 2))**2
+        mu += squared_distance
+
+    return mu
+
+# Analyze for different energy values
+relative_areas = []
+
+for E in energy_values:
+    # Generate initial conditions for the given energy level
+    initial_conditions = init.n_energy_part(pot.hh_potential, N=N_PART, E=E)
+    initial_conditions_eps = initial_conditions.copy()
+    initial_conditions_eps[0,1] += eps
+    n_regular = 0
+
+    # Loop through each initial condition
+    for i in range(N_PART):
+        P1 = initial_conditions[:, :, i]
+        # Perturb the initial point slightly
+        P1_prime = initial_conditions_eps[:,:,i]
+
+        # Classify the point based on the evolution
+        mu = classify_point(P1, P1_prime)
+        np.nan_to_num(mu, copy=False, nan=1.0)
+
+        print(f"the mu value is {mu}")
+        if mu < mu_c:
+            n_regular += 1
+
+    # Compute the relative area covered by regular points
+    relative_area = n_regular / N_PART
+    relative_areas.append(relative_area)
+
+# Plot the results
+plt.figure()
+plt.plot(energy_values, relative_areas, marker='o', color='C0')
+plt.xlabel("Energy (E)")
+plt.ylabel("Relative Area Covered by Regular Curves")
+plt.title("Transition from Ergodic to Stable Regime")
+plt.grid()
+plt.show()

Source/TODEL_equations.py

@@ -0,0 +1,17 @@
+from sympy import *
+
+# Forces
+X, Y, U, V = symbols("X, Y, U, V")
+POT = (X**2 + Y**2 + 2*X**2*Y - 2*Y**3/3)/2
+KIN = (U**2 + V**2)/2
+TOT = KIN + POT
+
+DX = diff(POT, X)
+DY = diff(POT, Y)
+
+FX = -DX
+FY = -DY
+
+print(FX)
+print(FY)
+

Source/TODO_chaos.py

@@ -0,0 +1,73 @@
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.integrate import quad
+
+import potentials as pot
+import integrator as itg
+import initial_conditions as init
+
+DEFAULT_N_iter = 30000
+DEFAULT_N_part = 100
+DEFAULT_h = 0.01
+DEFAULT_c = 1.7
+
+E_all = np.array([1/100, 1/12, 1/10, 1/8, 1/6])
+
+
+
+def compute_coordinates(E: float,
+                        N_iter: int = DEFAULT_N_iter,
+                        N_part: int = DEFAULT_N_part,
+                        h: float = DEFAULT_h) -> tuple:
+    W_part = init.n_energy_part(pot.hh_potential, N_part, E)
+
+    # Perform integration
+    t_part, coord_part = itg.rk4(0, W_part, h, N_iter, pot.hh_evolution)
+
+    # Extract positions and velocities
+    x_part = coord_part[:, 0, 0]
+    y_part = coord_part[:, 0, 1]
+    u_part = coord_part[:, 1, 0]
+    v_part = coord_part[:, 1, 1]
+
+    return t_part, x_part, y_part, u_part, v_part
+
+def x(t: float,
+      x_part: np.ndarray,
+      h: float = DEFAULT_h):
+    return x_part[t//h]
+
+def theta(t: float, 
+          x,
+          x_part: np.ndarray,
+          c: float = DEFAULT_c, 
+          h: float = DEFAULT_h,
+          ds: float = DEFAULT_h):
+    """theta(t) = c*t + int_0^t phi(x(s)) ds"""
+    I = quad(lambda s: x(s, x_part), 0, t)
+    return c*t + I
+
+def p(t: float,
+      x,
+      theta,
+      x_part: np.ndarray,
+      h: float = DEFAULT_h,
+      ds: float = DEFAULT_h):
+    I = quad(lambda s: x(s, x_part)*cos(theta(s, x, x_part)), 0, t)
+    return I
+
+def M(t,
+      p,
+      x_part: np.ndarray,
+      h: float = DEFAULT_h,
+      dtau: float = DEFAULT_h):
+    I = quad(lambda tau: (p(t+tau, x, theta, x_part) 
+                          - p(tau, x, theta, x_part)))**2, 0, np.max(t))
+    return I
+
+if __name__ == "__main__":
+    for i in range(len(E_all)):
+        E = E_all[i]
+        t_part, x_part, y_part, v_part, u_part = compute_coordinates(E)
+        I = M(t_part, x_part)
+

Source/YII_1.mplstyle

@@ -0,0 +1,78 @@
+### Axes
+axes.titlesize : large
+axes.labelsize : large
+
+### Lines
+lines.linewidth : 1.0
+hatch.linewidth: 0.25
+lines.markersize: 5
+
+patch.antialiased   : True
+
+### Ticks
+xtick.top: True
+xtick.bottom: True
+xtick.major.size: 4.0
+xtick.minor.size: 2.0
+xtick.major.width: 0.5
+xtick.minor.width: 0.5
+xtick.direction: in
+xtick.minor.visible: True
+xtick.major.top: True
+xtick.major.bottom: True
+xtick.minor.top: True
+xtick.minor.bottom: True
+xtick.major.pad: 5.0
+xtick.minor.pad: 5.0
+ytick.labelsize: large
+
+ytick.left: True
+ytick.right: True
+ytick.major.size: 4.0
+ytick.minor.size: 2.0
+ytick.major.width: 0.5
+ytick.minor.width: 0.5
+ytick.direction: in
+ytick.minor.visible: True
+ytick.major.left: True
+ytick.major.right: True
+ytick.minor.left: True
+ytick.minor.right: True
+ytick.major.pad: 5.0
+ytick.minor.pad: 5.0
+xtick.labelsize: large
+
+### Legend
+legend.frameon: True
+legend.fontsize: 9.0
+
+### Figure size
+figure.figsize: 5.9, 3.9
+figure.subplot.left: 0.1
+figure.subplot.bottom: 0.175
+figure.subplot.top: 0.90
+figure.subplot.right: 0.95
+figure.autolayout: True
+figure.dpi: 150
+
+### Fonts
+font.family: serif
+font.serif: Latin Modern Roman
+font.size: 9.0
+text.usetex: true
+
+### Saving figures
+savefig.format: pdf
+savefig.dpi: 300
+savefig.pad_inches: 0
+path.simplify: True
+
+### Cycler
+axes.prop_cycle: cycler(color=['06518E','E8BD0F','3B795E','ED1C24','C95BCF','009EDB'])
+
+### Colormap
+image.cmap: cividis
+
+
+### LaTeX 
+text.latex.preamble: \input{~/.config/matplotlib/stylelib/yii_plot.tex}

Source/YII_light_1.mplstyle

@@ -0,0 +1,74 @@
+### Axes
+axes.titlesize : large
+axes.labelsize : large
+
+### Lines
+lines.linewidth : 1.0
+hatch.linewidth: 0.25
+lines.markersize: 5
+
+patch.antialiased   : True
+
+### Ticks
+xtick.top: True
+xtick.bottom: True
+xtick.major.size: 4.0
+xtick.minor.size: 2.0
+xtick.major.width: 0.5
+xtick.minor.width: 0.5
+xtick.direction: in
+xtick.minor.visible: False
+xtick.major.top: True
+xtick.major.bottom: True
+xtick.minor.top: False
+xtick.minor.bottom: False
+xtick.major.pad: 5.0
+xtick.minor.pad: 5.0
+ytick.labelsize: large
+
+ytick.left: True
+ytick.right: True
+ytick.major.size: 4.0
+ytick.minor.size: 2.0
+ytick.major.width: 0.5
+ytick.minor.width: 0.5
+ytick.direction: in
+ytick.minor.visible: False
+ytick.major.left: True
+ytick.major.right: True
+ytick.minor.left: False
+ytick.minor.right: False
+ytick.major.pad: 5.0
+ytick.minor.pad: 5.0
+xtick.labelsize: large
+
+### Legend
+legend.frameon: True
+legend.fontsize: 9.0
+
+### Figure size
+figure.figsize: 5.9, 3.9
+figure.subplot.left: 0.1
+figure.subplot.bottom: 0.175
+figure.subplot.top: 0.90
+figure.subplot.right: 0.95
+figure.autolayout: True
+figure.dpi: 150
+
+### Fonts
+font.family: monospace
+font.serif: Latin Modern Mono
+font.size: 9.0
+#text.usetex: true
+
+### Saving figures
+savefig.format: pdf
+savefig.dpi: 300
+savefig.pad_inches: 0
+path.simplify: True
+
+### Cycler
+axes.prop_cycle: cycler(color=['06518E','E8BD0F','3B795E','ED1C24','C95BCF','009EDB'])
+
+### Colormap
+image.cmap: cividis

Source/energies.py

@@ -0,0 +1,20 @@
+import numpy as np
+
+def kinetic(W: np.ndarray) -> np.ndarray:
+    """Computes the kinetic energy.
+    :param W: Phase-space vectors
+    :returns T: Kinetic energy
+    """
+    U = W[1,0]
+    V = W[1,1]
+    # If U or V is not an array (or a list), but rather a scalar, then we
+    # create a list of one element so that it can work either way
+    if np.ndim(U) == 0: U = np.array([U])
+    if np.ndim(V) == 0: V = np.array([V])
+    
+    return (U**2 + V**2)/2
+
+def total(W: np.ndarray,
+          potential,
+          kinetic = kinetic) -> np.ndarray:
+    return potential(W) + kinetic(W)

Source/initial_conditions.py

@@ -0,0 +1,96 @@
+import numpy as np
+import matplotlib.pyplot as plt
+POS_MIN = -1
+POS_MAX = +1
+VEL_MIN = -1
+VEL_MAX = +1
+N_PART = 1
+
+def mesh_grid(N: int = N_PART,
+              xmin: float = POS_MIN, 
+              xmax: float = POS_MAX, 
+              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
+    """
+    X = np.linspace(xmin, xmax, N)
+    Y = np.linspace(ymin, ymax, N)
+    X,Y = np.meshgrid(X,Y, indexing="ij")
+    return np.array([X,Y])
+
+def one_part(x0: float = 0, 
+             y0: float = 0, 
+             u0: 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
+    """
+    X = x0
+    Y = y0
+    U = u0
+    V = v0
+    return np.array([[X,Y], [U,V]])
+
+def n_energy_part(potential,
+                  N: int = N_PART,
+                  E: float = 0,
+                  xmin: float = -1,
+                  xmax: float = +1,
+                  ymin: float = -0.5,
+                  ymax: float = +1):
+    X = []
+    Y = []
+    POT = []
+    U = []
+    V = []
+    while len(X) < N:
+        x = np.random.random()*(xmax-xmin)+xmin
+        y = np.random.random()*(ymax-ymin)+ymin
+        w = np.array([x, y])
+        pot = potential(w, position_only=True)[0]
+        if pot <= E:
+            X.append(x)
+            Y.append(y)
+            POT.append(pot)
+    X = np.array(X)
+    Y = np.array(Y)
+    POT = np.array(POT)
+    U = np.zeros_like(X)
+    V = np.zeros_like(Y)
+    C = np.sqrt(2 * (E - POT))
+    THETA = np.random.random(N)*2*np.pi
+    U = C*np.cos(THETA)
+    V = C*np.sin(THETA)
+    return np.array([[X, Y], [U, V]])
+
+def n_energy_2part(potential,
+                   N: int = N_PART,
+                   E: float = 0,
+                   sep: float = 1e-7,
+                   xmin: float = -1,
+                   xmax: float = +1,
+                   ymin: float = -0.5,
+                   ymax: float = +1):
+    W_1 = n_energy_part(potential, N, E)
+    W_2 = np.zeros_like(W_1)
+    alpha = np.random.uniform(0, 2*np.pi, N)
+    W_2[0, 0] = W_1[0, 0] + sep*np.cos(alpha)
+    W_2[0, 1] = W_1[0, 1] + sep*np.sin(alpha)
+    W_2[1, 0] = W_1[1, 0]
+    W_2[1, 1] = W_1[1, 1]
+    return (W_1, W_2)
+
+
+
+
+

Source/integrator.py

@@ -0,0 +1,87 @@
+import numpy as np
+
+def euler(x0, y0, h, n, func): # FIXME cannot be used with vectors
+    """
+    Euler 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
+    """
+    x_values = np.zeros(n)
+    y_values = np.zeros((n, 2, 2))  #  to accommodate the state vector
+
+    for i in range(n):
+        dydt = func(x0, y0)
+        y0 = y0 + h * dydt
+        x0 = x0 + h
+
+        x_values[i] = x0
+        y_values[i, :, :] = y0
+
+    return x_values, y_values
+
+# Updated RK2 integrator
+def rk2(x0, y0, h, n, func): # FIXME cannot be used with vectors
+    """
+    RK2 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
+    """
+    x_values = np.zeros(n)
+    y_values = np.zeros((n, 2, 2))  #  to accommodate the state vector
+
+    for i in range(n):
+        k1 = func(x0, y0)
+        k2 = func(x0 + h / 2., y0 + h / 2. * k1)
+
+        y0 = y0 + h * (k1 / 2. + k2 / 2.)
+        x0 = x0 + h
+
+        x_values[i] = x0
+        y_values[i, :, :] = y0
+
+    return x_values, y_values
+
+def rk4(t0: float, 
+        W0: np.ndarray, 
+        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
+    """
+    time = np.zeros(n)
+    W = np.zeros((n,) + np.shape(W0))
+    # to accommodate the state vector
+    t = t0
+    w = W0
+    for i in range(n):
+        k1 = func(t, w)
+        k2 = func(t + h / 2., w + h / 2. * k1)
+        k3 = func(t + h / 2., w + h / 2. * k2)
+        k4 = func(t + h, w + h * k3)
+
+        w = w + h * (k1 / 6. + k2 / 3. + k3 / 3. + k4 / 6.)
+        t = t + h
+
+        time[i] = t
+        W[i] = w
+
+    return t, W
+
+
+def integrator_type(x0, y0, h, n, func, int_type):
+    return int_type(x0, y0, h, n, func)

Source/main_area.py

@@ -0,0 +1,70 @@
+#!/usr/bin/env python
+"""
+Main: Compute Relative Area
+
+Computes the relative area covered bu the curves for different energies, to 
+study ordered and chaotic regimes.
+
+@ 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
+from scipy.optimize import curve_fit
+
+import potentials as pot
+import energies as ene
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+OUT_DIR = "./Output/"
+FILENAME_PREFIX = "phase_separation_"
+EXTENSION = ".csv"
+DEFAULT_N_iter = int(1e5)
+DEFAULT_N_part = 200
+DEFAULT_h = 0.005
+E_all = np.linspace(1/100, 1/6, 20)
+
+def compute_mu(E: float,
+               N_iter: int = DEFAULT_N_iter,
+               N_part: int = DEFAULT_N_part,
+               h: float = DEFAULT_h) -> tuple:
+    """
+    Computes the phase-space squared distances for particles of given energy E.
+    @params:
+        - E: the total energy of each particles
+        - N_iter: the number of iteration
+        - N_part: the number of particles
+        - h: integration steps
+    @returns:
+        - mu: phase-space squared distance
+    """
+    W_1, W_2 = init.n_energy_2part(pot.hh_potential, N_part, E)
+    t_1, positions_1 = itg.rk4(0, W_1, h, N_iter, pot.hh_evolution)
+    x_1 = positions_1[:, 0, 0]
+    y_1 = positions_1[:, 0, 1]
+    u_1 = positions_1[:, 1, 0]
+    v_1 = positions_1[:, 1, 1]
+    
+    t_2, positions_2 = itg.rk4(0, W_2, h, N_iter, pot.hh_evolution)
+    x_2 = positions_2[:, 0, 0]
+    y_2 = positions_2[:, 0, 1]
+    u_2 = positions_2[:, 1, 0]
+    v_2 = positions_2[:, 1, 1]
+    dist_sq = (x_2[-25:] - x_1[-25:])**2 \
+            + (y_2[-25:] - y_1[-25:])**2 \
+            + (u_2[-25:] - u_1[-25:])**2 \
+            + (v_2[-25:] - v_1[-25:])**2
+    
+    mu = np.sum(dist_sq, axis=0)
+    return mu
+
+if __name__ == "__main__":
+    mu_all = []
+    for i in range(len(E_all)):
+        mu = compute_mu(E_all[i])
+        filename = OUT_DIR + FILENAME_PREFIX\
+                 + str(i) + EXTENSION
+        np.savetxt(filename, mu)

Source/main_poincare_sections_linear.py

@@ -0,0 +1,86 @@
+#!/usr/bin/env python
+"""
+Main: Computes Poincaré Sections (Linear Algorithm)
+
+Computes the Poincaré Sections with a linear algorithm 
+(i.e. no parallel computing).
+
+@ 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 potentials as pot
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+# Parameters
+OUT_DIR = "./Output/"
+FILENAME_PREFIX = "poincare_sections_linear_"
+EXTENSION = ".csv"
+DEFAULT_N_iter = 30000
+DEFAULT_N_part = 100
+DEFAULT_h = 0.01
+E_all = np.array([1/100, 1/12, 1/10, 1/8, 1/6])
+
+text_E = ["1/100", "1/12", "1/10", "1/8", "1/6"]
+
+def compute_poincare_sections_linear(E: float,
+                                    N_iter: int = DEFAULT_N_iter,
+                                    N_part: int = DEFAULT_N_part,
+                                    h: float = DEFAULT_h) -> tuple:
+    """
+    Computes the Poincaré sections for a given energy E.
+    @params:
+        - E: the total energy of each particles
+        - N_iter: the number of iteration
+        - N_part: the number of particles
+        - h: integration steps
+    @returns:
+        - y_section, v_section: arrays containing the y and v coordinates of 
+          the Poincaré sections
+    """
+    W_all_part = init.n_energy_part(pot.hh_potential, N_part, E)
+    y_section = []
+    v_section = []
+    for i in range(N_part):
+        W_part = W_all_part[:,:,i]
+        
+        # Perform integration
+        t_part, coord_part = itg.rk4(0, W_part, h, N_inter, pot.hh_evolution)
+
+        # Extract positions and velocities
+        x_part = coord_part[:, 0, 0]
+        y_part = coord_part[:, 0, 1]
+        u_part = coord_part[:, 1, 0]
+        v_part = coord_part[:, 1, 1]
+
+        # Find Poincaré section points for the current initial condition
+        y_pcs, v_pcs = pcs.pcs_find_legacy(x_part, y_part, u_part, v_part)
+        # The legacy is important here, the algorithm is the same but the 
+        # data format is different...
+
+        # Append the current Poincaré section points to the overall lists
+        y_section += y_pcs
+        v_section += v_pcs
+    return y_section, v_section
+
+if __name__ == "__main__":
+    y_section_all = []
+    v_section_all = []
+    for i in range(len(E_all)):
+        y_section, v_section = compute_poincare_sections_linear(E_all[i])
+        section = np.array([y_section, v_section])
+        filename = OUT_DIR + FILENAME_PREFIX\
+                + str(text_E[i][2:]) + EXTENSION
+        np.savetxt(filename, section)
+
+
+
+
+
+    
+

Source/main_poincare_sections_parallel.py

@@ -0,0 +1,71 @@
+#!/usr/bin/env python
+"""
+Main: Computes Poincaré Sections (Parallel Algorithm)
+
+Computes the Poincaré Sections with a parallel algorithm (with Numpy)
+
+@ 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 potentials as pot
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+# Parameters
+OUT_DIR = "./Output/"
+FILENAME_PREFIX = "poincare_sections_parallel_"
+EXTENSION = ".csv"
+DEFAULT_N_iter = 30000
+DEFAULT_N_part = 100
+DEFAULT_h = 0.01
+E_all = np.array([1/100, 1/12, 1/10, 1/8, 1/6])
+
+text_E = ["1/100", "1/12", "1/10", "1/8", "1/6"]
+
+def compute_poincare_sections_numpy(E: float,
+                                    N_iter: int = DEFAULT_N_iter,
+                                    N_part: int = DEFAULT_N_part,
+                                    h: float = DEFAULT_h) -> tuple:
+    """
+    Computes the Poincaré sections for a given energy E.
+    @params:
+        - E: the total energy of each particles
+        - N_iter: the number of iteration
+        - N_part: the number of particles
+        - h: integration steps
+    @returns:
+        - y_section, v_section: arrays containing the y and v coordinates of 
+          the Poincaré sections
+    """
+    W_part = init.n_energy_part(pot.hh_potential, N_part, E)
+    y_section = []
+    v_section = []
+    
+    # Perform integration
+    t_part, coord_part = itg.rk4(0, W_part, h, N_iter, pot.hh_evolution)
+
+    # Extract positions and velocities
+    x_part = coord_part[:, 0, 0]
+    y_part = coord_part[:, 0, 1]
+    u_part = coord_part[:, 1, 0]
+    v_part = coord_part[:, 1, 1]
+
+    # Find Poincaré section points for the current initial condition
+    y_section, v_section = pcs.pcs_find(x_part, y_part, u_part, v_part)
+    return y_section, v_section
+
+if __name__ == "__main__":
+    y_section_all = []
+    v_section_all = []
+    for i in range(len(E_all)):
+        y_section, v_section = compute_poincare_sections_numpy(E_all[i])
+        section = np.array([y_section, v_section])
+        filename = OUT_DIR + FILENAME_PREFIX\
+                 + str(text_E[i][2:]) + EXTENSION
+        np.savetxt(filename, section)
+

Source/main_yael.py

@@ -0,0 +1,114 @@
+import time
+import numpy as np
+import matplotlib.pyplot as plt
+from scipy.optimize import curve_fit
+
+import potentials as pot
+import energies as ene
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+if "YII_light_1" in plt.style.available: 
+    plt.style.use("YII_light_1")
+# Loading matplotlib style...
+# PARAM
+N_grid = 1000
+N_inter = int(1e5)
+N_part = 200
+#all_E = np.array([1/100, 1/12, 1/10, 1/8, 1/6])
+all_E = np.linspace(1/100, 1/6, 20)
+h = 0.005
+
+mu_c = 1e-4
+d_12 = 1e-7
+
+# INIT
+W_grid = init.mesh_grid(N_grid, xmin=-1, xmax=1, ymin=-1, ymax=1)
+X_grid = W_grid[0]
+Y_grid = W_grid[1]
+potential = pot.hh_potential(W_grid, position_only=True)
+
+# MAIN
+def gen_mu(E):
+    pot_valid = np.ma.masked_where(potential > E, potential)
+    W_1 = init.n_energy_part(pot.hh_potential, N_part, E)
+    W_2 = np.zeros_like(W_1)
+    alpha = np.random.uniform(0, 2*np.pi, N_part)
+    W_2[0, 0] = W_1[0, 0] + d_12*np.cos(alpha)
+    W_2[0, 1] = W_1[0, 1] + d_12*np.sin(alpha)
+    W_2[1, 0] = W_1[1, 0]
+    W_2[1, 1] = W_1[1, 1]
+
+    t_1, positions_1 = itg.rk4(0, W_1, h, N_inter, pot.hh_evolution)
+    x_1 = positions_1[:, 0, 0]
+    y_1 = positions_1[:, 0, 1]
+    u_1 = positions_1[:, 1, 0]
+    v_1 = positions_1[:, 1, 1]
+    
+    t_2, positions_2 = itg.rk4(0, W_2, h, N_inter, pot.hh_evolution)
+    x_2 = positions_2[:, 0, 0]
+    y_2 = positions_2[:, 0, 1]
+    u_2 = positions_2[:, 1, 0]
+    v_2 = positions_2[:, 1, 1]
+    dist_sq = (x_2[-25:] - x_1[-25:])**2 \
+            + (y_2[-25:] - y_1[-25:])**2 \
+            + (u_2[-25:] - u_1[-25:])**2 \
+            + (v_2[-25:] - v_1[-25:])**2
+
+    mu = np.sum(dist_sq, axis=0)
+    return mu
+
+A = np.zeros_like(all_E)
+MU = []
+ALL_E = []
+for i in range(len(all_E)):
+    mu = gen_mu(all_E[i])
+    n_total = N_part
+    n_ergo = np.sum(mu < mu_c) # count how many times the condition is verified
+    A[i] = n_ergo/n_total
+    MU.append(mu)
+    ALL_E.append([all_E[i]]*len(mu))
+
+MU = np.array(MU)
+ALL_E = np.array(ALL_E)
+
+# def lin(x, a, b):
+#     return a*x + b
+
+# w_chaos = np.argwhere(A < 1).flatten()
+# X = all_E[w_chaos]
+# Y = A[w_chaos]
+
+# popt, pcov = curve_fit(lin, X, Y)
+# a, b = popt
+# da, db = np.sqrt(np.diag(pcov))
+
+fig, ax = plt.subplots(1)
+ax.scatter(ALL_E, MU, s=1, 
+           color="k", alpha=0.1, label="Data")
+ax.scatter(all_E, np.mean(MU, axis=1), s=10, 
+           color="C3", marker="*", label="Mean")
+ax.scatter(all_E, np.median(MU, axis=1), s=10, 
+           color="C2", marker="s", label="Median")
+ax.hlines(mu_c, np.min(all_E), np.max(all_E), color="C4")
+ax.legend()
+ax.set_yscale("log")
+ax.set_xlabel("energy E")
+ax.set_ylabel("quantity mu")
+
+fig.tight_layout()
+fig.savefig("mu")
+
+fig, ax = plt.subplots(1)
+ax.scatter(all_E, A, color="C0", s=5, marker="o", label="Data")
+ax.set_xlabel("energy E")
+ax.set_ylabel("area A")
+ax.legend()
+fig.tight_layout()
+fig.savefig("area")
+plt.show()
+
+
+
+    

Source/plot_poincare_sections.py

@@ -0,0 +1,88 @@
+#!/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")
+    axs = list(axs.values())
+    fig.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')
+    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()

Source/poincare_sections.py

@@ -0,0 +1,55 @@
+import numpy as np
+def pcs_find(pos_x, pos_y, vel_x, vel_y):
+    if np.ndim(pos_x) == 1: 
+        pos_x = np.array([pos_x])
+        pos_y = np.array([pos_y])
+        vel_x = np.array([vel_x])
+        vel_y = np.array([vel_y])
+    pcs_pos_y = []
+    pcs_vel_y = []
+    for j in range(len(pos_x[0])): # for each particle
+        i = 0
+        x = pos_x[:,j]
+        y = pos_y[:,j]
+        u = vel_x[:,j]
+        v = vel_y[:,j]
+        while i < len(x) - 1:
+            if x[i] * x[i+1] < 0:
+                y0 = y[i] \
+                        + (y[i+1] - y[i])/(x[i+1] - x[i]) \
+                        * (0 - x[i])
+                v0 = v[i] \
+                        + (v[i+1] - v[i])/(x[i+1] - x[i])\
+                        * (0 - x[i])
+                pcs_pos_y.append(y0)
+                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
+    """
+    if np.ndim(pos_x) == 1: 
+        pos_x = np.array([pos_x])
+        pos_y = np.array([pos_y])
+        vel_x = np.array([vel_x])
+        vel_y = np.array([vel_y])
+    pcs_pos_y = []
+    pcs_vel_y = []
+    for j in range(len(pos_x)): # for each particle
+        i = 0
+        x = pos_x[j]
+        y = pos_y[j]
+        u = vel_x[j]
+        v = vel_y[j]
+        while i < len(x) - 1:
+            if x[i] * x[i+1] < 0:
+                y0 = y[i] \
+                        + (y[i+1] - y[i])/(x[i+1] - x[i]) \
+                        * (0 - x[i])
+                v0 = v[i] \
+                        + (v[i+1] - v[i])/(x[i+1] - x[i])\
+                        * (0 - x[i])
+                pcs_pos_y.append(y0)
+                pcs_vel_y.append(v0)
+            i += 1
+    return pcs_pos_y, pcs_vel_y

Source/poincare_test.py

@@ -0,0 +1,58 @@
+import numpy as np
+import matplotlib.pyplot as plt
+import potentials as pot
+import energies as ene
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+h = 0.01
+N_inter = 30000
+eps = 1e-2
+
+x_0 = 0.0
+y_vals = np.linspace(0, 0.6, 50)
+E = 0.125  # 1/12
+
+# Initialize lists to store all Poincaré points across different initial conditions
+all_pcs_pos_y = []
+all_pcs_vel_y = []
+
+
+for i in range(len(y_vals)):
+    y0 = y_vals[i]
+    W_part = init.one_part(x_0, y0, 0, 0)
+    POT = pot.hh_potential(W_part)[0]
+    u_0 = np.sqrt(2 * (E - POT))
+    v_0 = 0.0
+
+
+    W_part[1, 0] = u_0
+    W_part[1, 1] = v_0
+
+    # Perform integration
+    pos_t, positions = itg.rk4(0, W_part, h, N_inter, pot.hh_evolution)
+
+    # Extract positions and velocities
+    pos_x = positions[:, 0, 0]
+    pos_y = positions[:, 0, 1]
+    vel_x = positions[:, 1, 0]
+    vel_y = positions[:, 1, 1]
+
+    # Find Poincaré section points for the current initial condition
+    pcs_pos_y, pcs_vel_y = pcs.pcs_find(pos_x, pos_y, vel_x, vel_y)
+
+    # Append the current Poincaré section points to the overall lists
+    all_pcs_pos_y.extend(pcs_pos_y)
+    all_pcs_vel_y.extend(pcs_vel_y)
+
+# Plotting all Poincaré section points on the same graph
+fig, ax = plt.subplots()
+ax.scatter(all_pcs_pos_y, all_pcs_vel_y, s=2, color="C4")
+
+# Set labels and legend
+ax.set_xlabel("coordinate y")
+ax.set_ylabel("velocity v")
+ax.set_title("Combined Poincaré Section for Multiple Initial y0 Values")
+
+plt.show()

Source/potentials.py

@@ -0,0 +1,65 @@
+import numpy as np
+
+MAX_VAL = 1e3
+
+def kepler_potential(W_grid: np.ndarray, 
+                     position_only: bool = False) -> np.ndarray:
+    """Computes the Kepler potential: V(R) = -G*m1*m2/R 
+    (assuming G = 1, m1 = 1, m2 = 1)
+    assuming the point mass at (x = 0, y = 0).
+    @params:
+        - W: Phase-space vector
+        - position_only: True if W is np.array([X, Y])
+    @returns:
+        - V: computed potential
+    """
+    if position_only: 
+        X = W_grid[0]
+        Y = W_grid[1]
+    else:
+        X = W_grid[0,0]
+        Y = W_grid[0,1]
+    # If X or Y is not an array (or a list), but rather a scalar, then we
+    # create a list of one element so that it can work either way
+    if np.ndim(X) == 0: X = np.array([X])
+    if np.ndim(Y) == 0: Y = np.array([Y])
+    R = np.sqrt(X**2 + Y**2)
+    return -1/R
+
+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
+    """
+    if position_only:
+        X = W_grid[0]
+        Y = W_grid[1]
+    else:
+        X = W_grid[0, 0]
+        Y = W_grid[0, 1]
+        
+    # If X or Y is not an array (or a list), but rather a scalar, then we
+    # create a list of one element so that it can work either way
+    if np.ndim(X) == 0: X = np.array([X])
+    if np.ndim(Y) == 0: Y = np.array([Y])
+
+    POT = (X**2 + Y**2 + 2*X**2*Y - 2*Y**3/3)/2
+    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
+    """
+    X = W[0 ,0]
+    Y = W[0, 1]
+    U = W[1, 0]
+    V = W[1, 1]
+    DX = U 
+    DY = V
+    DU = -(2*X*Y + X)
+    DV = -(X**2 - Y**2 + Y)
+
+    return np.array([[DX, DY], [DU, DV]])

Source/test_evolution.py

@@ -0,0 +1,69 @@
+#!/usr/bin/env python
+"""
+Test: Evolution
+
+Evolve a particle with a given energy in a potential and show the result path
+followed by the particle in a given time.
+
+@ 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
+import potentials as pot
+import energies as ene
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+if "YII_1" in plt.style.available: plt.style.use("YII_1")
+# Loading matplotlib style...
+
+# Constants
+N_grid = 5000
+h = 0.01
+N_inter = 5000
+eps = 1e-2
+
+x_0 = 0.0
+y_0 = 0.5
+E = 0.125
+
+# Main
+if __name__ == "__main__":
+    W_part = init.one_part(x_0, y_0,0,0)
+    POT = pot.hh_potential(W_part)[0]
+    u_0 = np.sqrt(2 * (E - POT))
+    v_0 = 0.0
+
+    W_part[1,0] = u_0
+    W_part[1,1] = v_0
+
+    pos_t, positions = itg.rk4(0, W_part, h, N_inter, pot.hh_evolution)
+
+    pos_x = positions[:,0,0]
+    pos_y = positions[:,0,1]
+    vel_x = positions[:,1,0]
+    vel_y = positions[:,1,1]
+
+    W_grid = init.mesh_grid(N_grid)
+    X_grid = W_grid[0]
+    Y_grid = W_grid[1]
+    potential = pot.hh_potential(W_grid, position_only=True)
+
+    fig, ax = plt.subplots(1)
+
+    pcm = ax.pcolormesh(X_grid, Y_grid, potential)
+    line = ax.plot(pos_x, pos_y, color="C3")
+    fig.colorbar(pcm, label="potential")
+
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$y$")
+    ax.set_aspect("equal")
+
+    plt.show(block=True)
+
+

Source/test_initial_E.py

@@ -0,0 +1,53 @@
+#!/usr/bin/env python
+"""
+Test: Initial Conditions With a Given Energy
+
+Sample random particles with the same given energy in a valid coordinates range.
+
+@ 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
+
+import potentials as pot
+import energies as ene
+import integrator as itg
+import initial_conditions as init
+import poincare_sections as pcs
+
+if "YII_1" in plt.style.available: plt.style.use("YII_1")
+# Loading matplotlib style...
+# parameters
+N_grid = 1000
+N_inter = 30000
+N_part = 100
+E = 1/12
+h = 0.01
+
+if __name__ == "__main__":
+    # Initial conditions
+    W_grid = init.mesh_grid(N_grid, xmin=-1, xmax=1, ymin=-1, ymax=1)
+    X_grid = W_grid[0]
+    Y_grid = W_grid[1]
+    potential = pot.hh_potential(W_grid, position_only=True)
+    pot_valid = np.ma.masked_where(potential > E, potential)
+
+    W_all_part = init.n_energy_part(pot.hh_potential, N_part, E)
+
+    # Plot
+    fig, ax = plt.subplots(1)
+
+    pcm = ax.pcolormesh(X_grid, Y_grid, pot_valid, vmin = 0)
+    sct = ax.scatter(W_all_part[0, 0], W_all_part[0, 1], s=1, color="C3")
+    fig.colorbar(pcm, label="potential")
+    ax.set_title("$E = {:.2f}$".format(E))
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$y$")
+    ax.set_aspect("equal")
+
+    plt.show(block=True)
+
+

Source/test_potentials.py

@@ -0,0 +1,56 @@
+#!/usr/bin/env python
+"""
+Test: Potential
+
+Draw the Kepler potential and the Hénon--Heils potential
+
+@ 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
+
+import potentials as pot
+import initial_conditions as init
+
+if "YII_1" in plt.style.available: plt.style.use("YII_1")
+# Loading matplotlib style...
+
+W = init.mesh_grid(300)
+
+def kepler(W):
+    """Plots the Kepler potential"""
+    X = W[0]
+    Y = W[1]
+    POT = pot.kepler_potential(W, position_only=True)
+    fig, ax = plt.subplots(1)
+    ax.set_title("Kepler potential")
+    pcm = ax.pcolormesh(X, Y, POT)
+    fig.colorbar(pcm, label="potential")
+    ax.set_aspect("equal")
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$y$")
+    return 0
+
+def hh(W):
+    """Plots the Hénon--Heils potential"""
+    X = W[0]
+    Y = W[1]
+    POT = pot.hh_potential(W, position_only=True)
+    fig, ax = plt.subplots(1)
+    ax.set_title("Hénon–Heils potential")
+    pcm = ax.pcolormesh(X, Y, POT)
+    fig.colorbar(pcm, label="potential")
+    ax.set_aspect("equal")
+    ax.set_xlabel("$x$")
+    ax.set_ylabel("$y$")
+    return 0
+
+if __name__ == "__main__":
+    kepler(W)
+    hh(W)
+
+    plt.show(block=True)
+

Source/time_poincare_sections.py

@@ -0,0 +1,52 @@
+#!/usr/bin/env python
+"""
+Time: Poincaré Sections (Linear and Parallel)
+
+Computes the running time between the linear and 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 time
+import numpy as np
+import main_poincare_sections_linear as lin
+import main_poincare_sections_parallel as par
+
+E_all = np.array([1/100, 1/12, 1/10, 1/8, 1/6])
+
+par_time = []
+lin_time = []
+
+print("\033[34m" + "Please wait..." + "\033[0m")
+for E in E_all:
+    t_0 = time.time()
+    par.compute_poincare_sections_numpy(E)
+    t_1 = time.time()
+    par_time.append(t_1-t_0)
+
+print("\033[34m" + "Still wait..." + "\033[0m")
+for E in E_all:
+    t_0 = time.time()
+    lin.compute_poincare_sections_linear(E)
+    t_1 = time.time()
+    lin_time.append(t_1-t_0)
+
+print("\033[34m" + "Done!" + "\033[0m")
+print("\033[36m" + "=== [ RESULTS ] ===" + "\033[0m")
+print("\033[36m" 
+      + "- Linear algorithm:   "
+      + "\033[0m"
+      + "({:07.4f} ± {:.4f}) s".format(np.mean(lin_time), np.std(lin_time))
+      + "\033[36m"
+      + " per energy iteration"
+      + "\033[0m")
+print("\033[36m" 
+      + "- Parallel algorithm: "
+      + "\033[0m"
+      + "({:07.4f} ± {:.4f}) s".format(np.mean(par_time), np.std(lin_time))
+      + "\033[36m"
+      + " per energy iteration"
+      + "\033[0m")