debrispy package

Submodules

debrispy.asd module

class debrispy.asd.ASD(kernel, sigma_a)

Bases: object

Computes the azimuthally averaged surface density profile, ASD, from a given semi-major axis surface density Sigma(a) and kernel Phi(r, a).

This class supports all kernel types.

For most use simple cases, the compute_quadvec method is the most convenient, and does not require any additional tuning of parameters, or knowledge of the integration limits.

For more complex cases, the compute_gl method is more flexible, this uses adaptive Gauss-Legendre quadrature, the user can specify the number of integration points, and the relative and absolute tolerances. Adaptive limits are supported and recommended for cases with very sharp features in the eccentricity profile.

Adaptive limits are not supported for calculations involving eccentricity distributions (only for unique eccentricity profiles).

Parameters:

• kernel (Kernel) – A kernel object used to compute Phi(r, a).

• sigma_a (SigmaA) – A surface density profile object defining Sigma(a).

compute_gl(r_vals=None, n_points: int = 64, tol_rel: float = 1e-08, tol_abs: float = 1e-08, max_level: int = 15, n_jobs: int = -1, show_progress: bool = True, pad: float = 0.05, rf: float = 5.0, adaptive_limits: bool = False, batch_size: int = 10, verbose: bool = True)

Compute the azimuthally averaged surface density profile, ASD, using adaptive Gauss-Legendre quadrature.

Integration is performed using recursive adaptive Gauss-Legendre quadrature with optional adaptive integration bounds for efficiency (recommended for unique eccentricity profiles).

Parameters:

• r_vals (array-like or None) – Array of radius values at which to compute the ASD. If None, defaults to 500 evenly spaced points between kernel.r_min and kernel.r_max.

• n_points (int, optional) – Number of Gauss-Legendre points per subinterval (default: 64).

• tol_rel (float, optional) – Relative tolerance for adaptive integration (default: 1e-8).

• tol_abs (float, optional) – Absolute tolerance for adaptive integration (default: 1e-8).

• max_level (int, optional) – Maximum recursion depth for adaptive integration (default: 15).

• n_jobs (int, optional) – Number of parallel jobs to use (-1 for all available CPUs).

• show_progress (bool, optional) – Whether to show a progress bar using tqdm (default: True).

• pad (float, optional) – Padding added to adaptive integration bounds as a fraction of a (default: 0.05).

• rf (float, optional) – Rescue factor: multiplier for n_points in rescue passes for suspicious results (default: 5.0).

• adaptive_limits (bool, optional) – Whether to use adaptive integration bounds based on the eccentricity profile (default: False). Only supported for UniqueEccentricity kernels.

• batch_size (int, optional) – Batch size for parallel integration jobs (default: 10).

• verbose (bool, optional) – Whether to print verbose output (default: True).

Raises:

• ValueError – If adaptive limits are enabled but the kernel does not use a unique eccentricity profile.

• Prints –

• ------ –

• - Progress bar (optional) and integration diagnostics. –

• - Warnings about any suspicious or unresolved zero results. –

Returns:

Results are stored in self._r_vals and self._sigma_r_vals.

Return type:

None

compute_quadvec(r_vals: float | ndarray[tuple[int, ...], dtype[float64]], tol_rel: float = 1e-08, tol_abs: float = 1e-08) → None

Compute and cache the azimuthally averaged surface density using vectorised adaptive quadrature via scipy.integrate.quad_vec.

This method uses non-adaptive integration limits.

Parameters:

• r_vals (float or array-like) – Radius or array of radii at which to compute ASD. The result is cached internally and accessible via self._sigma_r_vals.

• tol_rel (float, optional) – Relative tolerance for the integration.

• tol_abs (float, optional) – Absolute tolerance for the integration.

convolve(width: float, M: int = 2048) → ndarray

Apply a Gaussian convolution to the computed ASD profile.

Parameters:

• width (float) – The Gaussian kernel width (in the same units as r) to smooth ASD.

• M (int, optional) – Number of uniform samples used if r-grid is non-uniform. Default is 2048.

Returns:

The convolved ASD values at the original r grid.

Return type:

np.ndarray

get_values()

Return the cached ASD values.

Returns:

The cached azimuthally averaged surface density values.

Return type:

ndarray

Raises:

RuntimeError if compute_sigma_r() has not been called yet. –

integrand(a: float | ndarray[tuple[int, ...], dtype[float64]], r_vals: float | ndarray[tuple[int, ...], dtype[float64]]) → ndarray[tuple[int, ...], dtype[float64]]

Compute the 2D integrand over a grid of a and r values.

Parameters:

• a (float or array) – Semi-major axis value(s) at which to evaluate the integrand.

• r_vals (float or array) – Radial location(s) at which the azimuthally averaged profile is to be computed.

Returns:

The evaluated integrand values over all combinations of a and r_vals.

Return type:

(A, R) array

plot(a_vals: ndarray, plot_a: bool = True, x_lim: Tuple[float, float] | None = None, y_lim: Tuple[float, float] | None = None, save: bool = False, filename: str | None = None, title: str | None = None, grid_hist: bool = False) → None

Plot Sigma(a) and ASD on the same figure, optionally including a histogram of the r-grid point distribution.

Parameters:

• a_vals (array-like) – Array of semi-major axis values at which to evaluate Σ(a).

• plot_a (bool, default=True) – Whether to plot Sigma(a) alongside ASD.

• x_lim (tuple, optional) – Limits for the x-axis.

• y_lim (tuple, optional) – Limits for the y-axis.

• save (bool, default=False) – If True, saves the figure to file instead of displaying it.

• filename (str, optional) – Name of the output file (without extension) if saving the figure.

• title (str, optional) – Title for the plot.

• grid_hist (bool, default=False) – If True, adds a histogram of r grid points above the main plot.

plot_convolution(width: float | None = None, M: int | None = None, overlay_sigma_r: bool = True, overlay_sigma_a: bool = False, x_lim: Tuple[float, float] | None = None, y_lim: Tuple[float, float] | None = None, save: bool = False, filename: str | None = None, title: str | None = None) → None

Plot the original and convolved ASD profiles.

Parameters:

• width (float, optional) – The Gaussian kernel width (in the same units as r) to smooth ASD. Required if convolve() has not been called yet.

• M (int, optional) – Number of uniform samples used if r-grid is non-uniform. Required if convolve() has not been called yet.

• overlay_sigma_r (bool, optional) – If True, overlay the original ASD profile on the convolved profile.

• overlay_sigma_a (bool, optional) – If True, overlay the Sigma_a profile on the convolved profile.

Returns:

Displays the plot.

Return type:

None

