From 2a4a569511e782b65ef244074ba6278d3ad7fb83 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 08:45:35 +0100
Subject: [PATCH 01/18] refactor of the independent samplings to use the
 distributions and IMFs

---
 posydon/popsyn/IMFs.py               |  69 ++++
 posydon/popsyn/distributions.py      | 578 +++++++++++++++++++++++++++
 posydon/popsyn/independent_sample.py | 119 ++----
 posydon/popsyn/norm_pop.py           |  12 +-
 posydon/utils/common_functions.py    |  34 +-
 5 files changed, 715 insertions(+), 97 deletions(-)

diff --git a/posydon/popsyn/IMFs.py b/posydon/popsyn/IMFs.py
index d144470a75..62055864fb 100644
--- a/posydon/popsyn/IMFs.py
+++ b/posydon/popsyn/IMFs.py
@@ -166,6 +166,39 @@ def imf(self, m):
         valid = self._check_valid(m)
         return m ** (-self.alpha)
 
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the Salpeter IMF.
+
+        Uses analytical inverse transform sampling for efficiency.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random mass samples in solar masses.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Analytical inverse transform sampling
+        # For power law IMF: m^(-alpha), the inverse CDF is:
+        # m = (u * (m_max^(1-alpha) - m_min^(1-alpha)) + m_min^(1-alpha))^(1/(1-alpha))
+        normalization_constant = (1.0 - self.alpha) / (
+            self.m_max**(1.0 - self.alpha) - self.m_min**(1.0 - self.alpha)
+        )
+        u = rng.uniform(size=size)
+        masses = (u * (1.0 - self.alpha) / normalization_constant
+                  + self.m_min**(1.0 - self.alpha))**(1.0 / (1.0 - self.alpha))
+
+        return masses
+
+
 class Kroupa2001(IMFBase):
     """Initial Mass Function based on Kroupa (2001), which is
     defined as a broken power-law:
@@ -275,6 +308,42 @@ def imf(self, m):
         out[mask3] = const2 * (m[mask3] / self.m2break) ** (-self.alpha3)
         return out
 
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the Kroupa2001 IMF.
+
+        Uses inverse transform sampling with discretized PDF for the
+        broken power-law distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random mass samples in solar masses.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Import here to avoid circular dependency
+        from posydon.utils.common_functions import inverse_sampler
+
+        # Create discretized PDF for inverse sampling
+        # Use more points near the breaks for better accuracy
+        n_points = 2000
+        m_grid = np.linspace(self.m_min, self.m_max, n_points)
+        pdf_values = self.imf(m_grid)
+
+        # Sample using inverse transform method
+        masses = inverse_sampler(m_grid, pdf_values, size=size, rng=rng)
+
+        return masses
+
+
 class Chabrier2003(IMFBase):
     """Chabrier2003 Initial Mass Function (IMF), which is defined as
     a lognormal distribution for low-mass stars and a power-law distribution
diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index 36247403ab..82e721f804 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -5,6 +5,7 @@
 ]
 import numpy as np
 from scipy.integrate import quad
+from scipy.stats import truncnorm
 
 
 class FlatMassRatio:
