Defined Methods

ρ(d::AbstractDensity, r::Real)
ρ(uu::Unitful.DensityUnits, d::AbstractDensity, r::Real)
ρ(d::AbstractDensity, r::Unitful.Length)
ρ(uu::Unitful.DensityUnits, d::AbstractDensity, r::Unitful.Length)

Returns the density of profile d [M⊙ kpc^-3] at radius r [kpc].

invρ(d::AbstractDensity, x::Real) # Only valid if specialized method exists
invρ(d::AbstractDensity, x::Real,
     interval::NTuple{2,Real}=(scale_radius(d)/100,
                               100*scale_radius(d)); kws...) # Root-finding generic method
invρ(uu::Unitful.LengthUnits, d::AbstractDensity, x::Real)
invρ(d::AbstractDensity, x::Unitful.Density)
invρ(uu::Unitful.LengthUnits, d::AbstractDensity, x::Unitful.Density)

Solve for the radius r [kpc] at which the density [M⊙ kpc^-3] is x for density profile d. Requires x>0.

When this method is not specialized for d, it will use an interval bracketing method from Roots.jl, requiring that ρ(d,r) be defined. For this method, kws... are passed to Roots.find_zero.

∇ρ(d::AbstractDensity, r::Real)
∇ρ(uu::GalaxyProfiles.∇ρdimensionUnits, d::AbstractDensity, r::Real)
∇ρ(d::AbstractDensity, r::Unitful.Length)
∇ρ(uu::GalaxyProfiles.∇ρdimensionUnits, d::AbstractDensity, r::Unitful.Length)

Returns the radial derivative of the density [M⊙ kpc^-4] of the mass profile d evaluated at radius r [kpc].

ρmean(d::AbstractDensity, r::Real)
ρmean(uu::Unitful.DensityUnits, d::AbstractDensity, r::Real)
ρmean(d::AbstractDensity, r::Unitful.Length)
ρmean(uu::Unitful.DensityUnits, d::AbstractDensity, r::Unitful.Length)

Returns the average density [M⊙ kpc^-3] of density profile d inside r [kpc]. The generic method uses enclosed mass divided by spherical volume; M(d,r) * 3 / (4π*r^3).

invρmean(d::AbstractDensity, x::Real)
invρmean(d::AbstractDensity, x::Real,
         interval::NTuple{2,Real}=(scale_radius(d)/100,
                                   100*scale_radius(d)); kws...)
invρmean(uu::Unitful.LengthUnits, d::AbstractDensity, x::Real)
invρmean(d::AbstractDensity, x::Unitful.Density)
invρmean(uu::Unitful.LengthUnits, d::AbstractDensity, x::Unitful.Density)

Solve for the radius r [kpc] inside which the average density is x [M⊙ kpc^-3]. Requires x>0.

When this method is not specialized for d, it will use an interval bracketing method from Roots.jl, requiring that ρmean(d,r) or M(d,r) be defined. For this method, kws... are passed to Roots.find_zero.

Σ(d::AbstractMassProfile, r::Real)
Σ(uu::GalaxyProfiles.SurfaceDensityUnits, d::AbstractMassProfile, r::Real)
Σ(d::AbstractMassProfile, r::Unitful.Length)
Σ(uu::GalaxyProfiles.SurfaceDensityUnits, d::AbstractMassProfile, r::Unitful.Length)

Returns the surface density [M⊙ kpc^-2] of mass profile d at radius r. For 3D density profiles (i.e., subtypes of AbstractDensity), this will be the projected surface density, which, for spherical systems, is defined by the Abel integral

\[\Sigma(r) = 2 \int_R^\infty \rho(r) \frac{r}{ \sqrt{r^2-R^2} } \, dr\]

which has an inverse of

\[\rho(r) = -\frac{1}{\pi} \int_r^\infty \frac{d\Sigma(R)}{dR} \frac{dR}{\sqrt{R^2-r^2}}\]