plot_integrand(a_vals: ndarray[tuple[int, ...], dtype[float64]] | None = None, r_vals: ndarray[tuple[int, ...], dtype[float64]] | None = None, *, save: bool = False, filename: str | None = None, vmin: float | None = None, vmax: float | None = None, log: bool = True, cmap: str = 'viridis', colorbar: bool = True, r_lim: Tuple[float, float] | None = None, a_lim: Tuple[float, float] | None = None, **imshow_kwargs) → None

Plot the integrand as a 2D colourmap over (r, a). Useful for diagnosing sharp features or integration difficulties.

Parameters:

• a_vals (array-like, optional) – Semi-major axis values to evaluate the integrand over. If None, uses a default linspace.

• r_vals (array-like, optional) – Radius values to evaluate the integrand over. If None, uses a default linspace.

• save (bool, optional) – If True, the plot is saved to disk. Default is False.

• filename (str, optional) – File path to save the plot. Required if save is True.

• vmin (float, optional) – Minimum value for the colour scale. If None, uses the minimum non-zero value.

• vmax (float, optional) – Maximum value for the colour scale. If None, uses the maximum value.

• log (bool, optional) – If True, applies logarithmic colour scaling (ignoring non-positive values). Default is True.

• cmap (str, optional) – Colormap to use for the image. Default is “viridis”.

• colorbar (bool, optional) – If True, adds a colourbar to the plot. Default is True.

• r_lim (tuple of float, optional) – Limits for the x-axis (r-axis).

• a_lim (tuple of float, optional) – Limits for the y-axis (a-axis).

• **imshow_kwargs (dict) – Additional keyword arguments passed to imshow().

Raises:

ValueError – If save=True but filename is not provided.

refine(curvature_factor: float = 1.0, max_rounds: int = 3, subdiv: int = 2, n_jobs: int = -1, show_progress: bool = False, n_points: int = 64, tol_rel: float = 1e-08, tol_abs: float = 1e-08, max_level: int = 25, adaptive_limits: bool = False, rf: float = 10.0, pad: float = 0.05, batch_size: int = 10)

Refine the existing ASD grid by adding intermediate points based on curvature.

This method uses a curvature-driven refinement strategy: it identifies regions where the ASD changes rapidly and adds more sampling points in those regions. Optionally, it applies adaptive integration bounds based on the eccentricity profile and performs a rescue pass on suspicious zeros.

This uses the same adaptive Gauss-Legendre quadrature as the compute_sigma_r method. Currently, this is only available for unique eccentricity profiles.

Parameters:

• curvature_factor (float, optional) – Threshold multiplier for triggering refinement based on curvature.

• max_rounds (int, optional) – Maximum number of refinement rounds to perform. Each round may add new points.

• subdiv (int, optional) – Number of subdivisions to insert between points with high curvature. For example, subdiv=2 inserts one point between each flagged pair.

• n_jobs (int, optional) – Number of parallel jobs to use for integration (-1 = all available CPUs).

• show_progress (bool, optional) – Whether to show a progress bar during integration.

• n_points (int, optional) – Number of Gauss-Legendre nodes to use in initial integration.

• tol_rel (float, optional) – Relative tolerance for the adaptive integrator.

• tol_abs (float, optional) – Absolute tolerance for the adaptive integrator.

• max_level (int, optional) – Maximum recursion depth for the adaptive Gauss-Legendre integrator.

• adaptive_limits (bool, optional) – If True, use adaptive (localised) integration bounds based on the eccentricity profile. Requires that the kernel’s eccentricity profile is UniqueEccentricity.

• rf (float, optional) – Rescue factor: multiply n_points by this value during the rescue pass to resolve suspicious zeros.

• pad (float, optional) – Padding to apply to adaptive integration bounds (in units of semi-major axis).

• batch_size (int, optional) – Batch size for parallel job distribution.

Raises:

ValueError – If adaptive_limits=True but the kernel does not have a UniqueEccentricity profile.

Returns:

Updates the internal self._r_vals and self._sigma_r_vals arrays with the refined grid.

Return type:

None

debrispy.asd.adapt(a_lo: float, a_hi: float, level: int, x1: ndarray[tuple[int, ...], dtype[float64]], w1: ndarray[tuple[int, ...], dtype[float64]], x2: ndarray[tuple[int, ...], dtype[float64]], w2: ndarray[tuple[int, ...], dtype[float64]], tol_rel: float, tol_abs: float, max_level: int, kernel, sigma_a, r_j: float) → float

Perform Gauss-Legendre quadrature at two levels (lower and higher order) and adaptively refine the result until the desired tolerance is achieved or the maximum recursion depth is reached.

Parameters:

• a_lo (float) – Lower limit of integration.

• a_hi (float) – Upper limit of integration.

• level (int) – Current recursion depth.

• x1 ((N,) arrays) – Lower-order Gauss-Legendre nodes and weights.

• w1 ((N,) arrays) – Lower-order Gauss-Legendre nodes and weights.

• x2 ((N,) arrays) – Higher-order Gauss-Legendre nodes and weights.

• w2 ((N,) arrays) – Higher-order Gauss-Legendre nodes and weights.

• tol_rel (float) – Relative tolerance.

• tol_abs (float) – Absolute tolerance.

• max_level (int) – Maximum recursion depth.

• kernel (Kernel) – Kernel object.

• sigma_a (SigmaA) – Surface density profile.

• r_j (float) – Radius value at which the integral is being computed.

Returns:

Integral value of ASD.

Return type:

float

debrispy.asd.chunked(iterable: Iterable[T], batch_size: int) → Iterator[List[T]]

Return successive chunks of a given size from an iterable.

Parameters:

• iterable (Iterable[T]) – Input iterable to be split into chunks.

• batch_size (int) – Maximum size of each chunk.

Returns:

Batches of the original iterable as lists of size at most batch_size.

Return type:

Iterator[List[T]]

debrispy.asd.gl_quad(a_lo: float, a_hi: float, nodes: ndarray[tuple[int, ...], dtype[float64]], weights: ndarray[tuple[int, ...], dtype[float64]], kernel, sigma_a, r_j: float) → float

Perform Gauss-Legendre quadrature for ASD over [a_lo, a_hi].

Parameters:

• a_lo (float) – Lower integration bound.

• a_hi (float) – Upper integration bound.

• nodes ((N,) array) – Gauss-Legendre nodes in [-1, 1].

• weights ((N,) array) – Corresponding Gauss-Legendre weights.

• kernel (Kernel object) – Phi(r, a) evaluator.

• sigma_a (SigmaA object) – Sigma(a) surface density profile.

• r_j (float) – Evaluation point in radius.

Returns:

Integral value of ASD.

Return type:

float

debrispy.eccentricity module

