From 6fb5aa63081506d4cc67c5a4b6b6225c83db1eed Mon Sep 17 00:00:00 2001
From: Elleanor Lamb <elamb@cern.ch>
Date: Mon, 15 Sep 2025 11:33:06 +0200
Subject: [PATCH 01/14] Added 4D round q-Gaussian generator for q>1, and
 example examples/particles_generation/008...

---
 .../008_generate_q_gaussian.py                |  91 +++++++++++
 xpart/__init__.py                             |   1 +
 xpart/transverse_generators/__init__.py       |   1 +
 .../transverse_generators/q_gaussian_round.py | 145 ++++++++++++++++++
 4 files changed, 238 insertions(+)
 create mode 100644 examples/particles_generation/008_generate_q_gaussian.py
 create mode 100644 xpart/transverse_generators/q_gaussian_round.py

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
new file mode 100644
index 00000000..a650e43f
--- /dev/null
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -0,0 +1,91 @@
+# copyright ############################### #
+# This file is part of the Xpart Package.   #
+# Copyright (c) CERN, 2021.                 #
+# ######################################### #
+
+import json
+import numpy as np
+from matplotlib import pyplot as plt
+import xpart as xp
+import xtrack as xt
+
+
+def q_gaussian_1d(x, q, beta, normalize=False):
+    """
+
+    Args:
+        x:
+        q: q-parameter
+        beta: beta for q-Gaussian
+        normalize: if normalize area to 1
+
+    Returns:
+        q-Gaussian function defined on x
+
+    """
+    assert q < 3, "q must be less than 3 for normalizability"
+    # Compute the argument of the power
+    arg = 1 - (1 - q) * beta * x**2
+    # Set values outside domain to 0
+    f = np.where(arg > 0, arg**(1 / (1 - q)), 0)
+    if normalize:
+        dx = x[1] - x[0]
+        area = np.sum(f) * dx
+        f /= area
+    return f
+
+
+bunch_intensity = 1e11
+sigma_z = 22.5e-2
+n_part = int(5e5)
+nemitt_x = 2e-6
+nemitt_y = 2.5e-6
+
+filename = ('../../../xtrack/test_data/sps_w_spacecharge'
+            '/line_no_spacecharge_and_particle.json')
+with open(filename, 'r') as fid:
+    ddd = json.load(fid)
+line = xt.Line.from_dict(ddd['line'])
+line.particle_ref = xp.Particles.from_dict(ddd['particle'])
+
+line.build_tracker()
+
+q = 1.1
+beta = 0.6
+
+x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_qgaussian_normalised(q=q, beta=beta, n_part=int(1e6))
+
+
+x = np.linspace(-10, 10, 1000)
+f = q_gaussian_1d(x, q, beta, normalize=True)
+
+
+# PLOT normalised x against 1D q-Gaussian
+plt.plot(x, f, color='blue', label=f'1D q-Gaussian q={q}, beta={beta}')
+plt.hist(x_norm, bins=100, density=True, label=f'sampled q-Gaussian q={q}, beta={beta}')
+plt.legend()
+plt.show()
+
+
+
+particles = line.build_particles(
+                               zeta=0, delta=1e-3,
+                               x_norm=x_norm, # in sigmas
+                               px_norm=px_norm, # in sigmas
+                               y_norm=y_norm,
+                               py_norm=py_norm,
+                               nemitt_x=3e-6, nemitt_y=3e-6)
+
+# CHECKS
+
+y_rms = np.std(particles.y)
+py_rms = np.std(particles.py)
+x_rms = np.std(particles.x)
+px_rms = np.std(particles.px)
+
+
+
+
+
+
+
diff --git a/xpart/__init__.py b/xpart/__init__.py
index 5695e575..f978328c 100644
--- a/xpart/__init__.py
+++ b/xpart/__init__.py
@@ -22,6 +22,7 @@
 from .transverse_generators import generate_2D_gaussian
 from .transverse_generators import (generate_hypersphere_2D, generate_hypersphere_4D,
                                     generate_hypersphere_6D)
+from .transverse_generators import generate_round_4D_qgaussian_normalised
 
 from .longitudinal import generate_longitudinal_coordinates
 from .longitudinal.generate_longitudinal import _characterize_line
diff --git a/xpart/transverse_generators/__init__.py b/xpart/transverse_generators/__init__.py
index d1d6e467..dc10c8bd 100644
--- a/xpart/transverse_generators/__init__.py
+++ b/xpart/transverse_generators/__init__.py
@@ -8,3 +8,4 @@
 from .pencil import generate_2D_pencil, generate_2D_pencil_with_absolute_cut
 from .gaussian import generate_2D_gaussian
 from .hypersphere import generate_hypersphere_2D, generate_hypersphere_4D, generate_hypersphere_6D