@@ -129,6 +130,26 @@ def pdf(self, q):
         pdf_values[valid] = self.flat_mass_ratio(q[valid]) * self.norm
         return pdf_values
 
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the flat mass ratio distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        float or ndarray
+            Random samples from the distribution.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        return rng.uniform(self.q_min, self.q_max, size=size)
+
 
 class Sana12Period():
     """Period distribution from Sana et al. (2012).
@@ -310,6 +331,72 @@ def pdf(self, p, m1):
                              * self.norm(m1[valid]))
         return pdf_values
 
+    def rvs(self, size=1, m1=None, rng=None):
+        """Draw random samples from the Sana12 period distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        m1 : float or array_like, optional
+            Primary mass(es). If array, must have length equal to size.
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random period samples in days.
+
+        Raises
+        ------
+        ValueError
+            If m1 is None or has incorrect size.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        if m1 is None:
+            raise ValueError("m1 (primary mass) must be provided for Sana12Period sampling")
+
+        m1 = np.atleast_1d(m1)
+        if m1.size == 1:
+            m1 = np.full(size, m1[0])
+        elif m1.size != size:
+            raise ValueError(f"m1 must be a single value or have size={size}")
+
+        # Import here to avoid circular dependency
+        from posydon.utils.common_functions import rejection_sampler
+
+        periods = np.zeros(size)
+
+        # For low mass stars (m1 <= 15), use log-uniform distribution
+        low_mass_mask = m1 <= self.mbreak
+        n_low = np.sum(low_mass_mask)
+        if n_low > 0:
+            periods[low_mass_mask] = 10**rng.uniform(
+                np.log10(self.p_min),
+                np.log10(self.p_max),
+                size=n_low
+            )
+
+        # For high mass stars (m1 > 15), use rejection sampling
+        high_mass_mask = ~low_mass_mask
+        n_high = np.sum(high_mass_mask)
+        if n_high > 0:
+            # Create PDF function for rejection sampler
+            def pdf_high_mass(logp):
+                return self.sana12_period(logp, self.mbreak + 1)
+
+            periods[high_mass_mask] = 10**rejection_sampler(
+                size=n_high,
+                x_lim=[np.log10(self.p_min), np.log10(self.p_max)],
+                pdf=pdf_high_mass,
+                rng=rng
+            )
+
+        return periods
+
 
 class PowerLawPeriod():
     '''Power law period distribution with slope pi and boundaries m_min, m_max.
@@ -438,6 +525,38 @@ def pdf(self, p):
 
         return pdf_values
 
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the power law period distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random period samples in days.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Import here to avoid circular dependency
+        from posydon.utils.common_functions import inverse_sampler
+
+        # Create discretized PDF for inverse sampling
+        n_points = 1000
+        logp_grid = np.linspace(np.log10(self.p_min), np.log10(self.p_max), n_points)
+        pdf_values = self.power_law_period(logp_grid)
+
+        # Sample in log space
+        logp_samples = inverse_sampler(logp_grid, pdf_values, size=size, rng=rng)
+
+        # Convert back to linear space
+        return 10**logp_samples
+
 
 class LogUniform():
     """Log-uniform distribution between specified minimum and maximum values.
@@ -500,3 +619,462 @@ def pdf(self, x):
         pdf_values[valid] = self.norm
 
         return pdf_values
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the log-uniform distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random samples from the distribution.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Sample uniformly in log space
+        log_samples = rng.uniform(np.log10(self.min), np.log10(self.max), size=size)
+
+        # Convert back to linear space
+        return 10**log_samples
+
+
+class ThermalEccentricity:
+    """Thermal eccentricity distribution for binary star systems.
+
+    The thermal distribution follows pdf(e) = 2*e, which is the
+    distribution expected for binaries that have undergone significant
+    dynamical interactions or thermal relaxation.
+
+    Parameters
+    ----------
+    e_min : float, optional
+        Minimum eccentricity (default: 0.0). Must be in [0, 1).
+    e_max : float, optional
+        Maximum eccentricity (default: 1.0). Must be in (0, 1].
+
+    Raises
+    ------
+    ValueError
+        If e_min or e_max are not in [0, 1], or if e_min >= e_max.
+
+    Examples
+    --------
+    >>> dist = ThermalEccentricity()
+    >>> e_samples = dist.rvs(size=1000)
+    """
+
+    def __init__(self, e_min=0.0, e_max=1.0):
+        """Initialize the thermal eccentricity distribution.
+
+        Parameters
+        ----------
+        e_min : float, optional
+            Minimum eccentricity (default: 0.0). Must be in [0, 1).
+        e_max : float, optional
+            Maximum eccentricity (default: 1.0). Must be in (0, 1].
+
+        Raises
+        ------
+        ValueError
+            If e_min or e_max are not in [0, 1], or if e_min >= e_max.
+        """
+        if not (0 <= e_min < 1):
+            raise ValueError("e_min must be in [0, 1)")
+        if not (0 < e_max <= 1):
+            raise ValueError("e_max must be in (0, 1]")
+        if e_min >= e_max:
+            raise ValueError("e_min must be less than e_max")
+
+        self.e_min = e_min
+        self.e_max = e_max
+        self.norm = self._calculate_normalization()
+
+    def __repr__(self):
+        """Return string representation of the distribution."""
+        return f"ThermalEccentricity(e_min={self.e_min}, e_max={self.e_max})"
+
+    def _repr_html_(self):
+        """Return HTML representation for Jupyter notebooks."""
+        return (f"<h3>Thermal Eccentricity Distribution</h3>"
+                f"<p>e_min = {self.e_min}</p>"
+                f"<p>e_max = {self.e_max}</p>")
+
+    def _calculate_normalization(self):
+        """Calculate the normalization constant for the thermal eccentricity
+        distribution.
+
+        Returns
+        -------
+        float
+            The normalization constant ensuring the PDF integrates to 1.
+        """
+        # Integral of 2*e from e_min to e_max is e_max^2 - e_min^2
+        integral = self.e_max**2 - self.e_min**2
+        if integral == 0:
+            raise ValueError("Cannot normalize distribution: e_min == e_max")
+        return 1.0 / integral
+
+    def thermal_eccentricity(self, e):
+        """Compute the thermal eccentricity distribution value.
+
+        Parameters
+        ----------
+        e : float or array_like
+            Eccentricity value(s).
+
+        Returns
+        -------
+        float or ndarray
+            Distribution value (2*e).
+        """
+        return 2.0 * np.asarray(e)
+
+    def pdf(self, e):
+        """Probability density function of the thermal eccentricity distribution.
+
+        Parameters
+        ----------
+        e : float or array_like
+            Eccentricity value(s).
+
+        Returns
+        -------
+        float or ndarray
+            Probability density at eccentricity e.
+        """
+        e = np.asarray(e)
+        valid = (e >= self.e_min) & (e <= self.e_max)
+        pdf_values = np.zeros_like(e, dtype=float)
+        pdf_values[valid] = self.thermal_eccentricity(e[valid]) * self.norm
+        return pdf_values
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the thermal eccentricity distribution.
+
+        Uses the analytical inverse CDF: e = sqrt(u * (e_max^2 - e_min^2) + e_min^2)
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random eccentricity samples in [e_min, e_max].
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Inverse CDF: e = sqrt(u * (e_max^2 - e_min^2) + e_min^2)
+        u = rng.uniform(size=size)
+        return np.sqrt(u * (self.e_max**2 - self.e_min**2) + self.e_min**2)
+
+
+class UniformEccentricity:
+    """Uniform eccentricity distribution for binary star systems.
+
+    A flat distribution over eccentricities between e_min and e_max.
+
+    Parameters
+    ----------
+    e_min : float, optional
+        Minimum eccentricity (default: 0.0). Must be in [0, 1).
+    e_max : float, optional
+        Maximum eccentricity (default: 1.0). Must be in (0, 1].
+
+    Raises
+    ------
+    ValueError
+        If e_min or e_max are not in [0, 1], or if e_min >= e_max.
+
+    Examples
+    --------
+    >>> dist = UniformEccentricity()
+    >>> e_samples = dist.rvs(size=1000)
+    """
+
+    def __init__(self, e_min=0.0, e_max=1.0):
+        """Initialize the uniform eccentricity distribution.
+
+        Parameters
+        ----------
+        e_min : float, optional
+            Minimum eccentricity (default: 0.0). Must be in [0, 1).
+        e_max : float, optional
+            Maximum eccentricity (default: 1.0). Must be in (0, 1].
+
+        Raises
+        ------
+        ValueError
+            If e_min or e_max are not in [0, 1], or if e_min >= e_max.
+        """
+        if not (0 <= e_min < 1):
+            raise ValueError("e_min must be in [0, 1)")
+        if not (0 < e_max <= 1):
+            raise ValueError("e_max must be in (0, 1]")
+        if e_min >= e_max:
+            raise ValueError("e_min must be less than e_max")
+
+        self.e_min = e_min
+        self.e_max = e_max
+        self.norm = 1.0 / (e_max - e_min)
+
+    def __repr__(self):
+        """Return string representation of the distribution."""
+        return f"UniformEccentricity(e_min={self.e_min}, e_max={self.e_max})"
+
+    def _repr_html_(self):
+        """Return HTML representation for Jupyter notebooks."""
+        return (f"<h3>Uniform Eccentricity Distribution</h3>"
+                f"<p>e_min = {self.e_min}</p>"
+                f"<p>e_max = {self.e_max}</p>")
+
+    def pdf(self, e):
+        """Probability density function of the uniform eccentricity distribution.
+
+        Parameters
+        ----------
+        e : float or array_like
+            Eccentricity value(s).
+
+        Returns
+        -------
+        float or ndarray
+            Probability density at eccentricity e.
+        """
+        e = np.asarray(e)
+        valid = (e >= self.e_min) & (e <= self.e_max)
+        pdf_values = np.zeros_like(e, dtype=float)
+        pdf_values[valid] = self.norm
+        return pdf_values
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the uniform eccentricity distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random eccentricity samples in [e_min, e_max].
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        return rng.uniform(self.e_min, self.e_max, size=size)
+
+
+class ZeroEccentricity:
+    """Zero eccentricity distribution for circular binary orbits.
+
+    All samples are exactly zero (circular orbits). The PDF is a Dirac delta
+    at e=0.
+
+    Examples
+    --------
+    >>> dist = ZeroEccentricity()
+    >>> e_samples = dist.rvs(size=1000)  # All zeros
+    """
+
+    def __init__(self):
+        """Initialize the zero eccentricity distribution."""
+        pass
+
+    def __repr__(self):
+        """Return string representation of the distribution."""
+        return "ZeroEccentricity()"
+
+    def _repr_html_(self):
+        """Return HTML representation for Jupyter notebooks."""
+        return "<h3>Zero Eccentricity Distribution</h3><p>e = 0 (circular orbits)</p>"
+
+    def pdf(self, e):
+        """Probability density function of the zero eccentricity distribution.
+
+        This is formally a Dirac delta at e=0. Returns 1.0 at e=0, 0.0 elsewhere.
+
+        Parameters
+        ----------
+        e : float or array_like
+            Eccentricity value(s).
+
+        Returns
+        -------
+        float or ndarray
+            1.0 where e==0, 0.0 elsewhere.
+        """
+        e = np.asarray(e)
+        return np.where(e == 0, 1.0, 0.0)
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the zero eccentricity distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator (unused, included for API consistency).
+
+        Returns
+        -------
+        ndarray
+            Array of zeros with shape (size,).
+        """
+        return np.zeros(size)
+
+
+class LogNormalSeparation:
+    """Log-normal orbital separation distribution for binary star systems.
+
+    Orbital separations are drawn from a log-normal distribution in log10 space,
+    truncated between specified minimum and maximum values.
+
+    Parameters
+    ----------
+    mean : float, optional
+        Mean of the log10 distribution (default: 3.5, corresponding to ~3162 Rsun).
+    sigma : float, optional
+        Standard deviation of the log10 distribution (default: 2.3).
+    min : float, optional
+        Minimum orbital separation in solar radii (default: 5.0).
+    max : float, optional
+        Maximum orbital separation in solar radii (default: 1e5).
+
+    Raises
+    ------
+    ValueError
+        If min is not positive, max <= min, or sigma <= 0.
+
+    Examples
+    --------
+    >>> dist = LogNormalSeparation(mean=3.0, sigma=1.5, min=10, max=1e6)
+    >>> separations = dist.rvs(size=1000)
+
+    Notes
+    -----
+    Uses scipy.stats.truncnorm for efficient sampling from the truncated
+    normal distribution in log10 space.
+    """
+
+    def __init__(self, mean=3.5, sigma=2.3, min=5.0, max=1e5):
+        """Initialize the log-normal separation distribution.
+
+        Parameters
+        ----------
+        mean : float, optional
+            Mean of the log10 distribution (default: 3.5).
+        sigma : float, optional
+            Standard deviation of the log10 distribution (default: 2.3).
+        min : float, optional
+            Minimum orbital separation in solar radii (default: 5.0).
+        max : float, optional
+            Maximum orbital separation in solar radii (default: 1e5).
+
+        Raises
+        ------
+        ValueError
+            If min is not positive, max <= min, or sigma <= 0.
+        """
+        if min <= 0:
+            raise ValueError("min must be positive")
+        if max <= min:
+            raise ValueError("max must be greater than min")
+        if sigma <= 0:
+            raise ValueError("sigma must be positive")
+
+        self.mean = mean
+        self.sigma = sigma
+        self.min = min
+        self.max = max
+
+        # Compute truncation bounds for scipy.stats.truncnorm
+        self.a_low = (np.log10(min) - mean) / sigma
+        self.a_high = (np.log10(max) - mean) / sigma
+
+    def __repr__(self):
+        """Return string representation of the distribution."""
+        return (f"LogNormalSeparation(mean={self.mean}, sigma={self.sigma}, "
+                f"min={self.min}, max={self.max})")
+
+    def _repr_html_(self):
+        """Return HTML representation for Jupyter notebooks."""
+        return (f"<h3>Log-Normal Separation Distribution</h3>"
+                f"<p>mean (log10) = {self.mean}</p>"
+                f"<p>sigma (log10) = {self.sigma}</p>"
+                f"<p>min = {self.min} R☉</p>"
+                f"<p>max = {self.max} R☉</p>")
+
+    def pdf(self, a):
+        """Probability density function of the log-normal separation distribution.
+
+        Parameters
+        ----------
+        a : float or array_like
+            Orbital separation(s) in solar radii.
+
+        Returns
+        -------
+        float or ndarray
+            Probability density at separation a.
+        """
+        a = np.asarray(a)
+        valid = (a > 0) & (a >= self.min) & (a <= self.max)
+        pdf_values = np.zeros_like(a, dtype=float)
+
+        log_a = np.log10(a[valid])
+        # PDF of truncnorm in log space, transformed to linear space
+        # pdf(a) = pdf_log(log a) * |d(log a)/da| = pdf_log(log a) / (a * ln(10))
+        pdf_values[valid] = truncnorm.pdf(
+            log_a,
+            self.a_low, self.a_high,
+            loc=self.mean,
+            scale=self.sigma
+        ) / (a[valid] * np.log(10))
+
+        return pdf_values
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the log-normal separation distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random orbital separation samples in solar radii.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Sample from truncated normal in log10 space
+        log_separations = truncnorm.rvs(
+            self.a_low, self.a_high,
+            loc=self.mean,
+            scale=self.sigma,
+            size=size,
+            random_state=rng
+        )
+
+        # Convert back to linear space
+        return 10**log_separations
diff --git a/posydon/popsyn/independent_sample.py b/posydon/popsyn/independent_sample.py
index 13e99063c1..f75eec9d5a 100644
--- a/posydon/popsyn/independent_sample.py
+++ b/posydon/popsyn/independent_sample.py
@@ -16,6 +16,7 @@
 import numpy as np
 from scipy.stats import truncnorm
 
+from posydon.popsyn import IMFs, distributions
 from posydon.popsyn.Moes_distributions import Moe_17_PsandQs
 from posydon.utils.common_functions import rejection_sampler
 
@@ -130,45 +131,11 @@ def generate_orbital_periods(primary_masses,
 
     # Sana H., et al., 2012, Science, 337, 444
     if orbital_period_scheme == 'Sana+12_period_extended':
-        # compute periods as if all M1 <= 15Msun (where pi = 0.0)
-        orbital_periods_M_lt_15 = 10**RNG.uniform(
-            low=np.log10(orbital_period_min),
-            high=np.log10(orbital_period_max),
-            size=number_of_binaries)
-
-        # compute periods as if all M1 > 15Msun
-        def pdf(logp):
-            pi = 0.55
-            beta = 1 - pi
-            A = np.log10(10**0.15)**(-pi)*(np.log10(10**0.15)
-                                           - np.log10(orbital_period_min))
-            B = 1./beta*(np.log10(orbital_period_max)**beta
-                         - np.log10(10**0.15)**beta)
-            C = 1./(A + B)
-            pdf = np.zeros(len(logp))
-
-            for j, logp_j in enumerate(logp):
-                # for logP<=0.15 days, the pdf is uniform
-                if np.log10(orbital_period_min) <= logp_j and logp_j < 0.15:
-                    pdf[j] = C*0.15**(-pi)
-
-                # original Sana H., et al., 2012, Science, 337, 444
-                elif 0.15 <= logp_j and logp_j < np.log10(orbital_period_max):
-                    pdf[j] = C*logp_j**(-pi)
-
-                else:
-                    pdf[j] = 0.
-
-            return pdf
-
-        orbital_periods_M_gt_15 = 10**(rejection_sampler(
-            size=number_of_binaries,
-            x_lim=[np.log10(orbital_period_min), np.log10(orbital_period_max)],
-            pdf=pdf))
-
-        orbital_periods = np.where(primary_masses <= 15.0,
-                                   orbital_periods_M_lt_15,
-                                   orbital_periods_M_gt_15)
+        period_dist = distributions.Sana12Period(
+            p_min=orbital_period_min,
+            p_max=orbital_period_max
+        )
+        orbital_periods = period_dist.rvs(size=number_of_binaries, m1=primary_masses, rng=RNG)
     else:
         raise ValueError("You must provide an allowed orbital period scheme.")
 
@@ -216,10 +183,11 @@ def generate_orbital_separations(number_of_binaries=1,
                          "orbital separation scheme.")
 
     if orbital_separation_scheme == 'log_uniform':
-        orbital_separations = 10**RNG.uniform(
-            low=np.log10(orbital_separation_min),
-            high=np.log10(orbital_separation_max),
-            size=number_of_binaries)
+        sep_dist = distributions.LogUniform(
+            min=orbital_separation_min,
+            max=orbital_separation_max
+        )
+        orbital_separations = sep_dist.rvs(size=number_of_binaries, rng=RNG)
 
     if orbital_separation_max < orbital_separation_min:
         raise ValueError("`orbital_separation_max` must be "
@@ -233,20 +201,13 @@ def generate_orbital_separations(number_of_binaries=1,
                 "`log_orbital_separation_mean`, "
                 "`log_orbital_separation_sigma`.")
 
-        # Set limits for truncated normal distribution
-        a_low = (np.log10(orbital_separation_min)
-                 - log_orbital_separation_mean) / log_orbital_separation_sigma
-        a_high = (np.log10(orbital_separation_max)
-                  - log_orbital_separation_mean) / log_orbital_separation_sigma
-
-        # generate orbital separations from a truncted normal distribution
-        log_orbital_separations = truncnorm.rvs(
-            a_low, a_high,
-            loc=log_orbital_separation_mean,
-            scale=log_orbital_separation_sigma,
-            size=number_of_binaries,
-            random_state=RNG)
-        orbital_separations = 10**log_orbital_separations
+        sep_dist = distributions.LogNormalSeparation(
+            mean=log_orbital_separation_mean,
+            sigma=log_orbital_separation_sigma,
+            min=orbital_separation_min,
+            max=orbital_separation_max
+        )
+        orbital_separations = sep_dist.rvs(size=number_of_binaries, rng=RNG)
 
     else:
         pass
@@ -284,11 +245,14 @@ def generate_eccentricities(number_of_binaries=1,
         raise ValueError("You must provide an allowed eccentricity scheme.")
 
     if eccentricity_scheme == 'thermal':
-        eccentricities = np.sqrt(RNG.uniform(size=number_of_binaries))
+        ecc_dist = distributions.ThermalEccentricity()
+        eccentricities = ecc_dist.rvs(size=number_of_binaries, rng=RNG)
     elif eccentricity_scheme == 'uniform':
-        eccentricities = RNG.uniform(size=number_of_binaries)
+        ecc_dist = distributions.UniformEccentricity()
+        eccentricities = ecc_dist.rvs(size=number_of_binaries, rng=RNG)
     elif eccentricity_scheme == 'zero':
-        eccentricities = np.zeros(number_of_binaries)
+        ecc_dist = distributions.ZeroEccentricity()
+        eccentricities = ecc_dist.rvs(size=number_of_binaries, rng=RNG)
     else:
         # This should never be reached
         pass
@@ -331,12 +295,8 @@ def generate_primary_masses(number_of_binaries=1,
 
     # Salpeter E. E., 1955, ApJ, 121, 161
     if primary_mass_scheme == 'Salpeter':
-        alpha = 2.35
-        normalization_constant = (1.0-alpha) / (primary_mass_max**(1-alpha)
-                                                - primary_mass_min**(1-alpha))
-        random_variable = RNG.uniform(size=number_of_binaries)
-        primary_masses = (random_variable*(1.0-alpha)/normalization_constant
-                          + primary_mass_min**(1.0-alpha))**(1.0/(1.0-alpha))
+        imf = IMFs.Salpeter(alpha=2.35, m_min=primary_mass_min, m_max=primary_mass_max)
+        primary_masses = imf.rvs(size=number_of_binaries, rng=RNG)
 
     # Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
     elif primary_mass_scheme == 'Kroupa1993':
@@ -349,12 +309,8 @@ def generate_primary_masses(number_of_binaries=1,
 
     # Kroupa P., 2001, MNRAS, 322, 231
     elif primary_mass_scheme == 'Kroupa2001':
-        alpha = 2.3
-        normalization_constant = (1.0-alpha) / (primary_mass_max**(1-alpha)
-                                                - primary_mass_min**(1-alpha))
-        random_variable = RNG.uniform(size=number_of_binaries)
-        primary_masses = (random_variable*(1.0-alpha)/normalization_constant
-                          + primary_mass_min**(1.0-alpha))**(1.0/(1.0-alpha))
+        imf = IMFs.Kroupa2001(m_min=primary_mass_min, m_max=primary_mass_max)
+        primary_masses = imf.rvs(size=number_of_binaries, rng=RNG)
     else:
         pass
 
@@ -404,13 +360,18 @@ def generate_secondary_masses(primary_masses,
 
     # Generate secondary masses
     if secondary_mass_scheme == 'flat_mass_ratio':
-        mass_ratio_min = np.max([secondary_mass_min / primary_masses,
-                                 np.ones(len(primary_masses))*0.05], axis=0)
-        mass_ratio_max = np.min([secondary_mass_max / primary_masses,
-                                 np.ones(len(primary_masses))], axis=0)
-        secondary_masses = (
-            (mass_ratio_max - mass_ratio_min) * RNG.uniform(
-                size=number_of_binaries) + mass_ratio_min) * primary_masses
+        # Calculate mass ratio bounds for each primary mass
+        q_min = np.maximum(secondary_mass_min / primary_masses, 0.05)
+        q_max = np.minimum(secondary_mass_max / primary_masses, 1.0)
+
+        # Sample mass ratios using the distribution class
+        # For mass-dependent bounds, we need to sample individually
+        mass_ratios = np.zeros(number_of_binaries)
+        for i in range(number_of_binaries):
+            q_dist = distributions.FlatMassRatio(q_min=q_min[i], q_max=q_max[i])
+            mass_ratios[i] = q_dist.rvs(size=1, rng=RNG)[0]
+
+        secondary_masses = primary_masses * mass_ratios
 
     if secondary_mass_scheme == 'q=1':
         secondary_masses = primary_masses
diff --git a/posydon/popsyn/norm_pop.py b/posydon/popsyn/norm_pop.py
index 4f2946504e..da195213a6 100644
--- a/posydon/popsyn/norm_pop.py
+++ b/posydon/popsyn/norm_pop.py
@@ -100,6 +100,7 @@ def get_pdf_for_m1(m1):
                 [kwargs['secondary_mass_max'] / m1, np.ones(len(m1))],
                 axis=0)
 
+            # Use FlatMassRatio distribution class
             q_dist = lambda q: np.where((q >= minimum) & (q <= maximum),
                                         1/(maximum - minimum),
                                         0)
@@ -107,16 +108,17 @@ def get_pdf_for_m1(m1):
         q_pdf = lambda q, m1: get_pdf_for_m1(m1)(q)
     elif kwargs['secondary_mass_scheme'] == 'flat_mass_ratio':
         # flat mass ratio, where bounds are given
-        q_pdf = lambda q, m1=None: np.where(
-                                (q > kwargs['q_min']) & (q <= kwargs['q_max']),
-                                1/(kwargs['q_max'] - kwargs['q_min']),
-                                0)
+        from posydon.popsyn.distributions import FlatMassRatio
+        q_dist = FlatMassRatio(q_min=kwargs['q_min'], q_max=kwargs['q_max'])
+        q_pdf = lambda q, m1=None: q_dist.pdf(q)
 
     else:
         # default to a flat distribution
         Pwarn("The secondary_mass_scheme is not defined use a flat mass ratio "
               "distribution in (0,1].", "UnsupportedModelWarning")
-        q_pdf = lambda q, m1=None: np.where((q > 0.0) & (q<=1.0), 1, 0)
+        from posydon.popsyn.distributions import FlatMassRatio
+        q_dist = FlatMassRatio(q_min=0.0, q_max=1.0)
+        q_pdf = lambda q, m1=None: q_dist.pdf(q)
     return q_pdf
 
 def get_binary_fraction_pdf(kwargs):
diff --git a/posydon/utils/common_functions.py b/posydon/utils/common_functions.py
index 9716a20051..93883bdcb6 100644
--- a/posydon/utils/common_functions.py
+++ b/posydon/utils/common_functions.py
@@ -669,7 +669,7 @@ def bondi_hoyle(binary, accretor, donor, idx=-1, wind_disk_criteria=True,
     return np.squeeze(mdot_acc)
 
 
-def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None):
+def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None, rng=None):
     """Generate a sample from a 1d PDF using the acceptance-rejection method.
 
     Parameters
@@ -684,6 +684,8 @@ def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None):
         The boundary where the pdf function is defined if passed as pdf.
     pdf : func
         The pdf function
+    rng : numpy.random.Generator, optional
+        Random number generator. If None, uses np.random.default_rng().
 
     Returns
     -------
@@ -691,6 +693,8 @@ def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None):
         An array with `size` random numbers generated from the PDF.
 
     """
+    if rng is None:
+        rng = np.random.default_rng()
     if pdf is None:
         if ((x is None) or (y is None)):
             raise ValueError("x and y PDF values must be specified if no PDF"
@@ -702,32 +706,32 @@ def rejection_sampler(x=None, y=None, size=1, x_lim=None, pdf=None):
             idxs = np.argsort(x)
             pdf = interp1d(x.take(idxs), y.take(idxs))
 
-        x_rand = np.random.uniform(x.min(), x.max(), size)
-        y_rand = np.random.uniform(0, y.max(), size)
+        x_rand = rng.uniform(x.min(), x.max(), size)
+        y_rand = rng.uniform(0, y.max(), size)
         values = x_rand[y_rand <= pdf(x_rand)]
         while values.shape[0] < size:
             n = size - values.shape[0]
-            x_rand = np.random.uniform(x.min(), x.max(), n)
-            y_rand = np.random.uniform(0, y.max(), n)
+            x_rand = rng.uniform(x.min(), x.max(), n)
+            y_rand = rng.uniform(0, y.max(), n)
             values = np.hstack([values, x_rand[y_rand <= pdf(x_rand)]])
     else:
         if x_lim is None:
             raise ValueError("x_lim must be specified for passed PDF function"
                              " in rejection sampling")
-        x_rand = np.random.uniform(x_lim[0], x_lim[1], size)
-        pdf_max = max(pdf(np.random.uniform(x_lim[0], x_lim[1], 50000)))
-        y_rand = np.random.uniform(0, pdf_max, size)
+        x_rand = rng.uniform(x_lim[0], x_lim[1], size)
+        pdf_max = max(pdf(rng.uniform(x_lim[0], x_lim[1], 50000)))
+        y_rand = rng.uniform(0, pdf_max, size)
         values = x_rand[y_rand <= pdf(x_rand)]
         while values.shape[0] < size:
             n = size - values.shape[0]
-            x_rand = np.random.uniform(x_lim[0], x_lim[1], n)
-            y_rand = np.random.uniform(0, pdf_max, n)
+            x_rand = rng.uniform(x_lim[0], x_lim[1], n)
+            y_rand = rng.uniform(0, pdf_max, n)
             values = np.hstack([values, x_rand[y_rand <= pdf(x_rand)]])
 
     return values
 
 
-def inverse_sampler(x, y, size=1):
+def inverse_sampler(x, y, size=1, rng=None):
     """Sample from a PDF using the inverse-sampling method.
 
     Parameters
@@ -738,6 +742,8 @@ def inverse_sampler(x, y, size=1):
         The probablity density at `x` (or a scaled version of it).
     size : int
         Number of samples to generate.
+    rng : numpy.random.Generator, optional
+        Random number generator. If None, uses np.random.default_rng().
 
     Returns
     -------
@@ -745,6 +751,8 @@ def inverse_sampler(x, y, size=1):
         The sample drawn from the PDF.
 
     """
+    if rng is None:
+        rng = np.random.default_rng()
     # x should be sorted
     assert np.all(np.diff(x) > 0)
     # y should be above 0
@@ -757,7 +765,7 @@ def inverse_sampler(x, y, size=1):
     total_area = cumsum_areas[-1]
 
     # start the inverse sampling
-    u = np.random.uniform(size=size)
+    u = rng.uniform(size=size)
     # index of "bin" where each sampled value corresponds too
     u_indices = np.searchsorted(cumsum_areas, u * total_area)
     # the area that should be covered from the end of the previous bin
@@ -779,7 +787,7 @@ def inverse_sampler(x, y, size=1):
     if n_where_nan:
         assert np.all(dy_bins[where_nan] == 0)
         sample[where_nan] = x[u_indices][where_nan] + (
-            dx_bins[where_nan] * np.random.uniform(size=n_where_nan))
+            dx_bins[where_nan] * rng.uniform(size=n_where_nan))
 
     # make sure that everything worked as expected
     assert np.all(np.isfinite(sample))

From e948078f8da333350aac78743b1bc65569e01dce Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 09:09:02 +0100
Subject: [PATCH 02/18] add Kroupa1993 IMF

---
 posydon/popsyn/IMFs.py               | 98 ++++++++++++++++++++++++++++
 posydon/popsyn/independent_sample.py |  8 +--
 2 files changed, 100 insertions(+), 6 deletions(-)

diff --git a/posydon/popsyn/IMFs.py b/posydon/popsyn/IMFs.py
index 62055864fb..6b76096763 100644
--- a/posydon/popsyn/IMFs.py
+++ b/posydon/popsyn/IMFs.py
@@ -199,6 +199,104 @@ def rvs(self, size=1, rng=None):
         return masses
 
 
+class Kroupa1993(IMFBase):
+    """
+    Initial Mass Function based on Kroupa et al. (1993), which is defined as:
+
+        dN/dM = m^-2.7
+
+    References
+    ----------
+    Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
+    https://ui.adsabs.harvard.edu/abs/1993MNRAS.262..545K/abstract
+
+    Parameters
+    ----------
+    alpha : float, optional
+        The power-law index of the IMF (default is 2.7).
+    m_min : float, optional
+        The minimum allowable mass (default is 0.01) [Msun].
+    m_max : float, optional
+        The maximum allowable mass (default is 200.0) [Msun].
+
+    Attributes
+    ----------
+    alpha : float
+        Power-law index used in the IMF calculation.
+    m_min : float
+        Minimum stellar mass for the IMF [Msun].
+    m_max : float
+        Maximum stellar mass for the IMF [Msun].
+    """
+
+    def __init__(self, alpha=2.7, m_min=0.01, m_max=200.0):
+        self.alpha = alpha
+        super().__init__(m_min, m_max)
+
+    def __repr__(self):
+        return (f"Kroupa1993(alpha={self.alpha}, "
+                f"m_min={self.m_min}, "
+                f"m_max={self.m_max})")
+
+    def _repr_html_(self):
+        return (f"<h3>Kroupa (1993) IMF</h3>"
+                f"<p>alpha = {self.alpha}</p>"
+                f"<p>m_min = {self.m_min}</p>"
+                f"<p>m_max = {self.m_max}</p>")
+
+    def imf(self, m):
+        '''Computes the IMF value for a given mass or array of masses 'm'.
+        Raises a ValueError if any value in 'm' is less than or equal to zero.
+
+        Parameters
+        ----------
+        m : float or array_like
+            Stellar mass or array of stellar masses [Msun].
+
+        Returns
+        -------
+        float or ndarray
+            The IMF value for the given mass or masses.
+        '''
+        m = np.asarray(m)
+        if np.any(m <= 0):
+            raise ValueError("Mass must be positive.")
+        valid = self._check_valid(m)
+        return m ** (-self.alpha)
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the Kroupa1993 IMF.
+
+        Uses analytical inverse transform sampling for efficiency.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random mass samples in solar masses.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Analytical inverse transform sampling
+        # For power law IMF: m^(-alpha), the inverse CDF is:
+        # m = (u * (m_max^(1-alpha) - m_min^(1-alpha)) + m_min^(1-alpha))^(1/(1-alpha))
+        normalization_constant = (1.0 - self.alpha) / (
+            self.m_max**(1.0 - self.alpha) - self.m_min**(1.0 - self.alpha)
+        )
+        u = rng.uniform(size=size)
+        masses = (u * (1.0 - self.alpha) / normalization_constant
+                  + self.m_min**(1.0 - self.alpha))**(1.0 / (1.0 - self.alpha))
+
+        return masses
+
+
 class Kroupa2001(IMFBase):
     """Initial Mass Function based on Kroupa (2001), which is
     defined as a broken power-law:
diff --git a/posydon/popsyn/independent_sample.py b/posydon/popsyn/independent_sample.py
index f75eec9d5a..95d9c3c909 100644
--- a/posydon/popsyn/independent_sample.py
+++ b/posydon/popsyn/independent_sample.py
@@ -300,12 +300,8 @@ def generate_primary_masses(number_of_binaries=1,
 
     # Kroupa P., Tout C. A., Gilmore G., 1993, MNRAS, 262, 545
     elif primary_mass_scheme == 'Kroupa1993':
-        alpha = 2.7
-        normalization_constant = (1.0-alpha) / (primary_mass_max**(1-alpha)
-                                                - primary_mass_min**(1-alpha))
-        random_variable = RNG.uniform(size=number_of_binaries)
-        primary_masses = (random_variable*(1.0-alpha)/normalization_constant
-                          + primary_mass_min**(1.0-alpha))**(1.0/(1.0-alpha))
+        imf = IMFs.Kroupa1993(alpha=2.7, m_min=primary_mass_min, m_max=primary_mass_max)
+        primary_masses = imf.rvs(size=number_of_binaries, rng=RNG)
 
     # Kroupa P., 2001, MNRAS, 322, 231
     elif primary_mass_scheme == 'Kroupa2001':

From 419eda39efe63f8c2a148e7013dfb739b19440a7 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 10:46:08 +0100
Subject: [PATCH 03/18] fix q boundaries

---
 posydon/popsyn/distributions.py                 | 10 +++++-----
 posydon/unit_tests/popsyn/test_distributions.py | 14 +++++---------
 2 files changed, 10 insertions(+), 14 deletions(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index 82e721f804..aec4b4aa55 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -39,19 +39,19 @@ def __init__(self, q_min=0.05, q_max=1):
         Parameters
         ----------
         q_min : float, optional
-            Minimum mass ratio (default: 0.05). Must be in (0, 1].
+            Minimum mass ratio (default: 0.05). Must be in [0, 1).
         q_max : float, optional
             Maximum mass ratio (default: 1.0). Must be in (0, 1].
 
         Raises
         ------
         ValueError
-            If q_min or q_max are not in (0, 1], or if q_min >= q_max.
+            If q_min or q_max are not in valid range, or if q_min >= q_max.
         """
-        if not (0 < q_min <= 1):
-            raise ValueError("q_min must be in (0, 1)")
+        if not (0 <= q_min < 1):
+            raise ValueError("q_min must be in [0, 1)")
         if not (0 < q_max <= 1):
-            raise ValueError("q_max must be in (0, 1)")
+            raise ValueError("q_max must be in (0, 1]")
         if q_min >= q_max:
             raise ValueError("q_min must be less than q_max")
 
diff --git a/posydon/unit_tests/popsyn/test_distributions.py b/posydon/unit_tests/popsyn/test_distributions.py
index c2f80cbedf..b781209bd7 100644
--- a/posydon/unit_tests/popsyn/test_distributions.py
+++ b/posydon/unit_tests/popsyn/test_distributions.py
@@ -45,24 +45,20 @@ def test_initialization_custom(self, custom_flat_ratio):
 
     def test_initialization_invalid_parameters(self):
         """Test that initialization raises ValueError for invalid parameters."""
-        # Test q_min not in (0, 1]
-        with pytest.raises(ValueError, match="q_min must be in \\(0, 1\\)"):
-            FlatMassRatio(q_min=0.0, q_max=0.5)
-
-        with pytest.raises(ValueError, match="q_min must be in \\(0, 1\\)"):
+        with pytest.raises(ValueError, match="q_min must be in \\[0, 1\\)"):
             FlatMassRatio(q_min=-0.1, q_max=0.5)
 
-        with pytest.raises(ValueError, match="q_min must be in \\(0, 1\\)"):
+        with pytest.raises(ValueError, match="q_min must be in \\[0, 1\\)"):
             FlatMassRatio(q_min=1.5, q_max=2.0)
 
         # Test q_max not in (0, 1]
-        with pytest.raises(ValueError, match="q_max must be in \\(0, 1\\)"):
+        with pytest.raises(ValueError, match="q_max must be in \\(0, 1\\]"):
             FlatMassRatio(q_min=0.1, q_max=0.0)
 
-        with pytest.raises(ValueError, match="q_max must be in \\(0, 1\\)"):
+        with pytest.raises(ValueError, match="q_max must be in \\(0, 1\\]"):
             FlatMassRatio(q_min=0.1, q_max=-0.1)
 
-        with pytest.raises(ValueError, match="q_max must be in \\(0, 1\\)"):
+        with pytest.raises(ValueError, match="q_max must be in \\(0, 1\\]"):
             FlatMassRatio(q_min=0.1, q_max=1.5)
 
         # Test q_min >= q_max

From 6c0be9d76e44265fe95d61678e13f4cd6186aeae Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 10:52:20 +0100
Subject: [PATCH 04/18] revert back to exlusion of q=0 in output

---
 posydon/popsyn/distributions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index aec4b4aa55..ff0fdfd53b 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -125,7 +125,7 @@ def pdf(self, q):
             Probability density at mass ratio q.
         """
         q = np.asarray(q)
-        valid = (q >= self.q_min) & (q <= self.q_max)
+        valid = (q > self.q_min) & (q <= self.q_max)
         pdf_values = np.zeros_like(q, dtype=float)
         pdf_values[valid] = self.flat_mass_ratio(q[valid]) * self.norm
         return pdf_values

From 4d2915d887a52578cce95792235a809db3526b9f Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 10:59:24 +0100
Subject: [PATCH 05/18] update unit tests to use a mock RNG class instead of
 patching np.uniform

---
 .../unit_tests/utils/test_common_functions.py | 29 ++++++++++---------
 1 file changed, 16 insertions(+), 13 deletions(-)

diff --git a/posydon/unit_tests/utils/test_common_functions.py b/posydon/unit_tests/utils/test_common_functions.py
index ef405107b6..44786badf6 100644
--- a/posydon/unit_tests/utils/test_common_functions.py
+++ b/posydon/unit_tests/utils/test_common_functions.py
@@ -754,8 +754,10 @@ def mock_rand2(shape):
                approx(5.13970075150e-8, abs=6e-20)
 
     def test_rejection_sampler(self, monkeypatch):
-        def mock_uniform(low=0.0, high=1.0, size=1):
-            return np.linspace(low, high, num=size)
+        class MockRNG:
+            def uniform(self, low=0.0, high=1.0, size=1):
+                return np.linspace(low, high, num=size)
+        mock_rng = MockRNG()
         def mock_interp1d(x, y):
             if x[0] > x[-1]:
                 raise ValueError
@@ -786,7 +788,6 @@ def mock_pdf(x):
                                       +" function in rejection sampling"):
             totest.rejection_sampler(pdf=mock_pdf)
         # examples:
-        monkeypatch.setattr(np.random, "uniform", mock_uniform)
         monkeypatch.setattr(totest, "interp1d", mock_interp1d)
         tests = [(np.array([0.0, 1.0]), np.array([0.4, 0.6]), 5,\
                   np.array([0.0, 0.25, 0.5, 0.75, 1.0])),\
@@ -795,22 +796,24 @@ def mock_pdf(x):
                  (np.array([1.0, 0.0]), np.array([0.2, 0.8]), 6,\
                   np.array([0.0, 0.2, 0.4, 0.0, 0.5, 0.0]))]
         for (x, y, s, r) in tests:
-            assert np.allclose(totest.rejection_sampler(x=x, y=y, size=s),\
+            assert np.allclose(totest.rejection_sampler(x=x, y=y, size=s, rng=mock_rng),\
                                r)
         assert np.allclose(totest.rejection_sampler(x_lim=np.array([0.0,\
                                                                      1.0]),\
-                                                     pdf=mock_pdf, size=5),\
+                                                     pdf=mock_pdf, size=5, rng=mock_rng),\
                            np.array([0.0, 0.25, 0.0, 0.0, 0.0]))
         assert np.allclose(totest.rejection_sampler(x=np.array([1.0, 0.0]),\
                                                        y=np.array([0.2, 0.8]),\
                                                        x_lim=np.array([0.0,\
                                                                        1.0]),\
-                                                       pdf=mock_pdf, size=5),\
+                                                       pdf=mock_pdf, size=5, rng=mock_rng),\
                               np.array([0.0, 0.25, 0.0, 0.0, 0.0]))
 
     def test_inverse_sampler(self, monkeypatch):
-        def mock_uniform(low=0.0, high=1.0, size=1):
-            return np.linspace(low, high, num=size)
+        class MockRNG:
+            def uniform(self, low=0.0, high=1.0, size=1):
+                return np.linspace(low, high, num=size)
+        mock_rng = MockRNG()
         # missing argument
         with raises(TypeError, match="missing 2 required positional "\
                                      +"arguments: 'x' and 'y'"):
@@ -826,21 +829,21 @@ def mock_uniform(low=0.0, high=1.0, size=1):
             totest.inverse_sampler(x=np.array([0.0, 1.0]),\
                                    y=np.array([-0.4, 0.6]))
         # examples:
-        monkeypatch.setattr(np.random, "uniform", mock_uniform)
-        assert np.allclose(totest.inverse_sampler(x=np.array([0.0, 1.0]),\
+        output = totest.inverse_sampler(x=np.array([0.0, 1.0]),\
                                                   y=np.array([0.4, 0.6]),\
-                                                  size=5),\
+                                                  size=5, rng=mock_rng)
+        assert np.allclose(output,\
                            np.array([0.0, 0.29128785, 0.54950976, 0.78388218,\
                                      1.0]))
         assert np.allclose(totest.inverse_sampler(x=np.array([0.0, 1.0]),\
                                                   y=np.array([0.6, 0.4]),\
-                                                  size=4),\
+                                                  size=4, rng=mock_rng),\
                            np.array([0.0, 0.2919872, 0.61952386, 1.0]))
         with warns(RuntimeWarning, match="invalid value encountered in "\
                                          +"divide"):
             assert np.allclose(totest.inverse_sampler(x=np.array([0.0, 1.0]),\
                                                       y=np.array([0.5, 0.5]),\
-                                                      size=5),\
+                                                      size=5, rng=mock_rng),\
                                np.array([0.0, 0.25, 0.5, 0.75, 1.0]))
 
     def test_histogram_sampler(self, monkeypatch):

From feca7b76b2396f7ce35ba7540a05ea14b67dc94d Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:03:38 +0100
Subject: [PATCH 06/18] set q_min inclusive again

---
 posydon/popsyn/distributions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index ff0fdfd53b..aec4b4aa55 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -125,7 +125,7 @@ def pdf(self, q):
             Probability density at mass ratio q.
         """
         q = np.asarray(q)
-        valid = (q > self.q_min) & (q <= self.q_max)
+        valid = (q >= self.q_min) & (q <= self.q_max)
         pdf_values = np.zeros_like(q, dtype=float)
         pdf_values[valid] = self.flat_mass_ratio(q[valid]) * self.norm
         return pdf_values

From f220e10a87ce4389a9ac6fd86fee5becf2f7c108 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:11:58 +0100
Subject: [PATCH 07/18] fix q boundary

---
 posydon/popsyn/distributions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index aec4b4aa55..ff0fdfd53b 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -125,7 +125,7 @@ def pdf(self, q):
             Probability density at mass ratio q.
         """
         q = np.asarray(q)
-        valid = (q >= self.q_min) & (q <= self.q_max)
+        valid = (q > self.q_min) & (q <= self.q_max)
         pdf_values = np.zeros_like(q, dtype=float)
         pdf_values[valid] = self.flat_mass_ratio(q[valid]) * self.norm
         return pdf_values

From f807624efd439dc2e3bb92c7b552d741cc81773f Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:15:59 +0100
Subject: [PATCH 08/18] exlude lower bound

---
 posydon/unit_tests/popsyn/test_distributions.py | 1 +
 1 file changed, 1 insertion(+)

diff --git a/posydon/unit_tests/popsyn/test_distributions.py b/posydon/unit_tests/popsyn/test_distributions.py
index b781209bd7..3e48bf889c 100644
--- a/posydon/unit_tests/popsyn/test_distributions.py
+++ b/posydon/unit_tests/popsyn/test_distributions.py
@@ -110,6 +110,7 @@ def test_pdf_within_range(self, custom_flat_ratio):
         q_values = np.linspace(custom_flat_ratio.q_min, custom_flat_ratio.q_max, 10)
         pdf_values = custom_flat_ratio.pdf(q_values)
         expected_pdf = custom_flat_ratio.norm * np.ones_like(q_values)
+        expected_pdf[0] = 0.0  # q_values[0] is equal to q_min, which is outside the valid range
         np.testing.assert_allclose(pdf_values, expected_pdf)
 
     def test_pdf_outside_range(self, custom_flat_ratio):

From 531d4829210af659bb6a8183131b36895b408873 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:48:50 +0100
Subject: [PATCH 09/18] add rvs unit tests + add Chabrier2003 rvs function

---
 posydon/popsyn/IMFs.py                        |  40 +-
 posydon/unit_tests/popsyn/test_IMFs.py        | 179 ++++++
 .../unit_tests/popsyn/test_distributions.py   | 563 ++++++++++++++++++
 posydon/unit_tests/popsyn/test_norm_pop.py    |  26 +
 4 files changed, 805 insertions(+), 3 deletions(-)

diff --git a/posydon/popsyn/IMFs.py b/posydon/popsyn/IMFs.py
index 6b76096763..f4467b168c 100644
--- a/posydon/popsyn/IMFs.py
+++ b/posydon/popsyn/IMFs.py
@@ -9,6 +9,7 @@
 import numpy as np
 from scipy.integrate import quad
 
+from posydon.utils.common_functions import inverse_sampler, rejection_sampler
 from posydon.utils.posydonwarning import Pwarn
 
 
@@ -101,6 +102,10 @@ def imf(self, m): # pragma: no cover
         '''
         pass
 