class debrispy.eccentricity.EccentricityDistribution(a_min: float, a_max: float, distribution_func: Callable[[ndarray[tuple[int, ...], dtype[float64]], ndarray[tuple[int, ...], dtype[float64]]], ndarray[tuple[int, ...], dtype[float64]]], num_e_points: int = 1000, num_a_points: int = 1000, interpolation_method: str | None = None, auto_normalise: bool = False, grid_type: str | None = None, grid_spread: float | None = 1.0)

Bases: EccentricityProfile

Defines an eccentricity distribution ψ_e(e,a). Gridding and interpolation is required if auto_normalise is True, otherwise the distribution is determined directly via the distribution_func.

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• distribution_func (callable) – Function ψ_e(e,a) to sample; must accept array inputs for both e and a.

• num_e_points (int, optional) – Number of points in e (default 1000).

• num_a_points (int, optional) – Number of points in a (default 1000).

• interpolation_method (str, optional) – ‘nearest’, ‘linear’ or ‘cubic’.

• auto_normalise (bool, optional) – If True, normalise ψ_e along e for each a, this is done via gridding and interpolation.

• grid_type (str, optional) – ‘uniform’, ‘warped’ or ‘adaptive’. (If auto_normalise is True, this must be provided.)

• grid_spread (float, optional) – This is a parameter used in warped and adaptive gridding. The larger the grid_spread, the less concentrated the grid is around sharp features.

distribution(e: float | ndarray[tuple[int, ...], dtype[float64]], a: float | ndarray[tuple[int, ...], dtype[float64]]) → float | ndarray[tuple[int, ...], dtype[float64]]

Evaluate the distribution at the provided points.

If auto_normalise is True, the distribution is evaluated using the interpolator. If auto_normalise is False, the distribution is evaluated using analytic distribution function.

Parameters:

• e (float or array-like) – Eccentricity values.

• a (float or array-like) – Semi-major axis values.

Returns:

psi_e – Distribution values.

Return type:

float or array-like

eccentricity(a)

Return the eccentricity e(a) at given semi-major axis a (for unique eccentricity profiles). Should be overridden by subclasses if applicable.

get_sampled_distribution() → Tuple[ndarray[tuple[int, ...], dtype[float64]], ndarray[tuple[int, ...], dtype[float64]], ndarray[tuple[int, ...], dtype[float64]]]

Get the sampled distribution (e_grid, a_grid, psi_grid).

If auto_normalise is True, the sampled distribution is the same as the grid data. If auto_normalise is False, the sampled distribution is the distribution function evaluated on the grid.

Returns:

• e_grid (npt.NDArray[np.float64]) – Eccentricity grid.

• a_grid (npt.NDArray[np.float64]) – Semi-major axis grid.

• psi_grid (npt.NDArray[np.float64]) – Distribution values.

plot(save: bool = False, filename: str | None = None, log: bool = False, vmin: float | None = None, vmax: float | None = None, points: bool = False, point_size: float = 10, cmap: str = 'viridis') → None

Plot the eccentricity distribution.

Parameters:

• save (bool) – If True, save figure to filename.

• filename (str) – Path to save the figure.

• show_grid (bool) – Overlay the computational grid (only if auto_normalise=True).

• log (bool) – Use logarithmic color scale.

• vmin (float) – Color scale limits.

• vmax (float) – Color scale limits.

• points (bool) – Plot the distribution as points instead of a surface.

• point_size (float) – Size of the points if points=True.

• cmap (str) – Colormap to use.

plot_slice(*, fix_a: float | None = None, fix_e: float | None = None, num_points: int = 500, save: bool = False, filename: str | None = None, figsize: Tuple[int, int] = (8, 6), show: bool = True, ax: Axes | None = None, **plot_kwargs) → None

Plot a 1D slice of the 2D eccentricity distribution psi(e, a).

Parameters:

• fix_a (float, optional) – If provided, plots psi(e) at fixed semi-major axis a.

• fix_e (float, optional) – If provided, plots psi(a) at fixed eccentricity e.

• num_points (int, optional) – Number of points for plotting. Defaults to 500.

• save (bool, optional) – If True, saves the figure instead of displaying. Defaults to False.

• filename (str, optional) – Filename to save the figure. Required if save is True.

• figsize (tuple, optional) – Figure size in inches. Defaults to (8, 6).

• show (bool, optional) – Whether to show the plot. Defaults to True.

• ax (matplotlib.axes.Axes, optional) – Axes to plot on. If None, a new figure is created.

• **plot_kwargs (dict) – Additional keyword arguments for plt.plot.

Raises:

ValueError – If neither or both of fix_a and fix_e are provided.

class debrispy.eccentricity.EccentricityProfile

Bases: object

Abstract base class for eccentricity profiles.

distribution(e, a)

Return the eccentricity distribution ψ_e(e,a) (for distribution-based eccentricity profiles). Should be overridden by subclasses if applicable.

eccentricity(a)

Return the eccentricity e(a) at given semi-major axis a (for unique eccentricity profiles). Should be overridden by subclasses if applicable.

plot(*args, **kwargs)

Plot the eccentricity profile.

class debrispy.eccentricity.PowerLawEccentricity(a_min: float, a_max: float, zeta: float, lam: float | Callable[[float], float], num_e_points: int = 1000, num_a_points: int = 1000)

Bases: EccentricityDistribution

Eccentricity distribution psi(e,a) = (2*zeta + 1) * lambda(a)^(-(2*zeta + 1)) * e * (lambda(a)^2 - e^2)^(zeta - 1/2) defined for 0 ≤ e < lambda(a), with zeta > 0.

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• zeta (float) – Power-law shape parameter (must satisfy ζ > 0).

• lam (callable or float) – Function lambda(a). User can provide a callable or a constant lambda.

• num_e_points (int, optional) – Number of eccentricity grid points.

• num_a_points (int, optional) – Number of semi-major axis grid points.

class debrispy.eccentricity.RayleighEccentricity(a_min: float, a_max: float, sigma0: float | None = None, power: float | None = None, sigma_func: Callable[[float], float] | None = None, num_e_points: int = 1000, num_a_points: int = 1000)

Bases: EccentricityDistribution

Eccentricity distribution assuming a properly normalized Rayleigh distribution.

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• sigma0 (float, optional) – Amplitude of the power-law scale function sigma(a). Required if sigma_func is not provided.

• power (float, optional) – Power-law slope of sigma(a). Required if sigma_func is not provided.

• sigma_func (callable, optional) – Custom function sigma(a). If provided, sigma0 and power must be omitted.

• num_e_points (int) – Number of eccentricity grid points.

• num_a_points (int) – Number of semi-major axis grid points.

class debrispy.eccentricity.TopHatEccentricity(a_min: float, a_max: float, lam: float | Callable[[float], float], num_e_points: int = 1000, num_a_points: int = 1000)

Bases: EccentricityDistribution