+from .q_gaussian_round import generate_round_4D_qgaussian_normalised
diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
new file mode 100644
index 00000000..113b111e
--- /dev/null
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -0,0 +1,145 @@
+#################################################
+# This code randomly samples 4D                 #
+# q-Gaussian distributions (q>1) using the      #
+# methods of Batygin                            #
+# https://doi.org/10.1016/j.nima.2004.10.029    #
+# and the 4D q-Gaussian formula derived in      #
+# https://cds.cern.ch/record/2912366?ln=en      #
+#################################################
+
+
+import numpy as np
+from scipy.special import gamma
+from scipy.interpolate import interp1d
+
+
+def generate_radial_distribution(q, beta):
+    """
+    Compute the 4D radial distribution function for a round q-Gaussian.
+
+    Parameters:
+        q (float): Entropic index (q > 1).
+        beta (float): Scale parameter.
+
+    Returns:
+        tuple: (f_F, F) where f_F is the radial distribution, and F is the radial coordinate array.
+    """
+    assert q > 1, "q must be greater than 1"
+    F = np.linspace(0, 3000, 100000)
+    term1 = -(beta**2) * (q - 3) * (q**2 - 1) / 4 / np.pi**2
+    if q < 1.01:
+        term2 = -1 / (1 - q)
+    else:
+        term2 = gamma(q / (q - 1)) / gamma(1 / (q - 1))
+
+    term3 = (1 + beta * (q - 1) * F) ** (1 / (1 - q) - 3 / 2)
+    return term1 * term2 * term3, F
+
+
+def generate_PDF(f_F, F):
+    """
+    Compute the PDF g(F) from f(F) using the Abel transform in 4D.
+
+    Parameters:
+        f_F (np.ndarray): Distribution array.
+        F (np.ndarray): Radial coordinate array.
+
+    Returns:
+        np.ndarray: Transformed PDF g(F).
+    """
+    f_F[0] = 0
+    f_F[-1] = 0
+    g_F = np.pi**2 * f_F * F
+    return g_F
+
+
+def generate_CDF(g_F, F):
+    """
+    Compute the cumulative distribution function (CDF) of g(F).
+
+    Parameters:
+        g_F (np.ndarray): PDF values.
+        F (np.ndarray): Radial coordinates.
+
+    Returns:
+        np.ndarray: CDF of g(F).
+    """
+    return np.cumsum(
+        np.diff(np.insert(F, 0, 0)) * g_F
+    )  # fast cumulative trapezoid approx
+
+
+# Functions for random sampling in 4D
+def random_beta(F_G):
+    for i in range(len(F_G)):
+        beta_x = np.random.uniform(0, 2 * np.pi, 1)
+        beta_y = np.random.uniform(0, 2 * np.pi, 1)
+    return beta_x, beta_y
+
+
+def sample_from_inv_cdf(Np, cdf_g, F):
+    """
+    Sample F values from the inverse CDF of g(F).
+
+    Parameters:
+        Np (int): Number of particles to sample.
+        cdf_g (np.ndarray): CDF of g(F).
+        F (np.ndarray): Original F grid.
+
+    Returns:
+        np.ndarray: Sampled F values (F_G).
+    """
+    cdf_g /= cdf_g[-1]  # normalize
+    uniform_samples = np.random.uniform(0, 1, Np)
+    interpolator = interp1d(
+        cdf_g, F, kind="nearest", bounds_error=False, fill_value=(F[0], F[-1])
+    )
+    return interpolator(uniform_samples)
+
+
+def generate_random_A(F_G):
+    """
+    Generate A_x and A_y coordinates based on F_G distribution.
+
+    Parameters:
+        F_G (np.ndarray): Sampled F values.
+
+    Returns:
+        tuple: (A_x, A_y) arrays.
+    """
+    A_X_SQ = np.random.uniform(0, F_G)
+    A_x = np.sqrt(A_X_SQ)
+    A_y = np.sqrt(F_G - A_X_SQ)
+    return A_x, A_y
+
+
+# function to generate a round 4D q-Gaussian
+def generate_round_4D_qgaussian_normalised(q, beta, n_part):
+    """
+    Generate particles sampled from a 4D round q-Gaussian distribution.
+
+    Parameters:
+        q (float): q-Gaussian q parameter.
+        beta (float): Scale parameter.
+        n_part (int): Number of particles to sample.
+
+    Returns:
+        tuple: Arrays of positions and momenta (x, px, y, py).
+    """
+    f_F, F = generate_radial_distribution(q, beta)  # 4D distribution
+    g_F = generate_PDF(f_F, F)  # PDF of 4D distribution
+    cdf_g = generate_CDF(g_F, F)  # CDF
+    F_G = sample_from_inv_cdf(n_part, cdf_g, F)  # Inverse function
+    A_x, A_y = generate_random_A(F_G)  # random generator distributed like F_G
+
+    # Sample angles for all particles
+    beta_x = np.random.uniform(0, 2 * np.pi, n_part)
+    beta_y = np.random.uniform(0, 2 * np.pi, n_part)
+
+    # Compute positions and momenta
+    x = A_x * np.cos(beta_x)
+    px = -A_x * np.sin(beta_x)
+    y = -A_y * np.cos(beta_y)
+    py = -A_y * np.sin(beta_y)
+
+    return x, px, y, py

From f4d5aa9e210e62b29ecea860c4f26c6f8947cdf2 Mon Sep 17 00:00:00 2001
From: Elleanor Lamb <elamb@cern.ch>
Date: Mon, 15 Sep 2025 11:34:56 +0200
Subject: [PATCH 02/14] removed obsolete beta_generator

---
 xpart/transverse_generators/q_gaussian_round.py | 13 ++-----------
 1 file changed, 2 insertions(+), 11 deletions(-)

diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 113b111e..a0b696d1 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -18,7 +18,7 @@ def generate_radial_distribution(q, beta):
     Compute the 4D radial distribution function for a round q-Gaussian.
 
     Parameters:
-        q (float): Entropic index (q > 1).
+        q (float): q-parameter (q > 1).
         beta (float): Scale parameter.
 
     Returns:
@@ -66,16 +66,7 @@ def generate_CDF(g_F, F):
     """
     return np.cumsum(
         np.diff(np.insert(F, 0, 0)) * g_F
-    )  # fast cumulative trapezoid approx
-
-
-# Functions for random sampling in 4D
-def random_beta(F_G):
-    for i in range(len(F_G)):
-        beta_x = np.random.uniform(0, 2 * np.pi, 1)
-        beta_y = np.random.uniform(0, 2 * np.pi, 1)
-    return beta_x, beta_y
-
+    )
 
 def sample_from_inv_cdf(Np, cdf_g, F):
     """