+    def rvs(self, size=1, rng=None):
+        pass
+
+
 class Salpeter(IMFBase):
     """
     Initial Mass Function based on Salpeter (1955), which is defined as:
@@ -427,9 +432,6 @@ def rvs(self, size=1, rng=None):
         if rng is None:
             rng = np.random.default_rng()
 
-        # Import here to avoid circular dependency
-        from posydon.utils.common_functions import inverse_sampler
-
         # Create discretized PDF for inverse sampling
         # Use more points near the breaks for better accuracy
         n_points = 2000
@@ -529,3 +531,35 @@ def imf(self, m):
         C = (1.0 / (self.m_break * sqrt_2pi_sigma)) * np.exp(-log_term_break)
         powerlaw = C * (m / self.m_break) ** (-self.alpha)
         return np.where(m < self.m_break, lognormal, powerlaw)
+
+    def rvs(self, size=1, rng=None):
+        """Draw random samples from the Chabrier2003 IMF.
+
+        Uses inverse transform sampling with discretized PDF for the
+        combined lognormal and power-law distribution.
+
+        Parameters
+        ----------
+        size : int, optional
+            Number of samples to draw (default: 1).
+        rng : numpy.random.Generator, optional
+            Random number generator. If None, uses np.random.default_rng().
+
+        Returns
+        -------
+        ndarray
+            Random mass samples in solar masses.
+        """
+        if rng is None:
+            rng = np.random.default_rng()
+
+        # Create discretized PDF for inverse sampling
+        # Use more points near the transition for better accuracy
+        n_points = 2000
+        m_grid = np.linspace(self.m_min, self.m_max, n_points)
+        pdf_values = self.imf(m_grid)
+
+        # Sample using inverse transform method
+        masses = inverse_sampler(m_grid, pdf_values, size=size, rng=rng)
+
+        return masses
diff --git a/posydon/unit_tests/popsyn/test_IMFs.py b/posydon/unit_tests/popsyn/test_IMFs.py
index 40e29c4839..485268bab9 100644
--- a/posydon/unit_tests/popsyn/test_IMFs.py
+++ b/posydon/unit_tests/popsyn/test_IMFs.py
@@ -158,6 +158,141 @@ def test_imf_outside_range_warning(self, default_imf):
             # Verify the function still returns values (just with warning)
             assert len(result) == 1
 
+    def test_rvs_default(self, default_imf):
+        """Test random sampling from the Salpeter IMF."""
+        # Test single sample
+        rng = np.random.default_rng(42)
+        sample = default_imf.rvs(size=1, rng=rng)
+        assert len(sample) == 1
+        assert default_imf.m_min <= sample[0] <= default_imf.m_max
+
+        # Test multiple samples
+        rng = np.random.default_rng(42)
+        samples = default_imf.rvs(size=1000, rng=rng)
+        assert len(samples) == 1000
+        assert np.all(samples >= default_imf.m_min)
+        assert np.all(samples <= default_imf.m_max)
+
+    def test_rvs_without_rng(self, default_imf):
+        """Test random sampling without providing an RNG."""
+        samples = default_imf.rvs(size=100)
+        assert len(samples) == 100
+        assert np.all(samples >= default_imf.m_min)
+        assert np.all(samples <= default_imf.m_max)
+
+class TestKroupa1993IMF:
+    @pytest.fixture
+    def default_kroupa1993(self):
+        """Fixture for default Kroupa1993 instance."""
+        return IMFs.Kroupa1993()
+
+    @pytest.fixture
+    def custom_kroupa1993(self):
+        """Fixture for custom Kroupa1993 instance."""
+        return IMFs.Kroupa1993(alpha=2.5, m_min=0.05, m_max=150.0)
+
+    def test_initialization_default(self, default_kroupa1993):
+        """Test default initialization of Kroupa1993 IMF."""
+        assert default_kroupa1993.alpha == 2.7
+        assert default_kroupa1993.m_min == 0.01
+        assert default_kroupa1993.m_max == 200.0
+        integral, _ = quad(default_kroupa1993.imf,
+                           default_kroupa1993.m_min,
+                           default_kroupa1993.m_max)
+        assert np.isclose(integral*default_kroupa1993.norm, 1.0, rtol=1e-5)
+
+    def test_initialization_custom(self, custom_kroupa1993):
+        """Test custom initialization of Kroupa1993 IMF."""
+        assert custom_kroupa1993.alpha == 2.5
+        assert custom_kroupa1993.m_min == 0.05
+        assert custom_kroupa1993.m_max == 150.0
+        integral, _ = quad(custom_kroupa1993.imf,
+                           custom_kroupa1993.m_min,
+                           custom_kroupa1993.m_max)
+        assert np.isclose(integral*custom_kroupa1993.norm, 1.0, rtol=1e-5)
+
+    def test_repr(self, default_kroupa1993):
+        """Test string representation."""
+        rep_str = default_kroupa1993.__repr__()
+        assert "Kroupa1993(" in rep_str
+        assert "alpha=2.7" in rep_str
+        assert "m_min=0.01" in rep_str
+        assert "m_max=200.0" in rep_str
+
+    def test_repr_html(self, default_kroupa1993):
+        """Test HTML representation."""
+        html_str = default_kroupa1993._repr_html_()
+        assert "<h3>Kroupa (1993) IMF</h3>" in html_str
+        assert "alpha = 2.7" in html_str
+        assert "m_min = 0.01" in html_str
+        assert "m_max = 200.0" in html_str
+
+    def test_invalid_mass(self, default_kroupa1993):
+        """Test that the imf method raises ValueError for invalid mass values."""
+        with pytest.raises(ValueError, match="Mass must be positive."):
+            default_kroupa1993.imf(0)
+        with pytest.raises(ValueError, match="Mass must be positive."):
+            default_kroupa1993.imf(-1.0)
+        with pytest.raises(ValueError, match="Mass must be positive."):
+            default_kroupa1993.imf([0.0])
+
+    def test_imf(self, default_kroupa1993):
+        """Test the imf method for correct values."""
+        # Test with an array of mass values
+        m_values = np.array([1.0, 2.0, 5.0, 10.0])
+        expected = m_values ** (-default_kroupa1993.alpha)
+        computed = default_kroupa1993.imf(m_values)
+        assert np.allclose(computed, expected)
+
+        # Test with a single mass value
+        m = 3.0
+        expected = m ** (-default_kroupa1993.alpha)
+        computed = default_kroupa1993.imf(m)
+        assert np.isclose(computed, expected)
+
+    def test_pdf_within_range(self, default_kroupa1993):
+        """Test that PDF returns correct values within the mass range."""
+        m = np.linspace(default_kroupa1993.m_min, default_kroupa1993.m_max, 100)
+        pdf_values = default_kroupa1993.pdf(m)
+        expected_pdf = default_kroupa1993.imf(m) * default_kroupa1993.norm
+        assert np.allclose(pdf_values, expected_pdf)
+
+    def test_pdf_outside_range(self, default_kroupa1993):
+        """Test that PDF returns zero for masses outside the mass range."""
+        m = np.array([default_kroupa1993.m_min - 0.1, default_kroupa1993.m_max + 0.1])
+        pdf_values = default_kroupa1993.pdf(m)
+        assert np.allclose(pdf_values, 0.0)
+
+    def test_normalization(self, default_kroupa1993):
+        """Ensure that the integral of the PDF over the range is approximately 1."""
+        integral, _ = quad(default_kroupa1993.pdf,
+                           default_kroupa1993.m_min,
+                           default_kroupa1993.m_max)
+        assert np.isclose(integral, 1.0, rtol=1e-5)
+
+    def test_rvs_default(self, default_kroupa1993):
+        """Test random sampling from the Kroupa1993 IMF."""
+        # Test single sample
+        rng = np.random.default_rng(42)
+        sample = default_kroupa1993.rvs(size=1, rng=rng)
+        assert len(sample) == 1
+        assert default_kroupa1993.m_min <= sample[0] <= default_kroupa1993.m_max
+
+        # Test multiple samples
+        rng = np.random.default_rng(42)
+        samples = default_kroupa1993.rvs(size=1000, rng=rng)
+        assert len(samples) == 1000
+        assert np.all(samples >= default_kroupa1993.m_min)
+        assert np.all(samples <= default_kroupa1993.m_max)
+
+    def test_rvs_without_rng(self, default_kroupa1993):
+        """Test random sampling without providing an RNG."""
+        samples = default_kroupa1993.rvs(size=100)
+        assert len(samples) == 100
+        assert np.all(samples >= default_kroupa1993.m_min)
+        assert np.all(samples <= default_kroupa1993.m_max)
+
+
 class TestKroupa2001IMF:
     @pytest.fixture
     def default_kroupa(self):
@@ -277,6 +412,28 @@ def test_repr_html(self, default_kroupa):
         html_str = default_kroupa._repr_html_()
         assert "<h3>Kroupa (2001) IMF</h3>" in html_str
 
+    def test_rvs_default(self, default_kroupa):
+        """Test random sampling from the Kroupa2001 IMF."""
+        # Test single sample
+        rng = np.random.default_rng(42)
+        sample = default_kroupa.rvs(size=1, rng=rng)
+        assert len(sample) == 1
+        assert default_kroupa.m_min <= sample[0] <= default_kroupa.m_max
+
+        # Test multiple samples
+        rng = np.random.default_rng(42)
+        samples = default_kroupa.rvs(size=1000, rng=rng)
+        assert len(samples) == 1000
+        assert np.all(samples >= default_kroupa.m_min)
+        assert np.all(samples <= default_kroupa.m_max)
+
+    def test_rvs_without_rng(self, default_kroupa):
+        """Test random sampling without providing an RNG."""
+        samples = default_kroupa.rvs(size=100)
+        assert len(samples) == 100
+        assert np.all(samples >= default_kroupa.m_min)
+        assert np.all(samples <= default_kroupa.m_max)
+
 class TestChabrierIMF:
     @pytest.fixture
     def default_chabrier(self):
@@ -401,3 +558,25 @@ def test_repr(self, default_chabrier):
     def test_repr_html(self, default_chabrier):
         html_str = default_chabrier._repr_html_()
         assert "<h3>Chabrier IMF</h3>" in html_str
+
+    def test_rvs_default(self, default_chabrier):
+        """Test random sampling from the Chabrier2003 IMF."""
+        # Test single sample
+        rng = np.random.default_rng(42)
+        sample = default_chabrier.rvs(size=1, rng=rng)
+        assert len(sample) == 1
+        assert default_chabrier.m_min <= sample[0] <= default_chabrier.m_max
+
+        # Test multiple samples
+        rng = np.random.default_rng(42)
+        samples = default_chabrier.rvs(size=1000, rng=rng)
+        assert len(samples) == 1000
+        assert np.all(samples >= default_chabrier.m_min)
+        assert np.all(samples <= default_chabrier.m_max)
+
+    def test_rvs_without_rng(self, default_chabrier):
+        """Test random sampling without providing an RNG."""
+        samples = default_chabrier.rvs(size=100)
+        assert len(samples) == 100
+        assert np.all(samples >= default_chabrier.m_min)
+        assert np.all(samples <= default_chabrier.m_max)
diff --git a/posydon/unit_tests/popsyn/test_distributions.py b/posydon/unit_tests/popsyn/test_distributions.py
index 3e48bf889c..23b787461e 100644
--- a/posydon/unit_tests/popsyn/test_distributions.py
+++ b/posydon/unit_tests/popsyn/test_distributions.py
@@ -10,9 +10,13 @@
 
 from posydon.popsyn.distributions import (
     FlatMassRatio,
+    LogNormalSeparation,
     LogUniform,
     PowerLawPeriod,
     Sana12Period,
+    ThermalEccentricity,
+    UniformEccentricity,
+    ZeroEccentricity,
 )
 
 
@@ -150,6 +154,23 @@ def test_normalization_integral(self, default_flat_ratio):
         integral, _ = quad(default_flat_ratio.pdf, default_flat_ratio.q_min, default_flat_ratio.q_max)
         np.testing.assert_allclose(integral, 1.0, rtol=1e-10)
 
+    def test_rvs(self, custom_flat_ratio):
+        """Test random sampling."""
+        rng = np.random.default_rng(42)
+        samples = custom_flat_ratio.rvs(size=1000, rng=rng)
+
+        assert len(samples) == 1000
+        assert np.all(samples >= custom_flat_ratio.q_min)
+        assert np.all(samples <= custom_flat_ratio.q_max)
+
+    def test_rvs_without_rng(self, custom_flat_ratio):
+        """Test random sampling without providing an RNG."""
+        samples = custom_flat_ratio.rvs(size=100)
+
+        assert len(samples) == 100
+        assert np.all(samples >= custom_flat_ratio.q_min)
+        assert np.all(samples <= custom_flat_ratio.q_max)
+
 
 class TestSana12Period:
     """Test class for Sana12Period distribution."""
@@ -314,6 +335,64 @@ def test_calculate_normalization_integration(self, default_sana12):
         norm_high = default_sana12._calculate_normalization(m1_high)
         assert norm_high > 0
 
+    def test_rvs_with_m1_none(self, default_sana12):
+        """Test that rvs raises ValueError when m1 is None."""
+        rng = np.random.default_rng(42)
+
+        with pytest.raises(ValueError, match="m1 \\(primary mass\\) must be provided"):
+            default_sana12.rvs(size=10, m1=None, rng=rng)
+
+    def test_rvs_with_m1_wrong_size(self, default_sana12):
+        """Test that rvs raises ValueError when m1 has wrong size."""
+        rng = np.random.default_rng(42)
+        m1_wrong_size = np.array([10.0, 15.0])
+
+        with pytest.raises(ValueError, match="m1 must be a single value or have size="):
+            default_sana12.rvs(size=10, m1=m1_wrong_size, rng=rng)
+
+    def test_rvs_low_mass(self, default_sana12):
+        """Test random sampling for low mass stars."""
+        rng = np.random.default_rng(42)
+        m1 = 10.0  # Below mbreak
+
+        samples = default_sana12.rvs(size=100, m1=m1, rng=rng)
+
+        assert len(samples) == 100
+        assert np.all(samples >= default_sana12.p_min)
+        assert np.all(samples <= default_sana12.p_max)
+
+    def test_rvs_high_mass(self, default_sana12):
+        """Test random sampling for high mass stars."""
+        rng = np.random.default_rng(42)
+        m1 = 25.0  # Above mbreak
+
+        samples = default_sana12.rvs(size=100, m1=m1, rng=rng)
+
+        assert len(samples) == 100
+        assert np.all(samples >= default_sana12.p_min)
+        assert np.all(samples <= default_sana12.p_max)
+
+    def test_rvs_mixed_masses(self, default_sana12):
+        """Test random sampling with array of masses."""
+        rng = np.random.default_rng(42)
+        m1 = np.array([10.0, 15.0, 20.0, 25.0, 30.0])
+
+        samples = default_sana12.rvs(size=5, m1=m1, rng=rng)
+
+        assert len(samples) == 5
+        assert np.all(samples >= default_sana12.p_min)
+        assert np.all(samples <= default_sana12.p_max)
+
+    def test_rvs_without_rng(self, default_sana12):
+        """Test random sampling without providing an RNG."""
+        m1 = 20.0
+
+        samples = default_sana12.rvs(size=100, m1=m1)
+
+        assert len(samples) == 100
+        assert np.all(samples >= default_sana12.p_min)
+        assert np.all(samples <= default_sana12.p_max)
+
 
 class TestPowerLawPeriod:
     """Test class for PowerLawPeriod distribution."""
@@ -510,6 +589,24 @@ def test_normalization_consistency(self, default_power_law):
 
         np.testing.assert_allclose(pdf_values, expected_pdf)
 
+    def test_rvs(self, custom_power_law):
+        """Test random sampling."""
+        rng = np.random.default_rng(42)
+
+        samples = custom_power_law.rvs(size=1000, rng=rng)
+
+        assert len(samples) == 1000
+        assert np.all(samples >= custom_power_law.p_min)
+        assert np.all(samples <= custom_power_law.p_max)
+
+    def test_rvs_without_rng(self, custom_power_law):
+        """Test random sampling without providing an RNG."""
+        samples = custom_power_law.rvs(size=100)
+
+        assert len(samples) == 100
+        assert np.all(samples >= custom_power_law.p_min)
+        assert np.all(samples <= custom_power_law.p_max)
+
 
 class TestDistributionComparisons:
     """Test class for comparing distributions and edge cases."""
@@ -616,5 +713,471 @@ def test_initialization_invalid_parameters(self):
         with pytest.raises(ValueError, match="max must be greater than min"):
             LogUniform(min=1000.0, max=100.0)
 
+        with pytest.raises(ValueError, match="max must be greater than min"):
+            LogUniform(min=100.0, max=100.0)
+
+    def test_repr(self):
+        """Test string representation."""
+        log_uniform = LogUniform(min=10.0, max=1000.0)
+        rep_str = log_uniform.__repr__()
+        assert "LogUniform(" in rep_str
+        assert "min=10.0" in rep_str
+        assert "max=1000.0" in rep_str
+
+    def test_repr_html(self):
+        """Test HTML representation."""
+        log_uniform = LogUniform(min=10.0, max=1000.0)
+        html_str = log_uniform._repr_html_()
+        assert "<h3>Log-Uniform Distribution</h3>" in html_str
+        assert "min = 10.0" in html_str
+        assert "max = 1000.0" in html_str
+
+    def test_pdf_within_range(self):
+        """Test PDF within the valid range."""
+        log_uniform = LogUniform(min=10.0, max=1000.0)
+        x_values = np.array([10.0, 50.0, 100.0, 500.0, 1000.0])
+        pdf_values = log_uniform.pdf(x_values)
+        expected = log_uniform.norm / x_values
+        np.testing.assert_allclose(pdf_values, expected)
+
+    def test_pdf_outside_range(self):
+        """Test PDF outside the valid range."""
+        log_uniform = LogUniform(min=10.0, max=1000.0)
+        # Below range
+        x_below = np.array([1.0, 5.0])
+        pdf_below = log_uniform.pdf(x_below)
+        np.testing.assert_array_equal(pdf_below, np.zeros_like(x_below))
+
+        # Above range
+        x_above = np.array([2000.0, 5000.0])
+        pdf_above = log_uniform.pdf(x_above)
+        np.testing.assert_array_equal(pdf_above, np.zeros_like(x_above))
+
+    def test_rvs(self):
+        """Test random sampling."""
+        log_uniform = LogUniform(min=10.0, max=1000.0)
+        rng = np.random.default_rng(42)
+        samples = log_uniform.rvs(size=1000, rng=rng)
+
+        assert len(samples) == 1000
+        assert np.all(samples >= log_uniform.min)
+        assert np.all(samples <= log_uniform.max)
+
+    def test_rvs_without_rng(self):
+        """Test random sampling without providing an RNG."""
+        log_uniform = LogUniform(min=10.0, max=1000.0)
+        samples = log_uniform.rvs(size=100)
+
+        assert len(samples) == 100
+        assert np.all(samples >= log_uniform.min)
+        assert np.all(samples <= log_uniform.max)
+
+
+class TestThermalEccentricity:
+    """Test class for ThermalEccentricity distribution."""
+
+    @pytest.fixture
+    def default_thermal(self):
+        """Fixture for default ThermalEccentricity instance."""
+        return ThermalEccentricity()
+
+    @pytest.fixture
+    def custom_thermal(self):
+        """Fixture for custom ThermalEccentricity instance."""
+        return ThermalEccentricity(e_min=0.1, e_max=0.9)
+
+    def test_initialization_default(self, default_thermal):
+        """Test default initialization."""
+        assert default_thermal.e_min == 0.0
+        assert default_thermal.e_max == 1.0
+        assert hasattr(default_thermal, 'norm')
+        assert default_thermal.norm > 0
+
+    def test_initialization_custom(self, custom_thermal):
+        """Test custom initialization."""
+        assert custom_thermal.e_min == 0.1
+        assert custom_thermal.e_max == 0.9
+        assert hasattr(custom_thermal, 'norm')
+        assert custom_thermal.norm > 0
+
+    def test_initialization_invalid_parameters(self):
+        """Test that initialization raises ValueError for invalid parameters."""
+        # Test e_min not in [0, 1)
+        with pytest.raises(ValueError, match="e_min must be in \\[0, 1\\)"):
+            ThermalEccentricity(e_min=-0.1, e_max=0.5)
+
+        with pytest.raises(ValueError, match="e_min must be in \\[0, 1\\)"):
+            ThermalEccentricity(e_min=1.5, e_max=2.0)
+
+        # Test e_max not in (0, 1]
+        with pytest.raises(ValueError, match="e_max must be in \\(0, 1\\]"):
+            ThermalEccentricity(e_min=0.1, e_max=0.0)
+
+        with pytest.raises(ValueError, match="e_max must be in \\(0, 1\\]"):
+            ThermalEccentricity(e_min=0.1, e_max=-0.1)
+
+        with pytest.raises(ValueError, match="e_max must be in \\(0, 1\\]"):
+            ThermalEccentricity(e_min=0.1, e_max=1.5)
+
+        # Test e_min >= e_max
+        with pytest.raises(ValueError, match="e_min must be less than e_max"):
+            ThermalEccentricity(e_min=0.8, e_max=0.5)
+
+        with pytest.raises(ValueError, match="e_min must be less than e_max"):
+            ThermalEccentricity(e_min=0.5, e_max=0.5)
+
+    def test_repr(self, custom_thermal):
+        """Test string representation."""
+        rep_str = custom_thermal.__repr__()
+        assert "ThermalEccentricity(" in rep_str
+        assert "e_min=0.1" in rep_str
+        assert "e_max=0.9" in rep_str
+
+    def test_repr_html(self, custom_thermal):
+        """Test HTML representation."""
+        html_str = custom_thermal._repr_html_()
+        assert "<h3>Thermal Eccentricity Distribution</h3>" in html_str
+        assert "e_min = 0.1" in html_str
+        assert "e_max = 0.9" in html_str
+
+    def test_thermal_eccentricity_method(self, default_thermal):
+        """Test the thermal_eccentricity method."""
+        e_values = np.array([0.0, 0.25, 0.5, 0.75, 1.0])
+        result = default_thermal.thermal_eccentricity(e_values)
+        expected = 2.0 * e_values
+        np.testing.assert_allclose(result, expected)
+
+    def test_pdf_within_range(self, custom_thermal):
+        """Test PDF within the valid range."""
+        e_values = np.linspace(custom_thermal.e_min, custom_thermal.e_max, 10)
+        pdf_values = custom_thermal.pdf(e_values)
+
+        # All should be positive within range
+        assert np.all(pdf_values > 0)
+
+        # Check normalization
+        expected = custom_thermal.thermal_eccentricity(e_values) * custom_thermal.norm
+        np.testing.assert_allclose(pdf_values, expected)
+
+    def test_pdf_outside_range(self, custom_thermal):
+        """Test PDF outside the valid range."""
+        # Below range
+        e_below = np.array([0.0, 0.05])
+        pdf_below = custom_thermal.pdf(e_below)
+        np.testing.assert_array_equal(pdf_below, np.zeros_like(e_below))
+
+        # Above range
+        e_above = np.array([0.95, 1.0])
+        pdf_above = custom_thermal.pdf(e_above)
+        np.testing.assert_array_equal(pdf_above, np.zeros_like(e_above))
+
+    def test_pdf_scalar_input(self, default_thermal):
+        """Test PDF with scalar input."""
+        e = 0.5
+        pdf_value = default_thermal.pdf(e)
+        expected = default_thermal.thermal_eccentricity(e) * default_thermal.norm
+        np.testing.assert_allclose(pdf_value, expected)
+
+    def test_rvs(self, custom_thermal):
+        """Test random sampling."""
+        rng = np.random.default_rng(42)
+        samples = custom_thermal.rvs(size=1000, rng=rng)
+
+        assert len(samples) == 1000
+        assert np.all(samples >= custom_thermal.e_min)
+        assert np.all(samples <= custom_thermal.e_max)
+
+    def test_rvs_without_rng(self, custom_thermal):
+        """Test random sampling without providing an RNG."""
+        samples = custom_thermal.rvs(size=100)
+
+        assert len(samples) == 100
+        assert np.all(samples >= custom_thermal.e_min)
+        assert np.all(samples <= custom_thermal.e_max)
+
+
+class TestUniformEccentricity:
+    """Test class for UniformEccentricity distribution."""
+
+    @pytest.fixture
+    def default_uniform(self):
+        """Fixture for default UniformEccentricity instance."""
+        return UniformEccentricity()
+
+    @pytest.fixture
+    def custom_uniform(self):
+        """Fixture for custom UniformEccentricity instance."""
+        return UniformEccentricity(e_min=0.1, e_max=0.9)
+
+    def test_initialization_default(self, default_uniform):
+        """Test default initialization."""
+        assert default_uniform.e_min == 0.0
+        assert default_uniform.e_max == 1.0
+        assert hasattr(default_uniform, 'norm')
+        assert default_uniform.norm == 1.0
+
+    def test_initialization_custom(self, custom_uniform):
+        """Test custom initialization."""
+        assert custom_uniform.e_min == 0.1
+        assert custom_uniform.e_max == 0.9
+        assert hasattr(custom_uniform, 'norm')
+        expected_norm = 1.0 / (0.9 - 0.1)
+        np.testing.assert_allclose(custom_uniform.norm, expected_norm)
+
+    def test_initialization_invalid_parameters(self):
+        """Test that initialization raises ValueError for invalid parameters."""
+        # Test e_min not in [0, 1)
+        with pytest.raises(ValueError, match="e_min must be in \\[0, 1\\)"):
+            UniformEccentricity(e_min=-0.1, e_max=0.5)
+
+        with pytest.raises(ValueError, match="e_min must be in \\[0, 1\\)"):
+            UniformEccentricity(e_min=1.5, e_max=2.0)
+
+        # Test e_max not in (0, 1]
+        with pytest.raises(ValueError, match="e_max must be in \\(0, 1\\]"):
+            UniformEccentricity(e_min=0.1, e_max=0.0)
+
+        with pytest.raises(ValueError, match="e_max must be in \\(0, 1\\]"):
+            UniformEccentricity(e_min=0.1, e_max=-0.1)
+
+        with pytest.raises(ValueError, match="e_max must be in \\(0, 1\\]"):
+            UniformEccentricity(e_min=0.1, e_max=1.5)
+
+        # Test e_min >= e_max
+        with pytest.raises(ValueError, match="e_min must be less than e_max"):
+            UniformEccentricity(e_min=0.8, e_max=0.5)
+
+        with pytest.raises(ValueError, match="e_min must be less than e_max"):
+            UniformEccentricity(e_min=0.5, e_max=0.5)
+
+    def test_repr(self, custom_uniform):
+        """Test string representation."""
+        rep_str = custom_uniform.__repr__()
+        assert "UniformEccentricity(" in rep_str
+        assert "e_min=0.1" in rep_str
+        assert "e_max=0.9" in rep_str
+
+    def test_repr_html(self, custom_uniform):
+        """Test HTML representation."""
+        html_str = custom_uniform._repr_html_()
+        assert "<h3>Uniform Eccentricity Distribution</h3>" in html_str
+        assert "e_min = 0.1" in html_str
+        assert "e_max = 0.9" in html_str
+
+    def test_pdf_within_range(self, custom_uniform):
+        """Test PDF within the valid range."""
+        e_values = np.linspace(custom_uniform.e_min, custom_uniform.e_max, 10)
+        pdf_values = custom_uniform.pdf(e_values)
+
+        # All should equal the normalization constant
+        expected = custom_uniform.norm * np.ones_like(e_values)
+        np.testing.assert_allclose(pdf_values, expected)
+
+    def test_pdf_outside_range(self, custom_uniform):
+        """Test PDF outside the valid range."""
+        # Below range
+        e_below = np.array([0.0, 0.05])
+        pdf_below = custom_uniform.pdf(e_below)
+        np.testing.assert_array_equal(pdf_below, np.zeros_like(e_below))
+
+        # Above range
+        e_above = np.array([0.95, 1.0])
+        pdf_above = custom_uniform.pdf(e_above)
+        np.testing.assert_array_equal(pdf_above, np.zeros_like(e_above))
+
+    def test_pdf_scalar_input(self, default_uniform):
+        """Test PDF with scalar input."""
+        e = 0.5
+        pdf_value = default_uniform.pdf(e)
+        np.testing.assert_allclose(pdf_value, default_uniform.norm)
+
+    def test_rvs(self, custom_uniform):
+        """Test random sampling."""
+        rng = np.random.default_rng(42)
+        samples = custom_uniform.rvs(size=1000, rng=rng)
+
+        assert len(samples) == 1000
+        assert np.all(samples >= custom_uniform.e_min)
+        assert np.all(samples <= custom_uniform.e_max)
+
+    def test_rvs_without_rng(self, custom_uniform):
+        """Test random sampling without providing an RNG."""
+        samples = custom_uniform.rvs(size=100)
+
+        assert len(samples) == 100
+        assert np.all(samples >= custom_uniform.e_min)
+        assert np.all(samples <= custom_uniform.e_max)
+
+
+class TestZeroEccentricity:
+    """Test class for ZeroEccentricity distribution."""
+
+    @pytest.fixture
+    def zero_ecc(self):
+        """Fixture for ZeroEccentricity instance."""
+        return ZeroEccentricity()
+
+    def test_initialization(self, zero_ecc):
+        """Test initialization."""
+        # Should have no parameters
+        assert isinstance(zero_ecc, ZeroEccentricity)
+
+    def test_repr(self, zero_ecc):
+        """Test string representation."""
+        rep_str = zero_ecc.__repr__()
+        assert "ZeroEccentricity()" in rep_str
+
+    def test_repr_html(self, zero_ecc):
+        """Test HTML representation."""
+        html_str = zero_ecc._repr_html_()
+        assert "<h3>Zero Eccentricity Distribution</h3>" in html_str
+        assert "e = 0 (circular orbits)" in html_str
+
+    def test_pdf_at_zero(self, zero_ecc):
+        """Test PDF at e=0."""
+        pdf_value = zero_ecc.pdf(0.0)
+        assert pdf_value == 1.0
+
+    def test_pdf_away_from_zero(self, zero_ecc):
+        """Test PDF for non-zero eccentricities."""
+        e_values = np.array([0.1, 0.5, 0.9, 1.0])
+        pdf_values = zero_ecc.pdf(e_values)
+        np.testing.assert_array_equal(pdf_values, np.zeros_like(e_values))
+
+    def test_pdf_mixed(self, zero_ecc):
+        """Test PDF with mixture of zero and non-zero values."""
+        e_values = np.array([0.0, 0.1, 0.0, 0.5])
+        pdf_values = zero_ecc.pdf(e_values)
+        expected = np.array([1.0, 0.0, 1.0, 0.0])
+        np.testing.assert_array_equal(pdf_values, expected)
+
+    def test_rvs(self, zero_ecc):
+        """Test random sampling."""
+        rng = np.random.default_rng(42)
+        samples = zero_ecc.rvs(size=1000, rng=rng)
+
+        # All samples should be zero
+        assert len(samples) == 1000
+        np.testing.assert_array_equal(samples, np.zeros(1000))
+
+    def test_rvs_without_rng(self, zero_ecc):
+        """Test random sampling without providing an RNG."""
+        samples = zero_ecc.rvs(size=100)
+
+        # All samples should be zero
+        assert len(samples) == 100
+        np.testing.assert_array_equal(samples, np.zeros(100))
+
+
+class TestLogNormalSeparation:
+    """Test class for LogNormalSeparation distribution."""
+
+    @pytest.fixture
+    def default_lognormal(self):
+        """Fixture for default LogNormalSeparation instance."""
+        return LogNormalSeparation()
+
+    @pytest.fixture
+    def custom_lognormal(self):
+        """Fixture for custom LogNormalSeparation instance."""
+        return LogNormalSeparation(mean=1.0, sigma=0.5, min=10.0, max=1e4)
+
+    def test_initialization_default(self, default_lognormal):
+        """Test default initialization."""
+        assert default_lognormal.mean == 0.85
+        assert default_lognormal.sigma == 0.37
+        assert default_lognormal.min == 5.0
+        assert default_lognormal.max == 1e5
+
+    def test_initialization_custom(self, custom_lognormal):
+        """Test custom initialization."""
+        assert custom_lognormal.mean == 1.0
+        assert custom_lognormal.sigma == 0.5
+        assert custom_lognormal.min == 10.0
+        assert custom_lognormal.max == 1e4
+
+    def test_initialization_invalid_parameters(self):
+        """Test that initialization raises ValueError for invalid parameters."""
+        # Test min <= 0
+        with pytest.raises(ValueError, match="min must be positive"):
+            LogNormalSeparation(mean=1.0, sigma=0.5, min=0.0, max=1000.0)
+
+        with pytest.raises(ValueError, match="min must be positive"):
+            LogNormalSeparation(mean=1.0, sigma=0.5, min=-1.0, max=1000.0)
+
+        # Test max <= min
+        with pytest.raises(ValueError, match="max must be greater than min"):
+            LogNormalSeparation(mean=1.0, sigma=0.5, min=1000.0, max=100.0)
+
+        with pytest.raises(ValueError, match="max must be greater than min"):
+            LogNormalSeparation(mean=1.0, sigma=0.5, min=100.0, max=100.0)
+
+        # Test sigma <= 0
+        with pytest.raises(ValueError, match="sigma must be positive"):
+            LogNormalSeparation(mean=1.0, sigma=0.0, min=10.0, max=1000.0)
+
+        with pytest.raises(ValueError, match="sigma must be positive"):
+            LogNormalSeparation(mean=1.0, sigma=-0.5, min=10.0, max=1000.0)
+
+    def test_repr(self, custom_lognormal):
+        """Test string representation."""
+        rep_str = custom_lognormal.__repr__()
+        assert "LogNormalSeparation(" in rep_str
+        assert "mean=1.0" in rep_str
+        assert "sigma=0.5" in rep_str
+        assert "min=10.0" in rep_str
+        assert "max=10000.0" in rep_str
+
+    def test_repr_html(self, custom_lognormal):
+        """Test HTML representation."""
+        html_str = custom_lognormal._repr_html_()
+        assert "<h3>Log-Normal Separation Distribution</h3>" in html_str
+        assert "mean (log10) = 1.0" in html_str
+        assert "sigma (log10) = 0.5" in html_str
+
+    def test_pdf_within_range(self, custom_lognormal):
+        """Test PDF within the valid range."""
+        a_values = np.array([10.0, 50.0, 100.0, 500.0, 1000.0])
+        pdf_values = custom_lognormal.pdf(a_values)
+
+        # All should be positive within range
+        assert np.all(pdf_values > 0)
+
+    def test_pdf_outside_range(self, custom_lognormal):
+        """Test PDF outside the valid range."""
+        # Below range
+        a_below = np.array([1.0, 5.0])
+        pdf_below = custom_lognormal.pdf(a_below)
+        np.testing.assert_array_equal(pdf_below, np.zeros_like(a_below))
+
+        # Above range
+        a_above = np.array([2e4, 5e4])
+        pdf_above = custom_lognormal.pdf(a_above)
+        np.testing.assert_array_equal(pdf_above, np.zeros_like(a_above))
+
+    def test_pdf_zero_and_negative(self, custom_lognormal):
+        """Test PDF for zero and negative values."""
+        a_invalid = np.array([0.0, -10.0])
+        pdf_invalid = custom_lognormal.pdf(a_invalid)
+        np.testing.assert_array_equal(pdf_invalid, np.zeros_like(a_invalid))
+
+    def test_rvs(self, custom_lognormal):
+        """Test random sampling."""
+        rng = np.random.default_rng(42)
+        samples = custom_lognormal.rvs(size=1000, rng=rng)
+
+        assert len(samples) == 1000
+        assert np.all(samples >= custom_lognormal.min)
+        assert np.all(samples <= custom_lognormal.max)
+
+    def test_rvs_without_rng(self, custom_lognormal):
+        """Test random sampling without providing an RNG."""
+        samples = custom_lognormal.rvs(size=100)
+
+        assert len(samples) == 100
+        assert np.all(samples >= custom_lognormal.min)
+        assert np.all(samples <= custom_lognormal.max)
+
+
         with pytest.raises(ValueError, match="max must be greater than min"):
             LogUniform(min=100.0, max=100.0)
diff --git a/posydon/unit_tests/popsyn/test_norm_pop.py b/posydon/unit_tests/popsyn/test_norm_pop.py
index 1eb8f3fcce..a868d0d1b2 100644
--- a/posydon/unit_tests/popsyn/test_norm_pop.py
+++ b/posydon/unit_tests/popsyn/test_norm_pop.py
@@ -284,6 +284,32 @@ def test_mean_mass_without_q_bounds(self):
             'orbital_period_max': 6000,
         }
 
+        mean_mass = norm_pop.get_mean_mass(params)
+        assert mean_mass > 0
+
+    def test_computed_q_min_greater_than_q_max_error(self):
+        # Test the validation error when computed q_min > q_max
+        # This happens when secondary_mass_min/primary_mass_min > secondary_mass_max/primary_mass_max
+        params = {
+            'primary_mass_scheme': 'NonExistentIMF',
+            'primary_mass_min': 10,
+            'primary_mass_max': 20,
+            'secondary_mass_min': 15,  # Too high: 15/10 = 1.5
+            'secondary_mass_max': 5,   # Too low: 5/20 = 0.25
+            'secondary_mass_scheme': 'flat_mass_ratio',
+            'binary_fraction_scheme': 'const',
+            'binary_fraction_const': 0.5,
+            'orbital_scheme': 'period',
+            'orbital_period_scheme': 'Sana+12_period_extended',
+            'orbital_period_min': 0.35,
+            'orbital_period_max': 6000,
+        }
+
+        with pytest.raises(ValueError) as excinfo:
+            norm_pop.get_mean_mass(params)
+
+        assert "q_min must be less than q_max" in str(excinfo.value)
+
         # This should not raise an error and should return a valid mean mass
         result = norm_pop.get_mean_mass(params)
         assert isinstance(result, (float, np.floating))

From 878c98c4338ebb9d6fc5c9488328929951b00341 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:55:26 +0100
Subject: [PATCH 10/18] add repr and change logUniform pdf to include jacobian

---
 posydon/popsyn/distributions.py | 34 +++++++++++++++++++++++++++------
 1 file changed, 28 insertions(+), 6 deletions(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index ff0fdfd53b..1ae8c60d61 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -587,6 +587,28 @@ def __init__(self, min=5.0, max=1e5):
 
         self.norm = self._calculate_normalization()
 
+    def __repr__(self):
+        """Return string representation of the distribution.
+
+        Returns
+        -------
+        str
+            String representation showing the distribution parameters.
+        """
+        return f"LogUniform(min={self.min}, max={self.max})"
+
+    def _repr_html_(self):
+        """Return HTML representation for Jupyter notebooks.
+
+        Returns
+        -------
+        str
+            HTML string for rich display in notebooks.
+        """
+        return (f"<h3>Log-Uniform Distribution</h3>"
+                f"<p>min = {self.min}</p>"
+                f"<p>max = {self.max}</p>")
+
     def _calculate_normalization(self):
         """
         Calculate the normalization constant for the log-uniform distribution.