Eccentricity distribution \(\psi(e,a) = (1 / \lambda(a)) * e / \sqrt(\lambda(a)^2 - e^2)\), which yields a step-function kernel \(\phi(kappa,a) = (\pi / 2\lambda(a)) · 1_{kappa <= \lambda(a)}\).

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• lam (callable or float) – The user can either provide a callable lambda(a) or a constant lambda0.

• num_e_points (int, optional) – Number of eccentricity grid points.

• num_a_points (int, optional) – Number of semi-major axis grid points.

class debrispy.eccentricity.TriangularEccentricity(a_min: float, a_max: float, lam: float | Callable[[float], float], num_e_points: int = 1000, num_a_points: int = 1000)

Bases: EccentricityDistribution

Eccentricity distribution: \(\psi(e,a) = (2e / \lambda(a)^2) * \ln[ (\lambda(a) + \sqrt(\lambda(a)^2 - e^2)) / e ]\), valid for \(0 \leq e < \lambda(a)\)

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• lam (callable or float) – The user can either provide a callable lambda(a) or a constant lambda0.

• num_e_points (int, optional) – Number of eccentricity grid points.

• num_a_points (int, optional) – Number of semi-major axis grid points.

class debrispy.eccentricity.TruncGaussEccentricity(a_min: float, a_max: float, sigma: float | Callable[[float], float], lam: float | Callable[[float], float], num_e_points: int = 1000, num_a_points: int = 1000)

Bases: EccentricityDistribution

Truncated Gaussian eccentricity distribution with a normalisation term.

psi(e, a) = sqrt(2/pi) * C(a) * [

exp(-e² / (2sigma_k(a)²)) / sigma_k(a) * erf( sqrt((lambda(a)² - e²) / (2sigma_k(a)²)) ) + sqrt(2/pi) * exp(-lambda(a)² / (2sigma_k(a)²)) / sqrt(lambda(a)² - e²)

]

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• sigma (float or callable) – Constant sigma or a function sigma(a). Must be > 0.

• lam (float or callable) – Constant lambda or a function lambda(a). Must be 0 < lambda ≤ 1.

• num_e_points (int) – Number of eccentricity grid points.

• num_a_points (int) – Number of semi-major axis grid points.

class debrispy.eccentricity.UniqueEccentricity(a_min: float, a_max: float, eccentricity_func: Callable[[float | ndarray[tuple[int, ...], dtype[float64]]], ndarray[tuple[int, ...], dtype[float64]]] | None = None, e0: float = 1.0, power: float | None = None)

Bases: EccentricityProfile

Fixed eccentricity profile where eccentricity is a unique function of semi-major axis: e = e(a).

By default implements a power-law profile e(a) = e0 * (a_min/a)^{power}, but accepts custom user-defined functions.

Parameters:

• a_min (float) – Minimum semi-major axis.

• a_max (float) – Maximum semi-major axis.

• eccentricity_func (callable, optional) – Custom function e(a). If None, uses a default power-law form.

• e0 (float, optional) – Normalization constant. Default is 1.0.

• power (float, optional) – Power-law exponent for built-in power-law e(a) (required if eccentricity_func is None).

Raises:

ValueError – If the eccentricity function produces values outside the range [0, 1) across the specified semi-major axis range.

derivative(a: float | ndarray[tuple[int, ...], dtype[float64]]) → float | ndarray[tuple[int, ...], dtype[float64]]

Calculate the derivative de/da of the eccentricity profile using analytical or finite differences.

Parameters:

a (float or array-like) – Semi-major axis value(s).

Returns:

Derivative value(s) at the specified semi-major axis/axes.

Return type:

float or ndarray

distribution(e: float | ndarray[tuple[int, ...], dtype[float64]], a: float | ndarray[tuple[int, ...], dtype[float64]]) → None

Distribution function ψ_e(e,a) is not defined for deterministic profiles.

Raises:

NotImplementedError – Always raised since fixed eccentricity profiles do not have a distribution.

eccentricity(a: float | ndarray[tuple[int, ...], dtype[float64]]) → float | ndarray[tuple[int, ...], dtype[float64]]

Return eccentricity value(s) at given semi-major axis/axes.

Parameters:

a (float or array-like) – Semi-major axis value(s).

Returns:

Eccentricity value(s) corresponding to the input semi-major axis/axes.

Return type:

ndarray

plot(a_vals: ndarray[tuple[int, ...], dtype[float64]] | None = None, num_points: int = 500, save: bool = False, filename: str | None = None, figsize: Tuple[int, int] = (8, 6), ax: Axes | None = None, show: bool = True, **plot_kwargs)

Plot the eccentricity profile e(a) with flexible matplotlib customization.

Parameters:

• (array-like (a_vals) – If None, generate new ones. Defaults to None.

• optional) (Whether to display the plot. Defaults to True.) – If None, generate new ones. Defaults to None.

• (int (num_points) – Defaults to 500.

• optional) – Defaults to 500.

• (bool (show)

• optional)

• (str (filename) – will be generated. Defaults to None.

• optional) – will be generated. Defaults to None.

• (tuple (figsize) – Defaults to (8, 6).

• optional) – Defaults to (8, 6).

• (matplotlib.axes.Axes (ax) – a new figure and axes will be created. Defaults to None.

• optional) – a new figure and axes will be created. Defaults to None.

• (bool

• optional)

• **plot_kwargs (Additional keyword arguments passed to plt.plot() and ax.set().) –

  Examples include:

  • color: Color of the line

  • linestyle: Style of the line (‘-’, ‘–’, ‘-.’, ‘:’)

  • linewidth or lw: Width of the line

  • marker: Point marker style (‘o’, ‘s’, ‘^’, etc.)

  • alpha: Transparency of the line

  • label: Label for the legend

  • log (bool): Whether to use a logarithmic scale for the y-axis.

  • xlim (tuple): Limits for the x-axis.

  • ylim (tuple): Limits for the y-axis.

  • xlabel (str): Custom x-axis label.

  • ylabel (str): Custom y-axis label.

  • title (str): Plot title.

  • grid (bool): Whether to show the grid.

Raises:

ValueError – If save=True but no filename is provided.:

debrispy.geometry module

debrispy.geometry.apoapsis(a, e)

Compute apoapsis distance, r_a.

debrispy.geometry.kappa(r, a)

Compute kappa.

debrispy.geometry.periapsis(a, e)

Compute periapsis distance, r_p.

debrispy.geometry.r(a, e, phi)

Compute radius r given semi-major axis a, eccentricity e, and true anomaly phi.

debrispy.kernel module

class debrispy.kernel.Kernel(eccentricity_profile: EccentricityProfile, r_min: float, r_max: float, num_a_points: int = 500, num_r_points: int = 500)

Bases: object

Kernel class for easy conversion of a defined eccentrcity profile into the eccentricity kernel required for calculating the ASD.