From ee5ab355215937afd4620a0515206fdcb6b2fe4e Mon Sep 17 00:00:00 2001
From: Elleanor Lamb <elamb@cern.ch>
Date: Mon, 15 Sep 2025 11:56:29 +0200
Subject: [PATCH 03/14] removed obsolete beta_generator, updated naming of
 functions

---
 xpart/transverse_generators/q_gaussian_round.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index a0b696d1..427e9aea 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -1,7 +1,7 @@
 #################################################
 # This code randomly samples 4D                 #
 # q-Gaussian distributions (q>1) using the      #
-# methods of Batygin                            #
+# samplign methods of Batygin                   #
 # https://doi.org/10.1016/j.nima.2004.10.029    #
 # and the 4D q-Gaussian formula derived in      #
 # https://cds.cern.ch/record/2912366?ln=en      #

From b4bb7df4e2d9794826ce815e2a3e8c6554325f81 Mon Sep 17 00:00:00 2001
From: Elleanor Lamb <elamb@cern.ch>
Date: Mon, 15 Sep 2025 11:56:57 +0200
Subject: [PATCH 04/14] updated example, naming of functions of q-gaussian
 generator

---
 .../008_generate_q_gaussian.py                | 16 ++++----
 xpart/__init__.py                             |  2 +-
 xpart/transverse_generators/__init__.py       |  2 +-
 .../transverse_generators/q_gaussian_round.py | 37 ++++++++++---------
 4 files changed, 28 insertions(+), 29 deletions(-)

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index a650e43f..07e80070 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -12,7 +12,6 @@
 
 def q_gaussian_1d(x, q, beta, normalize=False):
     """
-
     Args:
         x:
         q: q-parameter
@@ -50,24 +49,22 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 
 line.build_tracker()
 
-q = 1.1
-beta = 0.6
+q = 1.2
+beta = 1
 
-x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_qgaussian_normalised(q=q, beta=beta, n_part=int(1e6))
+x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_q_gaussian_normalised(q=q, beta=beta, n_part=int(1e6))
 
 
-x = np.linspace(-10, 10, 1000)
+x = np.linspace(-5, 5, 4000)
 f = q_gaussian_1d(x, q, beta, normalize=True)
 
 
 # PLOT normalised x against 1D q-Gaussian
 plt.plot(x, f, color='blue', label=f'1D q-Gaussian q={q}, beta={beta}')
-plt.hist(x_norm, bins=100, density=True, label=f'sampled q-Gaussian q={q}, beta={beta}')
+plt.hist(x_norm, bins=200, density=True, label=f'sampled q-Gaussian q={q}, beta={beta}')
 plt.legend()
 plt.show()
 
-
-
 particles = line.build_particles(
                                zeta=0, delta=1e-3,
                                x_norm=x_norm, # in sigmas
@@ -77,12 +74,13 @@ def q_gaussian_1d(x, q, beta, normalize=False):
                                nemitt_x=3e-6, nemitt_y=3e-6)
 
 # CHECKS
-
 y_rms = np.std(particles.y)
 py_rms = np.std(particles.py)
 x_rms = np.std(particles.x)
 px_rms = np.std(particles.px)
 
+print('y rms: ', y_rms, 'py rms: ', py_rms,'x rms: ', x_rms, 'px rms: ', px_rms)
+
 
 
 
diff --git a/xpart/__init__.py b/xpart/__init__.py
index f978328c..8609d0ea 100644
--- a/xpart/__init__.py
+++ b/xpart/__init__.py
@@ -22,7 +22,7 @@
 from .transverse_generators import generate_2D_gaussian
 from .transverse_generators import (generate_hypersphere_2D, generate_hypersphere_4D,
                                     generate_hypersphere_6D)
-from .transverse_generators import generate_round_4D_qgaussian_normalised
+from .transverse_generators import generate_round_4D_q_gaussian_normalised
 
 from .longitudinal import generate_longitudinal_coordinates
 from .longitudinal.generate_longitudinal import _characterize_line
diff --git a/xpart/transverse_generators/__init__.py b/xpart/transverse_generators/__init__.py
index dc10c8bd..7740ee8a 100644
--- a/xpart/transverse_generators/__init__.py
+++ b/xpart/transverse_generators/__init__.py
@@ -8,4 +8,4 @@
 from .pencil import generate_2D_pencil, generate_2D_pencil_with_absolute_cut
 from .gaussian import generate_2D_gaussian
 from .hypersphere import generate_hypersphere_2D, generate_hypersphere_4D, generate_hypersphere_6D
-from .q_gaussian_round import generate_round_4D_qgaussian_normalised
+from .q_gaussian_round import generate_round_4D_q_gaussian_normalised
diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 427e9aea..8d2c6a7b 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -1,13 +1,12 @@
 #################################################
 # This code randomly samples 4D                 #
 # q-Gaussian distributions (q>1) using the      #
-# samplign methods of Batygin                   #
+# sampling methods of Batygin                   #
 # https://doi.org/10.1016/j.nima.2004.10.029    #
 # and the 4D q-Gaussian formula derived in      #
 # https://cds.cern.ch/record/2912366?ln=en      #
 #################################################
 
-
 import numpy as np
 from scipy.special import gamma
 from scipy.interpolate import interp1d
@@ -16,6 +15,7 @@
 def generate_radial_distribution(q, beta):
     """
     Compute the 4D radial distribution function for a round q-Gaussian.
+    This can be numerically unstable if extreme values of q, beta,
 
     Parameters:
         q (float): q-parameter (q > 1).