@@ -616,7 +638,7 @@ def pdf(self, x):
         valid = (x > 0) & (x >= self.min) & (x <= self.max)
         pdf_values = np.zeros_like(x, dtype=float)
 
-        pdf_values[valid] = self.norm
+        pdf_values[valid] = self.norm / x[valid]
 
         return pdf_values
 
@@ -948,9 +970,9 @@ class LogNormalSeparation:
     Parameters
     ----------
     mean : float, optional
-        Mean of the log10 distribution (default: 3.5, corresponding to ~3162 Rsun).
+        Mean of the log10 distribution (default: 0.85, corresponding to ~7.08 Rsun).
     sigma : float, optional
-        Standard deviation of the log10 distribution (default: 2.3).
+        Standard deviation of the log10 distribution (default: 0.37).
     min : float, optional
         Minimum orbital separation in solar radii (default: 5.0).
     max : float, optional
@@ -972,15 +994,15 @@ class LogNormalSeparation:
     normal distribution in log10 space.
     """
 
-    def __init__(self, mean=3.5, sigma=2.3, min=5.0, max=1e5):
+    def __init__(self, mean=0.85, sigma=0.37, min=5.0, max=1e5):
         """Initialize the log-normal separation distribution.
 
         Parameters
         ----------
         mean : float, optional