compute(eps: float = 1e-08, tol: float = 1e-10, method: str = 'gauss', rayleigh_approx: bool = False, adaptive_grid: bool = False, upper_limit: float | Callable | None = None, interpolation_method: str = 'linear', adaptive_integration: bool = False, split_points: list | None = None, n_points: int = 64, max_level: int = 25, n_jobs: int = 4) → None

Compute the kernel for the eccentricity distribution

Parameters:

• eps (float) – Epsilon precision parameter for the integration.

• tol (float) – Tolerance parameter for the integration.

• method (str) – Integration method to use. Options: ‘gauss’ (Gauss-Legendre Quadrature), ‘trapz’ (NumPy Trapezium), ‘quad’ (SciPy Quad Library)

• rayleigh_approx (bool) – Whether to use the Rayleigh approximation for the integration. (Only used for RayleighEccentricity)

• adaptive_grid (bool) – For Gauss-Legendre Quadrature: Whether to use an adaptive grid for the integration.

• upper_limit (float or callable) – Upper limit for the integration.

• interpolation_method (str) – Interpolation method to use. Options: ‘linear’ (linear interpolation), ‘cubic’ (cubic interpolation), ‘nearest’ (nearest neighbor interpolation)

• adaptive_integration (bool) – For Gauss-Legendre Quadrature: Whether to use an adaptive integration for the integration.

• split_points (list) – For Gauss-Legendre Quadrature: List of split points for the integration.

• n_points (int) – For Gauss-Legendre Quadrature: Number of points for the integration.

• max_level (int) – For Gauss-Legendre Quadrature: Maximum level for the integration.

• n_jobs (int) – Number of jobs to run in parallel.

Raises:

ValueError – If the integration method is not valid.:

compute_grad() → None

Compute the gradient of Phi(r,a) for a unique eccentricity distribution e = e(a), for the initialised grid. Uses the chain rule to compute the gradient.

get_values(a_vals, r_vals) → ndarray

num_adaptive_points() → int

Return the number of points in the adaptive Phi mesh.

phi_grid() → ndarray

Returns phi_grid if it has been computed.

phi_samples() → ndarray

Returns phi_samples if it has been computed.

plot(cmap: str = 'viridis', vmin: float | None = None, vmax: float | None = None, a_lim: tuple | None = None, r_lim: tuple | None = None, save: bool = False, filename: str | None = None, shading: str = 'auto', show_edges: bool = True, edgecolor: str = 'k', linewidth: float = 0.2, points: bool = False, point_size: float = 10) → None

Main plotting function for the Kernel.

This function plots the Phi(r,a) grid, or the Phi(r,a) samples, based on the chosen grid values.

Parameters:

• cmap (str) – The colormap to use.

• vmin (float) – The minimum value of the colormap.

• vmax (float) – The maximum value of the colormap.

• a_lim (tuple) – The limits of the semi-major axis.

• r_lim (tuple) – The limits of the radius.

• save (bool) – Whether to save the plot.

• filename (str) – The filename to save the plot to.

• shading (str) – The shading to use.

• show_edges (bool) – Whether to show the edges when using triangulation (unstructured grid).

• edgecolor (str) – The color of the edges (unstructured grid)

• linewidth (float) – The width of the edge (unstructured grid)

• points (bool) – Whether to plot a scatter plot instead of a color mesh.

• point_size (float) – The size of the points in the scatter plot.

plot_grad(type: str = 'norm', vmin: float | None = None, vmax: float | None = None, a_lim: tuple | None = None, r_lim: tuple | None = None, cmap: str = 'viridis', save: bool = False, filename: str | None = None) → None

Plot the gradient of the kernel.

Parameters:

• type (str) – The type of gradient to plot (‘norm’, ‘da’, ‘dr’, or ‘all’).

• vmin (float) – The minimum value of the gradient to plot.

• vmax (float) – The maximum value of the gradient to plot.

• a_lim (tuple) – The limits of the a-axis.

• r_lim (tuple) – The limits of the r-axis.

• cmap (str) – The colormap to use.

• save (bool) – Whether to save the plot.

• filename (str) – The filename to save the plot to.

plot_phi_kappa(cmap: str = 'viridis', save: bool = False, filename: str | None = None, a_slice: float | None = None) → None

Helper function to plot Phi(kappa, a) using sorted kappa and a 2D color plot, plus an optional secondary 1D plot of Phi(kappa) at a given a_slice (default: middle of a_grid).

Parameters:

• cmap (str) – The colormap to use.

• save (bool) – Whether to save the plot.

• filename (str) – Filename to save to (required if save=True).

• a_slice (float) – The a-value to slice at (default: middle of a_grid).

plot_slice(fix_a: float | None = None, fix_r: float | None = None, log_y: bool = False, log_x: bool = False, save: bool = False, filename: str | None = None, x_lim: tuple | None = None, y_lim: tuple | None = None) → None

Plot a 1D marginal slice of Phi(r, a) at fixed a or fixed r.

Parameters:

• fix_a (float) – Value of a at which to fix and vary r.

• fix_r (float) – Value of r at which to fix and vary a.

• log_y (bool) – Whether to plot the y-axis on a logarithmic scale.

• log_x (bool) – Whether to plot the x-axis on a logarithmic scale.

• save (bool) – Whether to save the plot.

• filename (str) – Filename to save to (required if save=True).

• x_lim (tuple) – Limits of the x-axis.

• y_lim (tuple) – Limits of the y-axis.

debrispy.kernel.compute_phi_row_gauss(i: int, a_val: float, r_grid: ndarray, psi_func: Callable, n_points: int = 64, eps: float = 1e-08, upper_limit: float | Callable | None = None, split_points: list[float | Callable] | None = None)

Compute Phi(r,a) for a general eccentricity distribution Psi(e,a) Piecewise, fixed-order Gauss-Legendre integration for a single row of the Phi(r,a) kernel.

Parameters:

• i (int) – Row index corresponding to a_val.

• a_val (float) – Semi-major axis value at this row.

• r_grid (np.ndarray) – Grid of radius values r.

• psi_func (Callable) – Function ψ(e, a) returning the eccentricity distribution.

• n_points (int) – Number of Gauss-Legendre quadrature points.

• eps (float) – Small epsilon to avoid sqrt singularities.

• upper_limit (float or callable (optional)) – Function of a (i.e., lambda a: … ) returning the upper limit of integration. If None, uses fixed upper limit of 1.0.

• split_points (list of float or callables, optional) – Breakpoints for piecewise integration. Each element may be a number or a function of a_val; only those strictly between e_min and e_max are used, in ascending order.

Returns:

• i (int) – The same row index passed in.

• phi_row (ndarray) – The computed Φ values for this a_val across r_grid.