invΣ(d::AbstractMassProfile, x::Real)
invΣ(d::AbstractMassProfile, x::Real,
     interval::NTuple{2,Real}=(scale_radius(d)/100,
                               100*scale_radius(d)); kws...)
invΣ(uu::Unitful.LengthUnits, d::AbstractMassProfile, x::Real)
invΣ(d::AbstractMassProfile, r::Unitful.Length)
invΣ(uu::Unitful.LengthUnits, d::AbstractMassProfile, r::Unitful.Length)

Solve for the radius r [kpc] at which the surface density [M⊙ kpc^-2] is x for profile d. Requires x>0.

When this method is not specialized for d, it will use an interval bracketing method from Roots.jl, requiring that Σ(d,r) be defined. For this method, kws... are passed to Roots.find_zero.

∇Σ(d::AbstractMassProfile, r::Real)
∇Σ(uu::Unitful.DensityUnits, d::AbstractMassProfile, r::Real)
∇Σ(d::AbstractMassProfile, r::Unitful.Length)
∇Σ(uu::Unitful.DensityUnits, d::AbstractMassProfile, r::Unitful.Length)

Returns the radial derivative of the surface density [M⊙ kpc^-3] of the mass profile d evaluated at radius r [kpc].

Σmean(d::AbstractMassProfile, r::Real)
Σmean(uu::GalaxyProfiles.SurfaceDensityUnits, d::AbstractMassProfile, r::Real)
Σmean(d::AbstractMassProfile, r::Unitful.Length)
Σmean(uu::GalaxyProfiles.SurfaceDensityUnits, d::AbstractMassProfile, r::Unitful.Length)

Returns the mean projected surface density [M⊙ kpc^-2] of the mass profile d inside the radius r [kpc]. The generic method uses projected enclosed mass divided by area: Mproj(d::AbstractMassProfile,r::Real) / (π * r^2).

M(d::AbstractDensity, r::Real)
M(uu::Unitful.MassUnits, d::AbstractDensity, r::Real)
M(d::AbstractDensity, r::Unitful.Length)
M(uu::Unitful.MassUnits, d::AbstractDensity, r::Unitful.Length)

Returns the total mass [M⊙] of the profile d enclosed within a spherical shell of radius r. For spherical systems this is given by the integral

\[M(\lt R) = 4\pi \int_0^R r^2 ρ(r) dr\]

When this method is not specialized for d, it will compute the numerical integral using QuadGK.quadgk.

invM(d::AbstractDensity, x::Real)
invM(d::AbstractDensity, x::Real,
     interval::NTuple{2,Real}=(scale_radius(d)/100,
                               100*scale_radius(d)); kws...)
invM(uu::Unitful.LengthUnits, d::AbstractDensity, x::Real)
invM(d::AbstractDensity, x::Unitful.Mass)
invM(uu::Unitful.LengthUnits, d::AbstractDensity, x::Unitful.Mass)

Solve for the radius r [kpc] at which the enclosed mass [M⊙] is x for profile d. Requires x>0.

When this method is not specialized for d, it will use an interval bracketing method from Roots.jl, requiring that M(d,r) be defined. For this method, kws... are passed to Roots.find_zero.

∇M(d::AbstractDensity, r::Real)
∇M(uu::GalaxyProfiles.∇mdimensionUnits, d::AbstractDensity, r::Real)
∇M(d::AbstractDensity, r::Unitful.Length)
∇M(uu::GalaxyProfiles.∇mdimensionUnits, d::AbstractDensity, r::Unitful.Length)

The radial derivative of the enclosed mass [M⊙ kpc^-1] of profile d evaluated at radius r.

Mtot(d::AbstractMassProfile)
Mtot(uu::Unitful.MassUnits, d::AbstractMassProfile)

Returns the total mass [kpc] of the profile d. For profiles which do not have a well-defined total mass M(d,Inf), this quantity should be defined as the limit $\lim_{r \to \infty} M(r)$.