@@ -25,7 +25,7 @@ def generate_radial_distribution(q, beta):
         tuple: (f_F, F) where f_F is the radial distribution, and F is the radial coordinate array.
     """
     assert q > 1, "q must be greater than 1"
-    F = np.linspace(0, 3000, 100000)
+    F = np.linspace(0, 3000, 100000) # can be unstable
     term1 = -(beta**2) * (q - 3) * (q**2 - 1) / 4 / np.pi**2
     if q < 1.01:
         term2 = -1 / (1 - q)
@@ -36,9 +36,10 @@ def generate_radial_distribution(q, beta):
     return term1 * term2 * term3, F
 
 
-def generate_PDF(f_F, F):
+def generate_pdf(f_F, F):
     """
     Compute the PDF g(F) from f(F) using the Abel transform in 4D.
+    Cleans up the boundaries and gives correct normalisation factor.
 
     Parameters:
         f_F (np.ndarray): Distribution array.
@@ -47,13 +48,12 @@ def generate_PDF(f_F, F):
     Returns:
         np.ndarray: Transformed PDF g(F).
     """
-    f_F[0] = 0
-    f_F[-1] = 0
-    g_F = np.pi**2 * f_F * F
-    return g_F
+    f_F = f_F.copy()
+    f_F[0] = f_F[-1] = 0  # Boundary cleanup
+    return np.pi**2 * f_F * F
 
 
-def generate_CDF(g_F, F):
+def generate_cdf(g_F, F):
     """
     Compute the cumulative distribution function (CDF) of g(F).
 
@@ -64,9 +64,11 @@ def generate_CDF(g_F, F):
     Returns:
         np.ndarray: CDF of g(F).
     """
-    return np.cumsum(
-        np.diff(np.insert(F, 0, 0)) * g_F
-    )
+    delta_F = np.diff(F, prepend=0)
+    cdf = np.cumsum(g_F * delta_F)
+    cdf = np.clip(cdf, 0, np.inf)
+    return cdf
+
 
 def sample_from_inv_cdf(Np, cdf_g, F):
     """
@@ -88,7 +90,7 @@ def sample_from_inv_cdf(Np, cdf_g, F):
     return interpolator(uniform_samples)
 
 
-def generate_random_A(F_G):
+def generate_random_a(F_G):
     """
     Generate A_x and A_y coordinates based on F_G distribution.
 
@@ -104,8 +106,7 @@ def generate_random_A(F_G):
     return A_x, A_y
 
 
-# function to generate a round 4D q-Gaussian
-def generate_round_4D_qgaussian_normalised(q, beta, n_part):
+def generate_round_4D_q_gaussian_normalised(q, beta, n_part):
     """
     Generate particles sampled from a 4D round q-Gaussian distribution.
 
@@ -118,10 +119,10 @@ def generate_round_4D_qgaussian_normalised(q, beta, n_part):
         tuple: Arrays of positions and momenta (x, px, y, py).
     """
     f_F, F = generate_radial_distribution(q, beta)  # 4D distribution
-    g_F = generate_PDF(f_F, F)  # PDF of 4D distribution
-    cdf_g = generate_CDF(g_F, F)  # CDF
+    g_F = generate_pdf(f_F, F)  # PDF of 4D distribution
+    cdf_g = generate_cdf(g_F, F)  # CDF
     F_G = sample_from_inv_cdf(n_part, cdf_g, F)  # Inverse function
-    A_x, A_y = generate_random_A(F_G)  # random generator distributed like F_G
+    A_x, A_y = generate_random_a(F_G)  # random generator distributed like F_G
 
     # Sample angles for all particles
     beta_x = np.random.uniform(0, 2 * np.pi, n_part)

From 9a80e14b46a5b55399b8edfc111bc11bcb104907 Mon Sep 17 00:00:00 2001
From: Elleanor Lamb <elamb@cern.ch>
Date: Mon, 15 Sep 2025 12:09:09 +0200
Subject: [PATCH 05/14] removed comma in docs

---
 xpart/transverse_generators/q_gaussian_round.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 8d2c6a7b..414355af 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -15,7 +15,7 @@
 def generate_radial_distribution(q, beta):
     """
     Compute the 4D radial distribution function for a round q-Gaussian.
-    This can be numerically unstable if extreme values of q, beta,
+    This can be numerically unstable if extreme values of q, beta.
 
     Parameters:
         q (float): q-parameter (q > 1).
@@ -135,3 +135,5 @@ def generate_round_4D_q_gaussian_normalised(q, beta, n_part):
     py = -A_y * np.sin(beta_y)
 
     return x, px, y, py
+
+

From a174dadf6d6d6b3e0ea66063b85df18bd0e5fa5c Mon Sep 17 00:00:00 2001
From: Elleanor Lamb <elamb@cern.ch>
Date: Mon, 15 Sep 2025 14:59:12 +0200
Subject: [PATCH 06/14] added test cases

---
 examples/particles_generation/008_generate_q_gaussian.py | 5 ++---
 tests/test_transverse_q_gaussian_4d.py                   | 0
 xpart/transverse_generators/q_gaussian_round.py          | 6 +++---
 3 files changed, 5 insertions(+), 6 deletions(-)
 create mode 100644 tests/test_transverse_q_gaussian_4d.py

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index 07e80070..5da820d6 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -55,11 +55,10 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_q_gaussian_normalised(q=q, beta=beta, n_part=int(1e6))
 
 
-x = np.linspace(-5, 5, 4000)
-f = q_gaussian_1d(x, q, beta, normalize=True)
-
 
 # PLOT normalised x against 1D q-Gaussian