debrispy.kernel.compute_phi_row_gauss_adaptive(i, a_val, r_grid, psi_func, n_points=64, eps=1e-08, upper_limit=None, tol=1e-10, max_level=25)

Adaptive Gauss-Legendre integration for one row of Phi(r,a).

This is particularly useful when there are discontinuities in the eccentricity distribution, but the user does not know where they exactly are (i.e. cannot provide split points).

Recursively splits [e_min,e_max] until the relative error tolerance is met or the maximum recursion depth is reached.

Parameters:

• i (int) – row index

• a_val (float) – semi-major axis

• r_grid (array_like) – radii to evaluate

• psi_func (callable) – psi(e,a) → eccentricity distribution

• n_points (int) – base number of Gauss-Legendre quadrature points

• eps (float) – small offset from kappa to avoid singularity

• upper_limit (None or callable) – if callable, upper_limit(a_val) → e_max; else e_max=1

• tol (float) – relative error tolerance (Default is 1e-10)

• max_level (int) – recursion depth limit (Default is 25)

Returns:

• i (int) – row index

• phi_row (np.ndarray) – computed Phi(r,a) values for this row

debrispy.kernel.compute_phi_row_quad(i: int, a_val: float, r_grid: ndarray, psi_func: Callable, eps: float = 1e-08, upper_limit: Callable | None = None)

Compute Phi(r,a) for an entire row using scipy.integrate.quad.

Parameters:

• i (int) – The index of the a_val in the a_grid.

• a_val (float) – The value of a to compute Phi(r,a) for (row)

• r_grid (np.ndarray) – The r_grid to compute Phi(r,a) for.

• psi_func (callable) – psi(e, a) → eccentricity distribution. The first argument is eccentricity (float), the second is semi-major axis (float).

• eps (float) – Small epsilon to avoid sqrt singularities.

• upper_limit (callable or None) – Function of a (i.e., lambda a: …) returning the upper limit of integration. If None, uses fixed upper limit of 1.0.

Returns:

phi_row – The computed Phi(r, a) values for the entire row.

Return type:

np.ndarray

debrispy.kernel.compute_phi_row_trapz(i: int, a_val: float, r_grid: ndarray, e_grid: ndarray, psi_col: ndarray, eps: float = 1e-05) → Tuple[int, ndarray]

Compute Phi(r,a) for a general eccentricity distribution Psi(e,a) This method is used to compute a single row of Phi(r,a) for a given a_val using the trapezium rule.

Parameters:

iint

The index of the a_val in the a_grid.

a_valfloat

The value of a to compute Phi(r,a) for (row)

r_gridnp.ndarray

The r_grid to compute Phi(r,a) for.

e_gridnp.ndarray

The eccentricity grid to compute Phi(r,a) from (integrate over).

psi_colnp.ndarray

The eccentricity distribution Psi(e,a) to compute Phi(r,a) from.

epsfloat

Epsilon precision parameter for the integration. (Default is 1e-5, lower is more accurate but can lead to aliasing/numerical instability when using the trapezium rule)

debrispy.kernel.compute_phi_single_gauss(pt: tuple, n_points: int, eps: float, psi_func: Callable, upper_limit: Callable | None, split_points: list[float | Callable] | None)

Compute Phi(r,a) for a single point using fixed-order Gauss-Legendre quadrature with piecewise integration. This is used for fixed-order Gauss-Legendre integration as part of the adaptive grid method.

Parameters:

• pt (tuple) – (a, r)

• n_points (int) – Number of Gauss-Legendre quadrature points per subinterval.

• eps (float) – Small epsilon to avoid sqrt singularities and overlapping bounds.

• psi_func (callable) – psi(e,a) → eccentricity distribution

• upper_limit (None or callable) – if callable, upper_limit(a) → e_max; else e_max=1

• split_points (list of float or callables, optional) – Breakpoints for piecewise integration. Each element may be a number or a function of a_val; only those strictly between e_min and e_max are used, in ascending order.

debrispy.kernel.compute_phi_single_gauss_adaptive(pt: tuple, n_points: int, eps: float, psi_func: Callable, upper_limit: float | Callable | None, tol: float = 1e-10, max_level: int = 25)

Compute Phi(r,a) for a single point using adaptive Gauss-Legendre quadrature. This is used for adaptive integration when using the adaptive grid method.

Parameters:

• pt (tuple) – (a, r)

• n_points (int) – Base number of Gauss-Legendre quadrature points.

• eps (float) – Small epsilon to avoid sqrt singularities.

• psi_func (callable) – psi(e,a) → eccentricity distribution

• upper_limit (None or callable) – if callable, upper_limit(a) → e_max; else e_max=1

• tol (float) – Relative error tolerance (Default is 1e-10).

• max_level (int) – Recursion depth limit (Default is 25).

debrispy.kernel.compute_phi_single_quad(pt: Tuple[float, float], eps: float, psi_func: Callable, upper_limit: Callable | None) → float

Compute Phi(r,a) for a single point using scipy.integrate.quad.

Parameters:

• pt (tuple) – (a, r)

• eps (float) – Small epsilon to avoid sqrt singularities.

• psi_func (callable) – psi(e, a) → eccentricity distribution. The first argument is eccentricity (float), the second is semi-major axis (float).

• upper_limit (callable or None) – Function of a (i.e., lambda a: …) returning the upper limit of integration. If None, uses fixed upper limit of 1.0.

Returns:

phi_val – The computed Phi(r, a) value.

Return type:

float

debrispy.montecarlo module

class debrispy.montecarlo.Histogram1D(edges: ndarray, values: ndarray, kind: Literal['a', 'r'], scaled: bool)

Bases: object

Container for a 1D histogram from MonteCarlo sampling.

edges

Bin edges, shape (N+1,).

Type:

np.ndarray

values

Histogram values after surface-density normalisation, shape (N,).

Type:

np.ndarray

kind

‘a’ = semi-major axis Sigma_a(a), ‘r’ = radial (ASD).

Type:

{‘a’, ‘r’}

scaled

True if the histogram was scaled to match the true area under Sigma_a(a), False if no scaling applied.

Type:

bool

property centers: ndarray

Bin centers.

edges: ndarray

get_values() → Tuple[ndarray, ndarray, ndarray]

(edges, centers, centre values) for plotting.

kind: Literal['a', 'r']

plot(ax=None, label=None, color=None, linestyle='-', **kwargs)

Plot the histogram.

Parameters:

• ax (matplotlib.axes.Axes, optional) – Axes to plot on. If None, a new figure is created.

• label (str, optional) – Label for the plot.

• color (str, optional) – Color for the plot.

• linestyle (str, optional) – Linestyle for the plot.

• **kwargs (dict, optional) – Additional keyword arguments for the plot.

Returns:

ax – Axes object with the plot.

Return type:

matplotlib.axes.Axes

scaled: bool

