debrispy.montecarlo

Functions

kepler_solver(M, e[, tol, max_iter])

Solve Kepler's equation for the eccentric anomaly E given the mean anomaly M and eccentricity e.

Classes

`Histogram1D`(edges, values, kind, scaled)	Container for a 1D histogram from MonteCarlo sampling.
`Histogram2D`(x_edges, y_edges, values, mode)	Container for a 2D histogram.
`MonteCarlo`(sigma_a, ecc_profile[, n_samples])	Monte Carlo sampler for generating particle positions in a debris disc.

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

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][source]: (edges, centers, centre values) for plotting.

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

plot(ax=None, label=None, color=None, linestyle='-', **kwargs)[source]

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'])[source]

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[source]

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][source]: (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[source]: 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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]

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)[source]: 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[source]

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[source]

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][source]

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[source]

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