profiles.density.CoredPowerLawDensityProfile.compute_surface_density#

CoredPowerLawDensityProfile.compute_surface_density(R: unyt_array | unyt_quantity, units: str | Unit | None = None, **kwargs) → unyt_array | unyt_quantity#

Numerically compute the projected surface density at radius R from the origin.

Mathematically, this is the operation

\[\Sigma(R) = 2 \int_0^\infty \rho \left( \sqrt{R^2 + z^2} \right) dz,\]

where \(z\) refers to the displacement along the line of sight.

For some value of \(R\), this method will perform the relevant integral of the profile and provide the calculated surface density.

Parameters:

  • R (unyt_array or unyt_quantity) – The projected radius from the origin in physical units with dimension length. R may be either an array of values or a scalar.

  • units (Unit or str, optional) – The units in which to return the calculated surface density. By default, the internal units are preserved.

  • kwargs – Additional keyword arguments to pass on to quad_vec().

Returns:

sigma – Surface density at each R, with correct units.

Return type:

unyt_array or unyt_quantity

Example

import numpy as np
from pisces.profiles.density import (
    NFWDensityProfile,
)

# Create a density profile with characteristic parameters
profile = NFWDensityProfile(rho_0=1.0, r_s=10.0)

# Define projected radii (can be a scalar or array)
R = np.linspace(0.1, 100.0, 50)  # kpc

# Compute the surface density at each radius
sigma = profile.compute_surface_density(
    R, units="Msun/pc**2"
)

print(sigma)