values: ndarray

property widths: ndarray

Bin widths.

class debrispy.montecarlo.Histogram2D(x_edges: ndarray, y_edges: ndarray, values: ndarray, mode: Literal['cartesian', 'polar'])

Bases: object

Container for a 2D histogram.

x_edges

Bin edges along x (Cartesian) or r (polar), shape (Nx+1,).

Type:

np.ndarray

y_edges

Bin edges along y (Cartesian) or phi (polar), shape (Ny+1,).

Type:

np.ndarray

values

2D array of histogram values, shape (Ny, Nx), suitable for:

plt.pcolormesh(x_edges, y_edges, values)

Type:

np.ndarray

mode

Coordinate system.

Type:

{‘cartesian’, ‘polar’}

convolve_gaussian(*, fwhm_x: float | None = None, fwhm_y: float | None = None, sigma_x: float | None = None, sigma_y: float | None = None, theta: float = 0.0, pad: float = 5.0) → Histogram2D

Convolve with a rotated Gaussian PSF.

Parameters:

• fwhm_x (float, optional) – FWHM in axis units If only fwhm_x is given, use circular PSF (fwhm_y = fwhm_x).

• fwhm_y (float, optional) – FWHM in axis units If only fwhm_x is given, use circular PSF (fwhm_y = fwhm_x).

• sigma_x (float, optional) – Sigma in axis units Mutually exclusive with FWHM. If only sigma_x is given, use circular PSF.

• sigma_y (float, optional) – Sigma in axis units Mutually exclusive with FWHM. If only sigma_x is given, use circular PSF.

• theta (float) – Rotation angle (radians, CCW).

• pad (float) – Padding margin to retain PSF wings.

Returns:

Convolved histogram.

Return type:

Histogram2D

get_values() → Tuple[ndarray, ndarray, ndarray]

(edges, centers, centre values) for plotting.

mode: Literal['cartesian', 'polar']

pad_to_limits(xlim: Tuple[float, float] | None = None, ylim: Tuple[float, float] | None = None, floor_value: float = 0.0) → Histogram2D

Return a histogram whose edges cover (xlim, ylim). Outside the original extent we create new bins (same mean bin size) and fill them with floor_value. If limits lie inside, we crop.

property phi_edges: ndarray

Alias for y_edges when mode=’polar’.

plot(*, log: bool = False, cmap: str = 'magma', shading: str = 'auto', vmin=None, vmax=None, floor_threshold=None, floor_value=None, xlim: Tuple[float, float] | None = None, ylim: Tuple[float, float] | None = None, ax=None, colorbar: bool = True, cbar_label: str = 'Counts per pixel', show_psf: bool = False, psf_scale: float = 1.0, psf_loc: Tuple[float, float] = (0.12, 0.12), psf_facecolor: str = 'white', psf_edgecolor: str = 'black', psf_alpha: float = 0.9, save: bool = False, filepath: str | None = None)

Plot with optional xlim/ylim. Regions outside the histogram are binned with same bin size and filled with floor_value (default=0).

Parameters:

• log (bool, optional) – If True, use a logarithmic scale.

• cmap (str, optional) – Colormap to use.

• shading (str, optional) – Shading to use.

• vmin (float, optional) – Minimum value to use for the colorbar.

• vmax (float, optional) – Maximum value to use for the colorbar.

• floor_threshold (float, optional) – Threshold value to use for flooring.

• floor_value (float, optional) – Value to use for flooring.

• xlim (tuple, optional) – x-axis limits.

• ylim (tuple, optional) – y-axis limits.

• ax (matplotlib.axes.Axes, optional) – Axes to plot on.

• colorbar (bool, optional) – If True, show the colorbar.

• cbar_label (str, optional) – Label for the colorbar.

• show_psf (bool, optional) – If True, show the PSF.

• psf_scale (float, optional) – Scale factor for the PSF.

• psf_loc (tuple, optional) – Location of the PSF.

• psf_facecolor (str, optional) – Facecolor of the PSF.

• psf_edgecolor (str, optional) – Edgecolor of the PSF.

• psf_alpha (float, optional) – Alpha of the PSF.

• save (bool, optional) – If True, save the figure.

• filepath (str, optional) – Filepath to save the figure.

Returns:

ax – Axes object with the plot.

Return type:

matplotlib.axes.Axes

property r_edges: ndarray

Alias for x_edges when mode=’polar’.

values: ndarray

x_edges: ndarray

y_edges: ndarray

class debrispy.montecarlo.MonteCarlo(sigma_a, ecc_profile, n_samples: int = 10000000)

Bases: object

Monte Carlo sampler for generating particle positions in a debris disc.

This class generates random samples of semi-major axis a, eccentricity e, and true anomaly f, then computes the corresponding radial positions r using orbital mechanics. The sampling is based on a given surface density profile with respect to semi-major axis (sigma_a) and an eccentricity profile (either unique or a function of semi-major axis)

sigma_a

The surface density profile used to sample semi-major axis values.

Type:

SigmaA

ecc_profile

The eccentricity profile used to sample eccentricities (can be unique or a function of semi-major axis)

Type:

EccentricityProfile

n_samples

The total number of Monte Carlo particles to generate.

Type:

int

a_samples

Cached array of sampled semi-major axis values after sampling.

Type:

np.ndarray or None

r_samples

Cached array of radial positions computed from a, e, f.

Type:

np.ndarray or None

e_samples

Cached array of eccentricity values if manually supplied or reused.

Type:

np.ndarray or None

f_samples

Cached array of true anomalies used in sampling.

Type:

np.ndarray or None

get_1d_histogram(bins: int = 500, scale: bool = True, verbose: bool = True)

Compute the 1D histogram of semi-major axis and radial positions.

This method computes the 1D histogram of semi-major axis and radial positions, optionally scaling the histogram to match the true area under the surface density profile.

Parameters:

• bins (int, optional) – Number of bins for the histogram.

• scale (bool, optional) – Whether to scale the histogram to match the true area under the surface density profile.

• verbose (bool, optional) – Whether to print progress messages.

Returns:

• histA (Histogram1D) – 1D Histogram object of semi-major axis values.

• histR (Histogram1D) – 1D Histogram object of radial positions.

get_cart_histogram(bins=500, varpi_func=None, verbose: bool = True, *, surface_density: bool = True)

Return a 2D histogram and edges in Cartesian (x, y) coordinates.

Returns:

hist_cart – 2D Histogram object in Cartesian coordinates (values shape: Ny x Nx).

Return type:

Histogram2D

get_polar_histogram(bins=500, varpi_func=None, verbose: bool = True, *, surface_density: bool = True)

Return a 2D histogram on a polar (r, phi) grid.

Returns:

hist_polar – 2D Histogram object in polar coordinates.

Return type:

Histogram2D