-            Mean of the log10 distribution (default: 3.5).
+            Mean of the log10 distribution (default: 0.85).
         sigma : float, optional
-            Standard deviation of the log10 distribution (default: 2.3).
+            Standard deviation of the log10 distribution (default: 0.37).
         min : float, optional
             Minimum orbital separation in solar radii (default: 5.0).
         max : float, optional

From 51219e9f7ba1a9542996c03fbc535acb870ba026 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:55:36 +0100
Subject: [PATCH 11/18] add repr and change logUniform pdf to include jacobian

---
 posydon/popsyn/distributions.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index 1ae8c60d61..85a13b1851 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -638,7 +638,7 @@ def pdf(self, x):
         valid = (x > 0) & (x >= self.min) & (x <= self.max)
         pdf_values = np.zeros_like(x, dtype=float)
 
-        pdf_values[valid] = self.norm / x[valid]
+        pdf_values[valid] = self.norm / x[valid]  # PDF is constant in log space, so divide by x for linear space
 
         return pdf_values
 

From 1806d719c7bf70fa7365ac5a6aa7b874a39031ca Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 11:59:11 +0100
Subject: [PATCH 12/18] remove duplicate test

---
 posydon/unit_tests/popsyn/test_norm_pop.py | 27 ----------------------
 1 file changed, 27 deletions(-)