Mproj(d::AbstractMassProfile, r::Real)
Mproj(uu::Unitful.MassUnits, d::AbstractMassProfile, r::Real)
Mproj(d::AbstractMassProfile, r::Unitful.Length)
Mproj(uu::Unitful.MassUnits, d::AbstractMassProfile, r::Unitful.Length)

Returns the total line-of-sight projected mass [M⊙] enclosed within a radius r [kpc] for the profile d. For spherical systems this is given by the integral

\[M_{\text{proj}}(\lt R) = 2\pi \int_0^R r \, \Sigma(r) dr = 2\pi \int_0^R r \left( 2 \int_r^\infty \rho(r^\prime) \frac{r^\prime}{ \sqrt{r^{\prime \, 2}-r^2} } \, dr^\prime \right) dr\]

∇Mproj(d::AbstractMassProfile, r::Real)
∇Mproj(uu::GalaxyProfiles.∇mdimensionUnits, d::AbstractMassProfile, r::Real)
∇Mproj(d::AbstractMassProfile, r::Unitful.Length)
∇Mproj(uu::GalaxyProfiles.∇mdimensionUnits, d::AbstractMassProfile, r::Unitful.Length)

The radial derivative of the line-of-sight projected mass [M⊙ kpc^-1] of profile d enclosed within a radius r [kpc].

invMproj(d::AbstractMassProfile, x::Real)
invMproj(d::AbstractMassProfile, x::Real,
         interval::NTuple{2,Real}=(scale_radius(d)/100,
                                   100*scale_radius(d)); kws...)
invMproj(d::AbstractMassProfile, x::Unitful.Mass)
invMproj(uu::Unitful.LengthUnits, d::AbstractMassProfile, x::Unitful.Mass)

Solve for the radius r [kpc] at which the line-of-sight projected enclosed mass is x [M⊙] for profile d. Requires x>0.

When this method is not specialized for d, it will use an interval bracketing method from Roots.jl, requiring that M(d,r) be defined. For this method, kws... are passed to Roots.find_zero.

dynamical_time(d::Union{AbstractDensity, AbstractMassProfile}, r::Real)
dynamical_time(d::Union{AbstractDensity, AbstractMassProfile}, r::Unitful.Length)
dynamical_time(uu::Unitful.TimeUnits, d::Union{AbstractDensity, AbstractMassProfile}, r::Unitful.Length)

Returns the dynamical time [yr] at radius r [kpc] in the provided density profile d. The dynamical time is a characteristic timescale associated with orbits in potentials. This implementation specifically calculates Equation 2.39 in Binney & Tremaine Galactic Dynamics 2E, with the replacement of the standard density ρ for the average density interior to r, ρmean, to better account for inhomogenous systems.

\[t_\text{dyn} \left( r \right) = \sqrt{ \frac{3\pi}{16 \, G \, \overline{\rho} \left( r \right)} }\]

The above equation is used when the argument d::AbstractDensity. As ρmean can be expressed as a function of the total mass interior to the orbit, this can also be written as

\[\begin{aligned} \overline{\rho} \left( r \right) &= \frac{3 \, M \left( r \right)}{4 \, \pi \, r^3} \newline t_\text{dyn} \left( r \right) &= \pi \, \sqrt{ \frac{r^3}{4 \, G \, M \left( r \right)} } \end{aligned}\]

where $M \left( r \right)$ is the mass enclosed inside radius $r$, provided by the method M. This equation is used when the argument d::AbstractMassProfile, with the substitution of the projected mass Mproj for the 3-D enclosed mass M, which is undefined for profiles that do not have ρ densities.

Alternative Definitions

Note that some authors prefer to omit the factor of $\sqrt{ \frac{1}{16} } = 1/4$ in the denominator of the first equation above; this is, for example, the convention used by galpy. Dynamical times thus defined will be four times larger than those calculated by this method.

cdf2D(d::AbstractMassProfile, r::Real)
cdf2D(d::AbstractMassProfile, r::Unitful.Quantity)

Returns the value of the cumulative distribution function of the mass profile d at r [kpc] in two dimensions (i.e., projected along a line of sight). This is defined as Mproj(d,r) / Mtot(d).