+x = np.linspace(-10, 10, 1000)
+f = q_gaussian_1d(x=x, q=q, beta=beta, normalize=True)
 plt.plot(x, f, color='blue', label=f'1D q-Gaussian q={q}, beta={beta}')
 plt.hist(x_norm, bins=200, density=True, label=f'sampled q-Gaussian q={q}, beta={beta}')
 plt.legend()
diff --git a/tests/test_transverse_q_gaussian_4d.py b/tests/test_transverse_q_gaussian_4d.py
new file mode 100644
index 00000000..e69de29b
diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 414355af..17412c89 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -27,7 +27,7 @@ def generate_radial_distribution(q, beta):
     assert q > 1, "q must be greater than 1"
     F = np.linspace(0, 3000, 100000) # can be unstable
     term1 = -(beta**2) * (q - 3) * (q**2 - 1) / 4 / np.pi**2
-    if q < 1.01:
+    if abs(q - 1) < 1e-2:
         term2 = -1 / (1 - q)
     else:
         term2 = gamma(q / (q - 1)) / gamma(1 / (q - 1))
@@ -65,7 +65,7 @@ def generate_cdf(g_F, F):
         np.ndarray: CDF of g(F).
     """
     delta_F = np.diff(F, prepend=0)
-    cdf = np.cumsum(g_F * delta_F)
+    cdf = np.cumsum(g_F * delta_F) # todo: rewrite with scipy.integrate.quad(g_F, -np.inf, np.inf)
     cdf = np.clip(cdf, 0, np.inf)
     return cdf
 
@@ -132,7 +132,7 @@ def generate_round_4D_q_gaussian_normalised(q, beta, n_part):
     x = A_x * np.cos(beta_x)
     px = -A_x * np.sin(beta_x)
     y = -A_y * np.cos(beta_y)
-    py = -A_y * np.sin(beta_y)
+    py = A_y * np.sin(beta_y)
 
     return x, px, y, py
 

From eb3264d7d0c5a6b047255f11735115d9b09233c5 Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Wed, 1 Oct 2025 17:53:37 +0200
Subject: [PATCH 07/14] added unit tests, improved doc strings

---
 .../008_generate_q_gaussian.py                |  4 +-
 tests/test_transverse_q_gaussian_4d.py        | 66 +++++++++++++++++++
 .../transverse_generators/q_gaussian_round.py | 12 ++--
 3 files changed, 74 insertions(+), 8 deletions(-)

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index 5da820d6..d35e4089 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -22,10 +22,8 @@ def q_gaussian_1d(x, q, beta, normalize=False):
         q-Gaussian function defined on x
 
     """
-    assert q < 3, "q must be less than 3 for normalizability"
-    # Compute the argument of the power
+    assert q < 5/3, "q must be less than 5/3"
     arg = 1 - (1 - q) * beta * x**2
-    # Set values outside domain to 0
     f = np.where(arg > 0, arg**(1 / (1 - q)), 0)
     if normalize:
         dx = x[1] - x[0]
diff --git a/tests/test_transverse_q_gaussian_4d.py b/tests/test_transverse_q_gaussian_4d.py
index e69de29b..4f05a432 100644
--- a/tests/test_transverse_q_gaussian_4d.py
+++ b/tests/test_transverse_q_gaussian_4d.py
@@ -0,0 +1,66 @@
+import numpy as np
+from scipy.optimize import curve_fit
+
+from xpart.transverse_generators.q_gaussian_round import (
+    generate_round_4D_q_gaussian_normalised)
+
+def test_transverse_q_gaussian_4d_sample_std():
+    """
+    test standard deviation is correct
+    """
+    q = 1.4
+    beta = 1
+    n_part = 10000000
+    x, _, _, _ = generate_round_4D_q_gaussian_normalised(q, beta, n_part)
+
+    # variance a function of q and beta
+    sample_variance = np.std(x)
+    expected_variance = np.sqrt(1 / (beta * (5 - 3 * q)))
+
+    assert np.isclose(sample_variance, expected_variance, atol=0.05), \
+        f"Sample variance {sample_variance} deviates from expected {expected_variance}"
+
+def q_gaussian_1d(x, q, beta, normalize=False):
+    """
+    Args:
+        x:
+        q: q-parameter
+        beta: beta for q-Gaussian
+        normalize: if normalize area to 1
+
+    Returns:
+        q-Gaussian function defined on x
+
+    """
+    assert q < 5/3, "q must be less than 5/3"
+    arg = 1 - (1 - q) * beta * x**2
+    f = np.where(arg > 0, arg**(1 / (1 - q)), 0)
+    if normalize:
+        dx = x[1] - x[0]
+        area = np.sum(f) * dx
+        f /= area
+    return f
+
+def test_transverse_q_gaussian_4d_sampler_returns_q0():
+    """
+    test that required q matches fitted q from scipy.optimize.curve fit
+    """
+    q = 1.4
+    beta = 1
+    n_part = 10000000
+    x, _, _, _ = generate_round_4D_q_gaussian_normalised(q, beta, n_part)
+    # Generate histogram
+    bins = np.linspace(-10, 10, 1000)  # 1000 bins between -10 and 10
+    counts, bin_edges = np.histogram(x, bins=bins, density=True)  # density=True to normalize to PDF
+
+    # Calculate bin centers for fitting
+    bin_centers = 0.5 * (bin_edges[:-1] + bin_edges[1:])
+    popt, pcov = curve_fit(q_gaussian_1d, bin_centers, counts, p0=[1.5, 1, 1.0])
+
+    assert np.isclose(popt[0], q, atol=0.05), \
+        f"Fitted q from samples {popt[0]} does not match requested q {q} "
+
+
+# test that variable F(q, beta) works?
+
+
diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 17412c89..62434f33 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -1,10 +1,11 @@
 #################################################
 # This code randomly samples 4D                 #