diff --git a/posydon/unit_tests/popsyn/test_norm_pop.py b/posydon/unit_tests/popsyn/test_norm_pop.py
index a868d0d1b2..db9f907850 100644
--- a/posydon/unit_tests/popsyn/test_norm_pop.py
+++ b/posydon/unit_tests/popsyn/test_norm_pop.py
@@ -287,33 +287,6 @@ def test_mean_mass_without_q_bounds(self):
         mean_mass = norm_pop.get_mean_mass(params)
         assert mean_mass > 0
 
-    def test_computed_q_min_greater_than_q_max_error(self):
-        # Test the validation error when computed q_min > q_max
-        # This happens when secondary_mass_min/primary_mass_min > secondary_mass_max/primary_mass_max
-        params = {
-            'primary_mass_scheme': 'NonExistentIMF',
-            'primary_mass_min': 10,
-            'primary_mass_max': 20,
-            'secondary_mass_min': 15,  # Too high: 15/10 = 1.5
-            'secondary_mass_max': 5,   # Too low: 5/20 = 0.25
-            'secondary_mass_scheme': 'flat_mass_ratio',
-            'binary_fraction_scheme': 'const',
-            'binary_fraction_const': 0.5,
-            'orbital_scheme': 'period',
-            'orbital_period_scheme': 'Sana+12_period_extended',
-            'orbital_period_min': 0.35,
-            'orbital_period_max': 6000,
-        }
-
-        with pytest.raises(ValueError) as excinfo:
-            norm_pop.get_mean_mass(params)
-
-        assert "q_min must be less than q_max" in str(excinfo.value)
-
-        # This should not raise an error and should return a valid mean mass
-        result = norm_pop.get_mean_mass(params)
-        assert isinstance(result, (float, np.floating))
-        assert result > 0
 
 class TestGetPdf:
     def test_single_star_pdf(self):