cdf3D(d::AbstractDensity, r::Real)
cdf3D(d::AbstractDensity, r::Unitful.Quantity)

Returns the value of the cumulative distribution function of the profile d at r [kpc] in three dimensions. This is defined as M(d,r) / Mtot(d).

ccdf2D(d::AbstractMassProfile, r::Real)
ccdf2D(d::AbstractMassProfile, r::Unitful.Quantity)

Returns the value of the complementary cumulative distribution function of the profile d at r [kpc] in two dimensions (i.e., along a line of sight). This is defined as 1 - cdf2D(d,r) = 1 - Σ(d,r) / Mtot(d).

ccdf3D(d::AbstractMassProfile, r::Real)
ccdf3D(d::AbstractMassProfile, r::Unitful.Quantity)

Returns the value of the complementary cumulative distribution function of the profile d at r [kpc] in three dimensions. This is defined as 1 - cdf3D(d,r) = 1 - M(d,r) / Mtot(d).

quantile2D(d::AbstractMassProfile, r::Real)
quantile2D(d::AbstractMassProfile, r::Unitful.Quantity)

Returns the value of the inverse cumulative distribution function of the profile d at r [kpc] in two dimensions (i.e., along a line of sight).

quantile3D(d::AbstractDensity, r::Real)
quantile3D(d::AbstractDensity, r::Unitful.Quantity)

Returns the value of the inverse cumulative distribution function of the profile d at r [kpc] in three dimensions.

cquantile2D(d::AbstractMassProfile, r::Real)
cquantile2D(d::AbstractMassProfile, r::Unitful.Quantity)

Returns the value of the complementary quantile of the profile d at r [kpc] in two dimensions (i.e., along a line of sight). This is defined as quantile2D(d, 1-r).

cquantile3D(d::AbstractDensity, r::Real)
cquantile3D(d::AbstractDensity, r::Unitful.Quantity)

Returns the value of the complementary quantile of the profile d at r [kpc] in three dimensions. This is defined as quantile3D(d, 1-r).

Vcirc(d::AbstractDensity, r::Real)
Vcirc(uu::Unitful.VelocityUnits, d::AbstractDensity, r::Real)
Vcirc(d::AbstractDensity, r::Unitful.Length)
Vcirc(uu::Unitful.VelocityUnits, d::AbstractDensity, r::Unitful.Length)

Returns the circular velocity [km/s] at r [kpc], defined as the speed of a particle of insignificant mass in a circular orbit at radius r. This is calculated as

\[v_c^2(r) = \frac{G M(r)}{r} = r \frac{d\Phi}{dr} = r\nabla\Phi(r)\]

By default uses G in units such that if rs and r are in kpc, the velocity ends up in km/s. For example, for GeneralIsothermal we have [G] = [kpc * km^2 / Msun / s^2] so that the velocity ends up in km/s. This method has a generic implementation of sqrt( GalaxyProfiles.constants.Gvelkpc * M(d::AbstractDensity,r) / r).

Vesc(d::AbstractDensity, r::Real)
Vesc(uu::Unitful.VelocityUnits, d::AbstractDensity, r::Real)
Vesc(d::AbstractDensity, r::Unitful.Length)
Vesc(uu::Unitful.VelocityUnits, d::AbstractDensity, r::Unitful.Length)

Returns the escape velocity [km/s] at r [kpc], calculated as

\[v^2_{\text{esc}}(r) = 2 |\Phi(r)|\]

if $\Phi \to 0$ for $r \to \infty$; see the note for Φ.

By default uses G in units such that if rs and r are in kpc, the velocity ends up in km/s. For example, for GeneralIsothermal we have [G] = [kpc * km^2 / Msun / s^2] so that the velocity ends up in km/s.

Vmax(d::AbstractDensity)
Vmax(uu::Unitful.VelocityUnits, d::AbstractDensity)

Returns the maximum circular velocity [km/s] of d and the corresponding radius [kpc]. Can be found by solving