-# q-Gaussian distributions (q>1) using the      #
+# q-Gaussian distributions (1<q<5/3) using the  #
 # sampling methods of Batygin                   #
 # https://doi.org/10.1016/j.nima.2004.10.029    #
 # and the 4D q-Gaussian formula derived in      #
 # https://cds.cern.ch/record/2912366?ln=en      #
+# These distributions are non-factorizable.     #
 #################################################
 
 import numpy as np
@@ -12,20 +13,21 @@
 from scipy.interpolate import interp1d
 
 
-def generate_radial_distribution(q, beta):
+def generate_radial_distribution(q:float, beta:float):
     """
     Compute the 4D radial distribution function for a round q-Gaussian.
-    This can be numerically unstable if extreme values of q, beta.
+    This can be numerically unstable if extreme values of q, beta as the
+    sample space needs to be increased (eg for very high q or beta)
 
     Parameters:
-        q (float): q-parameter (q > 1).
+        q (float): q-parameter (1 < q < 5/3).
         beta (float): Scale parameter.
 
     Returns:
         tuple: (f_F, F) where f_F is the radial distribution, and F is the radial coordinate array.
     """
     assert q > 1, "q must be greater than 1"
-    F = np.linspace(0, 3000, 100000) # can be unstable
+    F = np.linspace(0, 30000, 1000000) # can be unstable
     term1 = -(beta**2) * (q - 3) * (q**2 - 1) / 4 / np.pi**2
     if abs(q - 1) < 1e-2:
         term2 = -1 / (1 - q)

From 12b5d40379a76631c142fe09422d592b7bbde3b0 Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Wed, 1 Oct 2025 18:07:21 +0200
Subject: [PATCH 08/14] added user defined optional sample_space parameter for
 large q and beta

---
 .../transverse_generators/q_gaussian_round.py | 19 ++++++++++++++-----
 1 file changed, 14 insertions(+), 5 deletions(-)

diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 62434f33..76bf4ef8 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -13,7 +13,7 @@
 from scipy.interpolate import interp1d
 
 
-def generate_radial_distribution(q:float, beta:float):
+def generate_radial_distribution(q, beta, F):
     """
     Compute the 4D radial distribution function for a round q-Gaussian.
     This can be numerically unstable if extreme values of q, beta as the
@@ -21,13 +21,16 @@ def generate_radial_distribution(q:float, beta:float):
 
     Parameters:
         q (float): q-parameter (1 < q < 5/3).
-        beta (float): Scale parameter.
+        beta (float): Scale parameter, (beta > 1)
+        F: sample space parameter, can be user defined
 
     Returns:
         tuple: (f_F, F) where f_F is the radial distribution, and F is the radial coordinate array.
     """
     assert q > 1, "q must be greater than 1"
-    F = np.linspace(0, 30000, 1000000) # can be unstable
+    assert q < 5/3, "q must be less than 5/3"
+    assert beta > 0, "beta must be greater than 0"
+
     term1 = -(beta**2) * (q - 3) * (q**2 - 1) / 4 / np.pi**2
     if abs(q - 1) < 1e-2:
         term2 = -1 / (1 - q)
@@ -108,7 +111,7 @@ def generate_random_a(F_G):
     return A_x, A_y
 
 
-def generate_round_4D_q_gaussian_normalised(q, beta, n_part):
+def generate_round_4D_q_gaussian_normalised(q, beta, n_part, sample_space=None):
     """
     Generate particles sampled from a 4D round q-Gaussian distribution.
 
@@ -116,11 +119,17 @@ def generate_round_4D_q_gaussian_normalised(q, beta, n_part):
         q (float): q-Gaussian q parameter.
         beta (float): Scale parameter.
         n_part (int): Number of particles to sample.
+        sample_space: Default np.linspace(0, 30000, 1000000)
 
     Returns:
         tuple: Arrays of positions and momenta (x, px, y, py).
     """
-    f_F, F = generate_radial_distribution(q, beta)  # 4D distribution
+    if sample_space is None:
+        F = np.linspace(0, 30000, 1000000)
+    else:
+        F = sample_space
+
+    f_F, F = generate_radial_distribution(q, beta, F)  # 4D distribution
     g_F = generate_pdf(f_F, F)  # PDF of 4D distribution
     cdf_g = generate_cdf(g_F, F)  # CDF
     F_G = sample_from_inv_cdf(n_part, cdf_g, F)  # Inverse function

From 3b1b43870512cf419f0b32efbfc88713a3cd64f9 Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Wed, 1 Oct 2025 18:56:22 +0200
Subject: [PATCH 09/14] change plots

---
 .../008_generate_q_gaussian.py                | 33 ++++++++++++++-----
 1 file changed, 25 insertions(+), 8 deletions(-)

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index d35e4089..1bba3216 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -34,7 +34,7 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 
 bunch_intensity = 1e11
 sigma_z = 22.5e-2
-n_part = int(5e5)
+n_part = int(5e6)
 nemitt_x = 2e-6
 nemitt_y = 2.5e-6
 
@@ -47,20 +47,15 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 
 line.build_tracker()
 
-q = 1.2
+q = 1.3
 beta = 1
 
 x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_q_gaussian_normalised(q=q, beta=beta, n_part=int(1e6))
 
-
-
 # PLOT normalised x against 1D q-Gaussian
 x = np.linspace(-10, 10, 1000)
 f = q_gaussian_1d(x=x, q=q, beta=beta, normalize=True)