From b132392c4a9c179fe1f81de58ec3b33dd0b7397e Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 12:04:49 +0100
Subject: [PATCH 13/18] add pragma no cover

---
 posydon/popsyn/IMFs.py          | 3 ++-
 posydon/popsyn/distributions.py | 2 +-
 2 files changed, 3 insertions(+), 2 deletions(-)

diff --git a/posydon/popsyn/IMFs.py b/posydon/popsyn/IMFs.py
index f4467b168c..d09cf818a7 100644
--- a/posydon/popsyn/IMFs.py
+++ b/posydon/popsyn/IMFs.py
@@ -102,7 +102,8 @@ def imf(self, m): # pragma: no cover
         '''
         pass
 
-    def rvs(self, size=1, rng=None):
+    @abstractmethod
+    def rvs(self, size=1, rng=None): # pragma: no cover
         pass
 
 
diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index 85a13b1851..d8f1b66e5d 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -739,7 +739,7 @@ def _calculate_normalization(self):
         """
         # Integral of 2*e from e_min to e_max is e_max^2 - e_min^2
         integral = self.e_max**2 - self.e_min**2
-        if integral == 0:
+        if integral == 0: # pragma: no cover
             raise ValueError("Cannot normalize distribution: e_min == e_max")
         return 1.0 / integral
 

From 4fd1ba06012d2925209e4bc2d994328869899c3e Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Mon, 9 Feb 2026 12:07:06 +0100
Subject: [PATCH 14/18] add unite test from q_min>q_max from m2/m1 calculation

---
 posydon/unit_tests/popsyn/test_norm_pop.py | 23 ++++++++++++++++++++++
 1 file changed, 23 insertions(+)

diff --git a/posydon/unit_tests/popsyn/test_norm_pop.py b/posydon/unit_tests/popsyn/test_norm_pop.py
index db9f907850..e37c8e2b51 100644
--- a/posydon/unit_tests/popsyn/test_norm_pop.py
+++ b/posydon/unit_tests/popsyn/test_norm_pop.py
@@ -267,6 +267,29 @@ def test_q_min_greater_than_q_max_error(self):
 
         assert "q_min must be less than q_max" in str(excinfo.value)
 
+    def test_q_min_greater_than_q_max_computed_error(self):
+        # Test the validation error when computed q_min > q_max
+        # This happens when secondary_mass_min/primary_mass_min > secondary_mass_max/primary_mass_max
+        params = {
+            'primary_mass_scheme': 'NonExistentIMF',
+            'primary_mass_min': 5,
+            'primary_mass_max': 10,
+            'secondary_mass_min': 4,  # 4/5 = 0.8
+            'secondary_mass_max': 6,  # 6/10 = 0.6, so q_min (0.8) > q_max (0.6)
+            'secondary_mass_scheme': 'flat_mass_ratio',
+            'binary_fraction_scheme': 'const',
+            'binary_fraction_const': 0.5,
+            'orbital_scheme': 'period',
+            'orbital_period_scheme': 'Sana+12_period_extended',
+            'orbital_period_min': 0.35,
+            'orbital_period_max': 6000,
+        }
+
+        with pytest.raises(ValueError) as excinfo:
+            norm_pop.get_mean_mass(params)
+
+        assert "q_min must be less than q_max" in str(excinfo.value)
+
     def test_mean_mass_without_q_bounds(self):
         # Test the branch where q_min and q_max are computed from secondary masses
         params = {

From 85761a056296d47cdc3c9c02ef5d91345a63ec2d Mon Sep 17 00:00:00 2001
From: Max <14039563+maxbriel@users.noreply.github.com>
Date: Thu, 12 Mar 2026 14:52:13 +0100
Subject: [PATCH 15/18] Update posydon/popsyn/distributions.py

Co-authored-by: Seth Gossage <sethg45@gmail.com>
---
 posydon/popsyn/distributions.py | 3 ++-
 1 file changed, 2 insertions(+), 1 deletion(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index d8f1b66e5d..d7408b520a 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -363,7 +363,8 @@ def rvs(self, size=1, m1=None, rng=None):
         if m1.size == 1:
             m1 = np.full(size, m1[0])
         elif m1.size != size:
-            raise ValueError(f"m1 must be a single value or have size={size}")
+            raise ValueError(f"m1 must be a single value or have size={size}"
+                                      f"\n m1 = {m1}\n m1.size = {m1.size}")
 
         # Import here to avoid circular dependency
         from posydon.utils.common_functions import rejection_sampler

From 80c4410a427eaaab8635da5836e739684ee9c8bb Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Thu, 12 Mar 2026 15:03:32 +0100
Subject: [PATCH 16/18] remove unused import

---
 posydon/popsyn/IMFs.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/posydon/popsyn/IMFs.py b/posydon/popsyn/IMFs.py
index d09cf818a7..bfaace0213 100644
--- a/posydon/popsyn/IMFs.py
+++ b/posydon/popsyn/IMFs.py
@@ -9,7 +9,7 @@
 import numpy as np
 from scipy.integrate import quad
 
-from posydon.utils.common_functions import inverse_sampler, rejection_sampler
+from posydon.utils.common_functions import inverse_sampler
 from posydon.utils.posydonwarning import Pwarn
 
 

From 1f7a56b019689d81b59733dc0dba6a357b00bdfd Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Thu, 12 Mar 2026 16:06:22 +0100
Subject: [PATCH 17/18] make lower q boundary exclusive

---
 posydon/popsyn/norm_pop.py | 2 +-
 1 file changed, 1 insertion(+), 1 deletion(-)

diff --git a/posydon/popsyn/norm_pop.py b/posydon/popsyn/norm_pop.py
index da195213a6..fbaee17c32 100644
--- a/posydon/popsyn/norm_pop.py
+++ b/posydon/popsyn/norm_pop.py
@@ -101,7 +101,7 @@ def get_pdf_for_m1(m1):
                 axis=0)
 
             # Use FlatMassRatio distribution class
-            q_dist = lambda q: np.where((q >= minimum) & (q <= maximum),
+            q_dist = lambda q: np.where((q > minimum) & (q <= maximum),
                                         1/(maximum - minimum),
                                         0)
             return q_dist

From e1e2213e798288928d0f08a3d936ede9b8a4b8a4 Mon Sep 17 00:00:00 2001
From: Max Briel <max.briel@gmail.com>
Date: Thu, 19 Mar 2026 15:40:23 +0100
Subject: [PATCH 18/18] change description

---
 posydon/popsyn/distributions.py | 4 ++--
 1 file changed, 2 insertions(+), 2 deletions(-)

diff --git a/posydon/popsyn/distributions.py b/posydon/popsyn/distributions.py
index d7408b520a..8ac96ab764 100644
--- a/posydon/popsyn/distributions.py
+++ b/posydon/popsyn/distributions.py
@@ -13,12 +13,12 @@ class FlatMassRatio:
 
     A uniform distribution for mass ratios q = m2/m1 within specified bounds.
     This distribution assigns equal probability to all mass ratios within the
-    given range [q_min, q_max].
+    given range (q_min, q_max], exclusive bottom, inclusive top.
 
     Parameters
     ----------
     q_min : float, optional
-        Minimum mass ratio (default: 0.05). Must be in (0, 1].
+        Minimum mass ratio (default: 0.05). Must be in [0, 1).
     q_max : float, optional
         Maximum mass ratio (default: 1.0). Must be in (0, 1].