\[ \frac{d v_c(r)}{dr} = 0\]

for r, where $v_c$ is the circular velocity, then evaluating the circular velocity at r.

σr(d::AbstractDensity, r::Real, β::Real)

Returns the radial velocity dispersion [km/s] of d at radius r [kpc] for constant velocity anisotropy β given by Equation 4.216 in Binney & Tremaine Galactic Dynamics 2E,

\[\sigma_r^2 \left( R \right) = \frac{1}{R^{2\beta} \, \rho(R)} \int_R^\infty r^{2\beta} \, \rho(r) \, \frac{d\Phi}{dr} \, dr\]

and as $\frac{d\Phi}{dr} = -G M(r) / r^2$ we can alternatively write

\[\sigma_r^2 \left( R \right) = \frac{G}{R^{2\beta} \, \rho(R)} \int_R^\infty r^{2\left( \beta-1 \right)} \, \rho(r) \, M\left( r \right) \, dr\]

which is the generic method as we expect $M(r)$ to be more commonly available than $\frac{d\Phi}{dr}$.

σlos(d::AbstractDensity, r::Real, β::Real)

Returns the line-of-sight projected velocity dispersion [km/s] of d at projected radius r [kpc] for constant velocity anisotropy β given by

\[\sigma_{\text{LOS}}^2 \left( R \right) = \frac{2}{\Sigma \left( R \right)} \int_R^\infty \left(1 - \beta \frac{R^2}{r^2} \right) \frac{r \, \rho \left( r \right) \, \sigma_r^2 \left( r \right)}{\sqrt{r^2 - R^2}} \, dr\]

Φ(d::AbstractDensity, r::Real)
Φ(uu::GalaxyProfiles.ΦdimensionUnits, d::AbstractDensity, r::Real)
Φ(d::AbstractDensity, r::Unitful.Length)
Φ(uu::GalaxyProfiles.ΦdimensionUnits, d::AbstractDensity, r::Unitful.Length)

Returns the potential [km^2 s^-2] of the density distribution d at radius r [kpc].

This is typically defined as

\[\begin{aligned} \Phi(R) &= -4\pi G \left( \frac{1}{R} \int_0^R r^2 \rho(r) dr + \int_R^\infty r \rho(r) dr \right) \newline &= -G \left( \frac{M (R)}{R} - \int_R^\infty r \rho(r) dr \right) \newline % \Phi(R) &= -\frac{G}{R} \int_0^R dM(r) - G \int_R^\infty \frac{dM(r)}{r} \newline \end{aligned}\]

where $M(R)$ is the mass internal to radius $R$. The final expression is used in the generic implementation of Φ. Given that there is also a generic implementation of $M(R)$ that calculates it from the density ρ, the potential Φ can be calculated from a density profile by defining only the density ρ.

If an analytic expression for $M(r)$ exists, then an equivalent equation expressed in terms of the enclosed mass may be used,

\[\begin{aligned} \Phi(R) &= -4\pi G \left( \frac{1}{R} \int_0^R r^2 \rho(r) dr + \int_R^\infty r \rho(r) dr \right) \newline &= -G \left( \frac{M (R)}{R} - \int_R^\infty \frac{dM(r)}{r} \right) \newline &= -G \, \int_R^\infty \frac{M \left( r \right)}{r^2} dr \newline \end{aligned}\]

The above integrals are not finite for some mass distributions (e.g., GeneralIsothermal with some choices of α). In these cases, it is convention to define the potential ar r as the potential difference between r and the characteristic scale radius of the distribution; i.e.

\[\Phi(R) - \Phi(R_s) = G \int_{R_s}^R \frac{M(r)}{r^2} dr\]

such that $\Phi(R_s)\equiv0$.

∇Φ(d::AbstractDensity, r::Real)
∇Φ(uu::u.AccelerationUnits, d::AbstractDensity, r::Real)
∇Φ(d::AbstractDensity, r::Unitful.Length)
∇Φ(uu::u.AccelerationUnits, d::AbstractDensity, r::Unitful.Length)