plot_1d(bins: int = 500, save: bool = False, filepath: str | None = None, overlay: bool = False, scale: bool = True, asd=None, x_lim: tuple[float, float] | None = None, y_lim: tuple[float, float] | None = None)

Plot the 1D histogram of semi-major axis and radial positions.

Parameters:

• bins (int, optional) – Number of bins for the histogram.

• save (bool, optional) – Whether to save the figure.

• filepath (str, optional) – Path to save the figure.

• overlay (bool, optional) – Whether to overlay the histogram with the analytic ASD.

• scale (bool, optional) – Whether to scale the histogram to match the true area under the surface density profile.

• asd (ASD, optional) – ASD object to use for the overlay.

• x_lim (tuple[float, float], optional) – Limits for the x-axis.

• y_lim (tuple[float, float], optional) – Limits for the y-axis.

plot_2d(varpi_func=None, bins=500, log=True, mode='cartesian', save=False, filepath=None, surface_density=True, **plot_kwargs)

Thin wrapper around Histogram2D.plot(). Extra kwargs are forwarded to Histogram2D.plot (e.g., cmap, shading, vmin, vmax, colorbar=False).

sample_a(use_jacobian: bool = True) → ndarray

Sample semi-major axis values from the surface density profile. This function uses batched and vectorised rejection sampling.

Parameters:

use_jacobian (bool, optional) – Whether to use the Jacobian of the surface density profile in the sampling process. If True, the sampling is weighted by the product of the surface density and semi-major axis: Sigma(a)*a If False, the sampling is uniform in the surface density: Sigma(a)

Returns:

a_samples – Array of sampled semi-major axis values.

Return type:

np.ndarray

sample_eccentricities(a_samples: ndarray) → ndarray

Sample eccentricities using proper rejection sampling conditioned on each input semi-major axis a_i.

Parameters:

a_samples (np.ndarray) – Semi-major axis values. Each e_i will be drawn from Psi(e | a_i).

Returns:

e_samples – Eccentricity values corresponding to each a_i.

Return type:

np.ndarray

sampler(use_jacobian: bool = True, verbose: bool = True, return_samples: bool = True) → tuple[ndarray, ndarray, ndarray, ndarray]

Perform the Monte Carlo sampling of semi-major axis, eccentricities, and true anomalies, and then compute the corresponding radial positions.

This method orchestrates the entire sampling process, including: - Sampling semi-major axis values - Sampling (or calculating) eccentricities - Solving Kepler’s equation for the eccentric anomaly - Computing radial positions

Parameters:

• use_jacobian (bool, optional) – Whether to use the Jacobian of the surface density profile in the sampling process. If True, the sampling is weighted by the product of the surface density and semi-major axis: Sigma(a)*a If False, the sampling is uniform in the surface density: Sigma(a)

• verbose (bool, optional) – Whether to print progress messages.

• return_samples (bool, optional) – Whether to return the samples. If False, the samples are cached internally but not returned directly.

Returns:

• a_samples (np.ndarray) – Array of sampled semi-major axis values.

• r_samples (np.ndarray) – Array of radial positions computed from a, e, f.

• e_samples (np.ndarray) – Array of eccentricities corresponding to the given semi-major axes.

• f_samples (np.ndarray) – Array of true anomalies corresponding to the given semi-major axes and eccentricities.

debrispy.montecarlo.kepler_solver(M: float, e: float, tol: float = 1e-10, max_iter: int = 100) → ndarray

Solve Kepler’s equation for the eccentric anomaly E given the mean anomaly M and eccentricity e. This function uses Newton-Raphson method to solve the equation.

Parameters:

• M (float or array-like) – Mean anomaly.

• e (float or array-like) – Eccentricity.

• tol (float, optional) – Tolerance for the solution.

• max_iter (int, optional) – Maximum number of iterations.

Returns:

E – Eccentric anomaly.

Return type:

float or array-like

debrispy.sigma_a module

class debrispy.sigma_a.SigmaA(a_min: float, a_max: float, sigma0: float = 1.0, profile_type: str = 'power_law', **kwargs)

Bases: object

Represents a surface density profile, Sigma(a), as a function of semi-major axis.

a_min(float)

Type:

Minimum value of the semi-major axis range.

a_max(float)

Type:

Maximum value of the semi-major axis range.

sigma0(float)

Type:

Base amplitude (normalization) of the profile. Defaults to 1.0.

profile_type(str)

Type:

Type of profile (‘power_law’, ‘gaussian’, ‘step_up’, ‘step_down’, ‘custom’).

VALID_PROFILES = ['power_law', 'gaussian', 'step_up', 'step_down', 'custom']

compute_area() → float

Compute the integral Sigma(a) over [a_min, a_max].

Returns:

float

Return type:

Total area under the surface density curve.

get_values(a_vals: float | ndarray[tuple[int, ...], dtype[float64]]) → ndarray[tuple[int, ...], dtype[float64]]

Calculate surface density at the given semi-major axis values.

Parameters:

array-like) (a_vals (float or)

Returns:

ndarray

Return type:

Surface density values at the specified semi-major axis values.

plot(a_vals: ndarray[tuple[int, ...], dtype[float64]] | None = None, num_points: int = 500, save: bool = False, filename: str | None = None, figsize: Tuple[int, int] = (8, 6), show: bool = True, ax: Axes | None = None, **plot_kwargs) → Tuple[Figure, Axes]

Plot the surface density profile with flexible matplotlib customization.

Parameters:

• (array-like (a_vals) – If None, use cached values or generate new ones. Defaults to None.

• optional) (Existing axes to plot on. If None,) – If None, use cached values or generate new ones. Defaults to None.

• (int (num_points) – Defaults to 500.

• optional) – Defaults to 500.

• (bool (show)

• optional)

• (str (filename) – will be generated. Defaults to None.

• optional) – will be generated. Defaults to None.

• (tuple (figsize) – Defaults to (8, 6).

• optional) – Defaults to (8, 6).

• (bool – Defaults to True.

• optional) – Defaults to True.

• (matplotlib.axes.Axes (ax) – a new figure and axes will be created. Defaults to None.

• optional) – a new figure and axes will be created. Defaults to None.

• **plot_kwargs (Additional keyword arguments passed to plt.plot().) –

  Examples include:

  • color: Color of the line

  • linestyle: Style of the line (‘-’, ‘–’, ‘-.’, ‘:’)

  • linewidth or lw: Width of the line

  • marker: Point marker style (‘o’, ‘s’, ‘^’, etc.)

  • alpha: Transparency of the line

  • label: Label for the legend.

Returns:

tuple

Return type:

Figure and axes objects for further customization if needed.

Module contents

DebrisPy: A Python package to compute the azimuthally averaged surface density (ASD) of debris discs given an eccentricity profile.