-plt.plot(x, f, color='blue', label=f'1D q-Gaussian q={q}, beta={beta}')
-plt.hist(x_norm, bins=200, density=True, label=f'sampled q-Gaussian q={q}, beta={beta}')
-plt.legend()
-plt.show()
+
 
 particles = line.build_particles(
                                zeta=0, delta=1e-3,
@@ -78,6 +73,28 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 
 print('y rms: ', y_rms, 'py rms: ', py_rms,'x rms: ', x_rms, 'px rms: ', px_rms)
 
+plt.close('all')
+fig1 = plt.figure(1, figsize=(6.4, 7))
+ax21 = fig1.add_subplot(3,1,1)
+ax22 = fig1.add_subplot(3,1,2)
+ax23 = fig1.add_subplot(3,1,3)
+ax21.plot(particles.x*1000, particles.px, '.', markersize=1)
+ax21.set_xlabel(r'x [mm]')
+ax21.set_ylabel(r'px [-]')
+ax22.plot(particles.y*1000, particles.py, '.', markersize=1)
+ax22.set_xlabel(r'y [mm]')
+ax22.set_ylabel(r'py [-]')
+ax23.plot(x, f, color='k', label=f'1D q-Gaussian q={q}, beta={beta}')
+ax23.hist(x_norm, bins=200, density=True, label=f'normalised sampled q-Gaussian q={q}, beta={beta}')
+ax23.set_xlabel(r'normalised x')
+ax23.set_xlabel(r'normalised amplitude')
+ax23.legend()
+
+fig1.subplots_adjust(bottom=.08, top=.93, hspace=.33, left=.18,
+                     right=.96, wspace=.33)
+plt.show()
+
+
 
 
 

From 9ded2ad16b8042d6825a8ddddab6765aef0d8347 Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Thu, 2 Oct 2025 09:17:24 +0200
Subject: [PATCH 10/14] change plots in q-gaussian example to have x and y
 projection

---
 examples/particles_generation/008_generate_q_gaussian.py | 4 +++-
 1 file changed, 3 insertions(+), 1 deletion(-)

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index 1bba3216..9b5e330f 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -85,7 +85,9 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 ax22.set_xlabel(r'y [mm]')
 ax22.set_ylabel(r'py [-]')
 ax23.plot(x, f, color='k', label=f'1D q-Gaussian q={q}, beta={beta}')
-ax23.hist(x_norm, bins=200, density=True, label=f'normalised sampled q-Gaussian q={q}, beta={beta}')
+ax23.hist(x_norm, bins=200, density=True,alpha=0.5, label=f'normalised sampled q-Gaussian q={q}, beta={beta}, $x$ projection')
+ax23.hist(y_norm, bins=200, density=True,alpha=0.5, label=f'normalised sampled q-Gaussian q={q}, beta={beta} $y$ projection')
+
 ax23.set_xlabel(r'normalised x')
 ax23.set_xlabel(r'normalised amplitude')
 ax23.legend()

From 563c9cf075807a26acbe0f8e41158bd01e320864 Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Thu, 2 Oct 2025 09:40:55 +0200
Subject: [PATCH 11/14] reformat example

---
 .../008_generate_q_gaussian.py                | 22 +++++++++++--------
 1 file changed, 13 insertions(+), 9 deletions(-)

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index 9b5e330f..5c2c40ee 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -50,12 +50,11 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 q = 1.3
 beta = 1
 
-x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_q_gaussian_normalised(q=q, beta=beta, n_part=int(1e6))
-
-# PLOT normalised x against 1D q-Gaussian
-x = np.linspace(-10, 10, 1000)
-f = q_gaussian_1d(x=x, q=q, beta=beta, normalize=True)
-
+x_norm, px_norm, y_norm, py_norm = xp.generate_round_4D_q_gaussian_normalised(
+    q=q,
+    beta=beta,
+    n_part=int(1e6)
+)
 
 particles = line.build_particles(
                                zeta=0, delta=1e-3,
@@ -73,6 +72,10 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 
 print('y rms: ', y_rms, 'py rms: ', py_rms,'x rms: ', x_rms, 'px rms: ', px_rms)
 
+# plot 1D q-gaussian against projection
+x = np.linspace(-10, 10, 1000)
+f = q_gaussian_1d(x=x, q=q, beta=beta, normalize=True)
+
 plt.close('all')
 fig1 = plt.figure(1, figsize=(6.4, 7))
 ax21 = fig1.add_subplot(3,1,1)
@@ -85,9 +88,10 @@ def q_gaussian_1d(x, q, beta, normalize=False):
 ax22.set_xlabel(r'y [mm]')
 ax22.set_ylabel(r'py [-]')
 ax23.plot(x, f, color='k', label=f'1D q-Gaussian q={q}, beta={beta}')
-ax23.hist(x_norm, bins=200, density=True,alpha=0.5, label=f'normalised sampled q-Gaussian q={q}, beta={beta}, $x$ projection')
-ax23.hist(y_norm, bins=200, density=True,alpha=0.5, label=f'normalised sampled q-Gaussian q={q}, beta={beta} $y$ projection')
-
+ax23.hist(x_norm, bins=200, density=True, alpha=0.5,
+          label=f'normalised sampled q-Gaussian q={q}, beta={beta}, $x$ projection')
+ax23.hist(y_norm, bins=200, density=True, alpha=0.5,
+          label=f'normalised sampled q-Gaussian q={q}, beta={beta} $y$ projection')
 ax23.set_xlabel(r'normalised x')
 ax23.set_xlabel(r'normalised amplitude')
 ax23.legend()

From f2bbe9e9f891fb9978606c606938e8ba51c07f5b Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Thu, 2 Oct 2025 09:42:13 +0200
Subject: [PATCH 12/14] reformat example

---
 .../particles_generation/008_generate_q_gaussian.py    | 10 ----------
 1 file changed, 10 deletions(-)

diff --git a/examples/particles_generation/008_generate_q_gaussian.py b/examples/particles_generation/008_generate_q_gaussian.py
index 5c2c40ee..0aa9aff8 100644
--- a/examples/particles_generation/008_generate_q_gaussian.py
+++ b/examples/particles_generation/008_generate_q_gaussian.py
@@ -9,7 +9,6 @@
 import xpart as xp
 import xtrack as xt
 
-
 def q_gaussian_1d(x, q, beta, normalize=False):
     """
     Args:
@@ -31,7 +30,6 @@ def q_gaussian_1d(x, q, beta, normalize=False):
         f /= area
     return f
 
-
 bunch_intensity = 1e11
 sigma_z = 22.5e-2
 n_part = int(5e6)
@@ -64,14 +62,6 @@ def q_gaussian_1d(x, q, beta, normalize=False):
                                py_norm=py_norm,
                                nemitt_x=3e-6, nemitt_y=3e-6)
 
-# CHECKS
-y_rms = np.std(particles.y)
-py_rms = np.std(particles.py)
-x_rms = np.std(particles.x)
-px_rms = np.std(particles.px)
-
-print('y rms: ', y_rms, 'py rms: ', py_rms,'x rms: ', x_rms, 'px rms: ', px_rms)
-
 # plot 1D q-gaussian against projection
 x = np.linspace(-10, 10, 1000)
 f = q_gaussian_1d(x=x, q=q, beta=beta, normalize=True)

From 0fbaff0d8c75565ebe24b0ee2c2ffd4052edef2f Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Thu, 2 Oct 2025 09:58:54 +0200
Subject: [PATCH 13/14] updated docstring for
 generate_round_4D_q_gaussian_normalised

---
 .../transverse_generators/q_gaussian_round.py | 30 ++++++++++++++-----
 1 file changed, 22 insertions(+), 8 deletions(-)

diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index 76bf4ef8..ed117bfe 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -116,16 +116,30 @@ def generate_round_4D_q_gaussian_normalised(q, beta, n_part, sample_space=None):
     Generate particles sampled from a 4D round q-Gaussian distribution.
 
     Parameters:
-        q (float): q-Gaussian q parameter.
-        beta (float): Scale parameter.
-        n_part (int): Number of particles to sample.
-        sample_space: Default np.linspace(0, 30000, 1000000)
-
-    Returns:
-        tuple: Arrays of positions and momenta (x, px, y, py).
+    -----------
+    q : float
+        The q parameter of the q-Gaussian distribution. Must satisfy 1 < q < 5/3.
+    beta : float
+        Scale parameter (analogous to beta function) for the distribution.
+    n_part : int
+        Number of particles to generate.
+    sample_space : array_like, optional
+        1D array of radius values at which to evaluate the radial PDF/CDF.
+        If None, defaults to `np.linspace(0, 300, 1000000)`.
+
+    Returns
+    -------
+    x : np.ndarray
+        Horizontal position coordinates.
+    px : np.ndarray
+        Horizontal momentum coordinates.
+    y : np.ndarray
+        Vertical position coordinates.
+    py : np.ndarray
+        Vertical momentum coordinates.
     """
     if sample_space is None:
-        F = np.linspace(0, 30000, 1000000)
+        F = np.linspace(0, 3e2, 1e6)
     else:
         F = sample_space
 

From 6ba814bbaa33473773c643ea69b08689bab3acc4 Mon Sep 17 00:00:00 2001
From: elamb <elamb@cern.ch>
Date: Thu, 2 Oct 2025 10:05:28 +0200
Subject: [PATCH 14/14] updated docstring for
 generate_round_4D_q_gaussian_normalised

---
 .../transverse_generators/q_gaussian_round.py | 29 +++++--------------
 1 file changed, 8 insertions(+), 21 deletions(-)

diff --git a/xpart/transverse_generators/q_gaussian_round.py b/xpart/transverse_generators/q_gaussian_round.py
index ed117bfe..ce34e578 100644
--- a/xpart/transverse_generators/q_gaussian_round.py
+++ b/xpart/transverse_generators/q_gaussian_round.py
@@ -116,27 +116,14 @@ def generate_round_4D_q_gaussian_normalised(q, beta, n_part, sample_space=None):
     Generate particles sampled from a 4D round q-Gaussian distribution.
 
     Parameters:
-    -----------
-    q : float
-        The q parameter of the q-Gaussian distribution. Must satisfy 1 < q < 5/3.
-    beta : float
-        Scale parameter (analogous to beta function) for the distribution.
-    n_part : int
-        Number of particles to generate.
-    sample_space : array_like, optional
-        1D array of radius values at which to evaluate the radial PDF/CDF.
-        If None, defaults to `np.linspace(0, 300, 1000000)`.
-
-    Returns
-    -------
-    x : np.ndarray
-        Horizontal position coordinates.
-    px : np.ndarray
-        Horizontal momentum coordinates.
-    y : np.ndarray
-        Vertical position coordinates.
-    py : np.ndarray
-        Vertical momentum coordinates.
+        q (float): q-Gaussian shape parameter. Must satisfy 1 < q < 5/3.
+        beta (float): Scale parameter.
+        n_part (int): Number of particles to generate.
+        sample_space (array-like, optional): 1D array of radius values used for sampling.
+            Defaults to np.linspace(0, 3000, 100000) if None.
+
+    Returns:
+        tuple: Arrays of normalised transverse coordinates (x, px, y, py).
     """
     if sample_space is None:
         F = np.linspace(0, 3e2, 1e6)