Returns the radial derivative of the potential Φ(d,r) [km s^-2] of profile d at radius r [kpc].

For a spherical mass distribution, the gravitational acceleration term in Newton's second law $F=m*a$ is $a = -G M(r) / r^2 = -\nabla\Phi$ by definition. This is derived from the definition of Φ below.

Derivation from Φ

Recall that the potential can be written as

\[\Phi(R) = -G \, \int_R^\infty \frac{M \left( r \right)}{r^2} dr.\]

so that our derivative is

\[\begin{aligned} \nabla\Phi(R) &= \frac{\partial \Phi (R)}{\partial R} \newline &= \frac{\partial}{\partial R} \left[ -G \, \int_R^\infty \frac{M \left( r \right)}{r^2} dr \right] \newline \end{aligned}\]

We are thus taking a radial derivative of a radial integral. By applying the fundamental thereom of calculus and remembering that we desire $\lim_{R \to \infty} \Phi(R) = 0$, we can simply write

\[\begin{aligned} \nabla\Phi(R) &= \frac{G M \left( r \right)}{r^2} \end{aligned}\]

∇∇Φ(d::AbstractDensity, r::Real)
∇∇Φ(uu::GalaxyProfiles.∇∇ΦdimensionUnits, d::AbstractDensity, r::Real)
∇∇Φ(d::AbstractDensity, r::Unitful.Length)
∇∇Φ(uu::GalaxyProfiles.∇∇ΦdimensionUnits, d::AbstractDensity, r::Unitful.Length)

Returns the second radial derivative of the potential Φ(d,r) [km s^-2 kpc^-1] evaluated at radius r [kpc].

As the first radial derivative $\nabla\Phi(R) = \frac{G M \left( r \right)}{r^2}$, by the product rule the second radial derivative is

\[\nabla\nabla\Phi(R) = G \left( \frac{\nabla M(R)}{r^2} - \frac{2 M(R)}{r^3} \right)\]

which is used as the generic implementation through the methods M and ∇M.

Note that this is not the same as the Laplacian operator that appears the Poisson equation $\nabla^2 \Phi(R) = 4 \pi G \rho(R)$. In spherical coordinates, the radial component of the left hand side of Poisson's equation expands to

\[\nabla^2 \Phi(R) = \frac{1}{R^2} \frac{\partial}{\partial R} \left( R^2 \frac{\partial \Phi(R)}{\partial R} \right)\]

which is not equivalent to the second radial gradient $\frac{\partial^2 \Phi(R)}{\partial R^2}$, which is what this method returns.

plummer_a_to_rh(a::Real)

Returns the 3-D half-light (or half-mass) radius given the Plummer scale radius a. This is equivalent to quantile3D(d::Plummer, 0.5) but faster because the argument x=0.5 is known at compile-time. Note that in projection (i.e., along a line-of-sight) the half-light radius is simply equal to the Plummer scale radius a, verifiable by quantile2D(d::Plummer, 0.5) == scale_radius(d).

plummer_angular_avalue(absmag, sb, DM)

This is an observational utility designed for use with the astronomical magnitude system. Returns the Plummer scale radius a given a desired absolute magnitude absmag, average surface brightness within the half-light radius sb, and distance modulus DM. If the units of sb are, for example, [mags/arcsec^2], and the units of absmag and DM should always be magnitudes, then this will return the Plummer scale radius a in [arcsec]. Mathematically, we solve the equation for the surface brightness S as a function of magnitude m = absmag + DM and area A=πr^2, with r being the half-light radius:

\[\begin{aligned} S &= m + 2.5 \times \log(A) = m + 2.5 \times \log(πr^2) = m + 2.5 \times ( \log(π) + 2\log(r) ) \newline r &= \text{exp10} \left( \frac{S - m - 2.5 \, \log(π)}{5} \right) = \text{exp10} \left( \frac{S - m}{5} \right) \times π^{1/2} \end{aligned}\]

We then only need to convert the half-light radius r to the scale radius a, which is just the inverse of GalaxyProfiles.plummer_a_to_rh.

Examples

For a Plummer profile with a total magnitude of -5 (flux = -2.5*log10(mag) = 100), an average surface brightness within the half-light radius of 25 mag/arcsec^2, and a distance of 1 Mpc (DM = 5*log10(1e6 [pc] / 10 [pc]) = 25), we can compute the scale radius of the corresponding Plummer profile in arcseconds as

julia> isapprox( GalaxyProfiles.plummer_angular_avalue(-5, 25, 25), 4.32; rtol=1e-2)
true

sample2D(d::AbstractMassProfile, center::SkyCoords.AbstractSkyCoords, distance; 
         dims..., rng=Random.default_rng())

Sample 2D sky coordinates from the mass profile d. Requires that sample2D_r is implemented for the profile d.

Arguments

• d::AbstractMassProfile: The mass profile to sample from
• center::SkyCoords.AbstractSkyCoords: The center position of the galaxy
• distance: Distance to the galaxy in the same physical units as the profile scale radius

Keyword Arguments

• dims...: Optional dimensions for the output array (e.g., 100 for 100 samples, or (10, 10) for a 10×10 array)
• rng: Random number generator (default: Random.default_rng())

Returns

A vector or array of SkyCoords objects representing the positions of sampled points on the sky.

Details

The function works by:

1. Sampling projected radii from the surface density profile of d
2. Converting physical radii to angular separations using the distance
3. Randomly distributing points uniformly in position angle (0 to 2π)
4. Using SkyCoords.offset to compute sky coordinates relative to the galaxy center

Example

using GalaxyProfiles, SkyCoords

# Create a Plummer profile with scale radius 500 pc  
plummer = Plummer(1000.0, 0.5)

# Galaxy center at RA=270°, Dec=66° (converted to radians)
center = ICRSCoords(deg2rad(270.0), deg2rad(66.0))

# Distance to galaxy: 100 kpc
distance = 100.0

# Sample 1000 positions on the sky
coords = sample2D(plummer, center, distance, 1000)

# These coordinates are of the same type as the input `center`,
# in radians. We can separate out into RA and Dec and convert to deg2rad
# if needed.
ra, dec = rad2deg.(SkyCoords.lon.(coords)), rad2deg.(SkyCoords.lat.(coords))

sample2D_r!([rng::Random.AbstractRNG=Random.default_rng()], d::AbstractMassProfile, x::AbstractArray{<:Real})

Draw multiple samples of the radius following the surface density profile of d and write them to x in place. By default, falls back to quantile2D.(d, rand!(rng, x)). See also sample2D_r.

sample2D_r([rng::Random.AbstractRNG=Random.default_rng()], d::AbstractMassProfile [, dims...])

Draws samples of the radius following the surface density profile of d with shape [dims...], which can either be a Vararg{Integer, N} (e.g., sample2D_r(rng, d, 2, 2) to generate a 2x2 Matrix) or a Dims{N}, which is an NTuple of N Ints (e.g., sample2D_r(rng, d, (2,2)). By default, falls back to quantile2D.(d, rand(rng [, dims...])). See also sample2D_r!.

sample3D_r!([rng::Random.AbstractRNG=Random.default_rng()], d::AbstractDensity, x::AbstractArray{<:Real})

Draw multiple samples of the radius following the density profile of d and write them to x in place. By default, falls back to quantile3D.(d, rand!(rng, x)). See also sample3D_r.

sample3D_r([rng::Random.AbstractRNG=Random.default_rng()], d::AbstractDensity [, dims...])

Sample random points following the density profile of d with shape [dims...], which can either be a Vararg{Int,N} (e.g., sample3D_r(rng, d, 2, 2) to generate a 2x2 Matrix) or a Dims{N}, which is an NTuple of N Ints (e.g., sample3D_r(rng, d, (2,2)). By default, falls back to quantile3D.(d, rand(rng [, dims...])). See also sample3D_r!.