This page was generated from a Jupyter notebook. You can download it here.

Particle-spray Stream Modeling with streamspraydf

The particle-spray technique generates tidal streams by “spraying” particles from the Lagrange points of a progenitor cluster as it orbits in a gravitational potential. Each particle is given a small velocity kick and then integrated forward in the potential. galpy provides two implementations (described in Qian et al. 2022):

• chen24spraydf: the spray model from Chen et al. (2024)

• fardal15spraydf: the spray model from Fardal et al. (2015)

These are faster to set up than streamdf and work in any potential.

%matplotlib inline
import numpy
import matplotlib.pyplot as plt
from astropy import units as u
from galpy.df import chen24spraydf, fardal15spraydf
from galpy.potential import LogarithmicHaloPotential
from galpy.orbit import Orbit
import warnings

warnings.filterwarnings("ignore", category=RuntimeWarning)
warnings.filterwarnings("ignore", category=UserWarning)

Setup

We set up a potential and a progenitor orbit. The particle-spray classes only need the potential and the progenitor orbit (no action-angle machinery).

lp = LogarithmicHaloPotential(normalize=1.0, q=0.9)
prog = Orbit(
    [1.56148083, 0.35081535, -1.15481504, 0.88719443, -0.47713334, 0.12019596],
    ro=8.0,
    vo=220.0,
)

Initialize chen24spraydf

We create a spray model for the leading tail. The key parameters are:

• progenitor_mass: mass of the progenitor in solar masses

• progenitor: the progenitor orbit

• pot: the gravitational potential

• tdisrupt: time since disruption began (in Gyr when using physical units)

• tail: which tail to model ('leading', 'trailing', or 'both')

spdf_lead = chen24spraydf(
    progenitor_mass=2 * 10**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="leading",
)

Sample stream particles

The sample method generates stream particles by spraying n particles from the progenitor’s Lagrange points and integrating them forward.

numpy.random.seed(1)
lead_stream = spdf_lead.sample(n=500)

Let’s visualize the sampled stream particles:

plt.figure(figsize=(8, 4))
plt.subplot(1, 2, 1)
plt.plot(lead_stream.R(), lead_stream.z(), "k.", ms=1)
plt.xlabel(r"$R$")
plt.ylabel(r"$z$")
plt.title("Leading arm (R, z)")
plt.subplot(1, 2, 2)
plt.plot(lead_stream.R(), lead_stream.vR(), "k.", ms=1)
plt.xlabel(r"$R$")
plt.ylabel(r"$v_R$")
plt.title("Leading arm (R, vR)")
plt.tight_layout();

Leading and trailing arms with tail=’both’

Using tail='both' generates both arms in one call. This is more efficient than creating separate leading and trailing instances.

numpy.random.seed(1)
spdf_both = chen24spraydf(
    progenitor_mass=2 * 10**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="both",
)
both_stream = spdf_both.sample(n=1000)

Plotting both arms in Cartesian coordinates:

plt.figure(figsize=(6, 6))
plt.plot(both_stream.x(), both_stream.y(), "k.", ms=1)
plt.plot(prog.x(), prog.y(), "ro", ms=5)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.title("Both tails in Cartesian coordinates")
plt.axis("equal");

Using fardal15spraydf

fardal15spraydf provides an alternative spray model based on Fardal et al. (2015). The interface is the same.

numpy.random.seed(1)
spdf_fardal = fardal15spraydf(
    progenitor_mass=2 * 10**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="both",
)
fardal_stream = spdf_fardal.sample(n=1000)

Comparing the two spray models side by side:

plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(both_stream.x(), both_stream.y(), "k.", ms=1)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.title("chen24spraydf")
plt.axis("equal")
plt.subplot(1, 2, 2)
plt.plot(fardal_stream.x(), fardal_stream.y(), "k.", ms=1)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.title("fardal15spraydf")
plt.axis("equal")
plt.tight_layout();

Progenitor self-gravity with progpot

The streamspraydf models can include the progenitor’s gravity using the progpot parameter:

numpy.random.seed(1)
from galpy.potential import PlummerPotential

progpot = PlummerPotential(2 * 10.0**4.0 * u.Msun, 4.0 * u.pc)
spdf_with_prog = chen24spraydf(
    progenitor_mass=2 * 10.0**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="both",
    progpot=progpot,
)
stream_with_prog = spdf_with_prog.sample(n=500)
plt.plot(stream_with_prog.x(), stream_with_prog.z(), "k.", ms=1)
plt.xlabel(r"$x$")
plt.ylabel(r"$z$")
plt.title("chen24spraydf with progenitor gravity");

Time-dependent progenitor mass

Real cluster progenitors lose mass over their lifetimes. The progenitor_mass argument accepts either a constant (as above) or a callable M(t). The callable receives the progenitor-time coordinate t — t=0 is now and t<0 is the past — and may take/return astropy Quantity values, which are auto-detected.

Here we compare a constant-mass run to one in which the progenitor mass declines by a factor of 10 from start of disruption to the present (so the cluster is more massive in the past and is now nearly dissolved). Earlier strippings happen at the higher mass and therefore a wider tidal radius.

from astropy import units as u


def M_declining(t):
    # t = Quantity time, t=0 now, t<0 past
    # M0 at t=-4.5 Gyr, 0.1*M0 at t=0
    return 2e4 * u.Msun * (0.1 - 0.9 * t / (4.5 * u.Gyr))


spdf_evol = chen24spraydf(
    progenitor_mass=M_declining,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="both",
)
spdf_const = chen24spraydf(
    progenitor_mass=2e4 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="both",
)

numpy.random.seed(0)
stream_evol = spdf_evol.sample(n=1000)
numpy.random.seed(0)
stream_const = spdf_const.sample(n=1000)

plt.figure(figsize=(12, 4))
plt.subplot(1, 2, 1)
plt.plot(stream_const.x(), stream_const.y(), "k.", ms=1)
plt.plot(prog.x(), prog.y(), "ro", ms=5)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.title("Constant mass")
plt.axis("equal")
plt.subplot(1, 2, 2)
plt.plot(stream_evol.x(), stream_evol.y(), "k.", ms=1)
plt.plot(prog.x(), prog.y(), "ro", ms=5)
plt.xlabel(r"$x$")
plt.ylabel(r"$y$")
plt.title("Declining mass (1.0 → 0.1 M0)")
plt.axis("equal")
plt.tight_layout();

Integrating sampled orbits

The sample method can return integrated orbits (to the present time) using integrate=True. You can also get the stripping times with returndt=True. The returned Orbit object can be further integrated or analyzed:

numpy.random.seed(1)
# Sample with integration and stripping times
orbs_integrated, dt = spdf_both.sample(n=500, returndt=True, integrate=True)

# Integrate the progenitor for context
ts_prog = numpy.linspace(0.0, 3.0, 301)
prog.integrate(ts_prog, lp)
prog.integrate(-ts_prog, lp)

# Plot the integrated stream with progenitor orbit
prog.plot(d1="R", d2="z", color="k", lw=1)
plt.plot(orbs_integrated.R(), orbs_integrated.z(), "r.", ms=1, alpha=0.5)
plt.xlabel(r"$R$")
plt.ylabel(r"$z$")
plt.title("Integrated stream particles with progenitor orbit");

You can also obtain the stream points at the time of stripping by using integrate=False:

numpy.random.seed(1)
# Sample with integration and stripping times
orbs_not_integrated, dt = spdf_both.sample(n=500, returndt=True, integrate=False)

# Integrate the progenitor backwards for context
ts_prog = numpy.linspace(0.0, -4.5, 3001) * u.Gyr
prog.integrate(ts_prog, lp)

# Plot the integrated stream with progenitor orbit
prog.plot(d1="R", d2="z", color="k", lw=1, xrange=[0.0, 32.0], yrange=[-16.0, 16.0])
plt.scatter(
    orbs_not_integrated.R(), orbs_not_integrated.z(), c=dt, marker="o", s=20, alpha=0.5
)
plt.xlabel(r"$R$")
plt.ylabel(r"$z$")
plt.title("Stream particles at the time of stripping with progenitor orbit");

Different potential inputs

When initializing a particle-spray DF, you can also specify a different potential for computing the tidal radius and velocity distribution of the tidal debris, which can be useful when the overall potential contains pieces that are irrelevant for computing the tidal radius and that don’t allow the tidal radius to be computed (using the rtpot= option). If you want to generate a stream around a moving object, for example, a stream created within a satellite galaxy of the Milky Way, you can specify an orbit for the center of the satellite (center=) and the stream will be generated around this center rather than around the center of the total potential (this was used in Qian et al. 2022); the center orbit is integrated in centerpot, which can also differ from the potential that stream stars are integrated in (e.g., the stream stars may feel the potential from the satellite itself and/or the satellite could be experiencing dynamical friction which the stream stars do not feel).

Smooth stream track from spray samples

After drawing a set of stream particles with .sample(), we can construct a smooth phase-space track through them using .streamTrack(). The track is parameterized by an affine curve parameter tp along the progenitor’s orbit: for the leading arm tp runs from 0 (progenitor today) up to a positive tp_hi; for the trailing arm tp runs from a negative tp_lo up to 0. The exact range is auto-fit from the particle distribution and is much smaller than tdisrupt (typically a small fraction of an orbital period).

At each tp the track gives the mean galactocentric phase-space position of the corresponding stream particles, plus a 6×6 covariance that captures the scatter around that mean. The returned object is a StreamTrack — or a StreamTrackPair with .leading / .trailing attributes when tail='both' — that exposes the same coordinate accessors as Orbit: x, y, z, vx, vy, vz, R, vR, vT, phi, plus heliocentric (ra, dec, dist, ll, bb, pmra, pmdec, pmll, pmbb, vlos).

Accessors return NaN for tp values outside the track’s valid range (rather than silent cubic-spline extrapolation); use track.tp_grid() to see what range is in scope.

numpy.random.seed(0)
# Build the track for both arms
track = spdf_both.streamTrack(n=3000, tail="both")
# Also draw a particle cloud to overlay
orbs = spdf_both.sample(n=3000)
fig, ax = plt.subplots(figsize=(6, 6))
ax.scatter(orbs.x(), orbs.y(), s=1, alpha=0.2, color="0.6", label="samples")
track.plot(d1="x", d2="y", spread=1, color="C0")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend();

The shaded band shows a +/- 1 sigma envelope using the projected covariance at each point along the track. Pass spread=0 to suppress it.

# Heliocentric view: track overlaid on sample (ra, dec)
fig, ax = plt.subplots(figsize=(6, 5))
ax.scatter(orbs.ra(), orbs.dec(), s=1, alpha=0.2, color="0.6")
track.leading.plot(d1="ra", d2="dec", color="C0", label="leading")
track.trailing.plot(d1="ra", d2="dec", color="C3", label="trailing")
ax.set_xlabel("RA [deg]")
ax.set_ylabel("Dec [deg]")
ax.legend();

Controlling the smoothness and iteration

The smoothing keyword sets the spline smoothness; None (default) picks an automatic Generalized-Cross-Validation (GCV) value from the binned scatter via scipy.interpolate.make_smoothing_spline — for galpy’s normalized-internal units the GCV optimizer is run on a y-rescaled problem so it is well-behaved at the small absolute scales of the covariance entries. A float scales it; a dict keyed by 'x','y','z','vx','vy','vz' sets per-coordinate values; an array of length 6 (or 27 with order=2) supplies per-spline s values directly (e.g. from a previous track.smoothing_s).

smoothing_factor (default 1.0) is a second, simpler knob: it multiplies whatever s was selected (by GCV or by the user) by a constant, with >1 smoother and <1 rougher. The default GCV usually produces a good fit, but this dial lets you trade bias for variance if you have a particular reason to (e.g. a noisy short-stream sample where you’d rather oversmooth, or a long well-sampled stream where you want every wiggle of the predicted track resolved). See below for an example.

niter runs additional refinement iterations that reassign each particle to the closest point on the current track before refitting. For typical streams, 1-2 iterations suffice.

# Compare default fit with niter=2 (closest-point reassignment) on a single sample
spdf_both._tail = "leading"
tr_default = spdf_both.streamTrack(n=3000, tail="leading")
xv_dt = tr_default.particles  # save once, reuse below to skip the resample
tr_iter2 = spdf_both.streamTrack(particles=xv_dt, niter=2)
fig, ax = plt.subplots(figsize=(7, 5))
ax.scatter(orbs.x(), orbs.y(), s=1, alpha=0.2, color="0.6")
tr_default.plot(d1="x", d2="y", color="C0", label="default")
tr_iter2.plot(d1="x", d2="y", color="C2", label="niter=2")
ax.set_xlabel("x")
ax.set_ylabel("y")
ax.legend();

Pulling the 6×6 covariance with cov()

Beyond the mean track, track.cov(tp, basis=...) returns the 6×6 covariance of the particle distribution at the requested tp. Five output bases are supported: 'galcenrect' (default; native), 'galcencyl', 'sky' (equatorial), 'galsky' (Galactic), and 'customsky' (requires a custom_sky_transform). The non-native bases use analytical Jacobians — no finite differences — so the result is a self-consistent transform of the stored galactocentric covariance.

The diagonal entries can be used directly to draw ±σ bands; the cross-terms encode correlations (e.g. between pmra and pmdec):

tps = numpy.linspace(*track.leading.tp_grid()[[0, -1]], 5)
for tp in tps:
    C = track.leading.cov(tp, basis="sky")
    sigma_ra = numpy.sqrt(C[0, 0])
    sigma_dec = numpy.sqrt(C[1, 1])
    sigma_dist = numpy.sqrt(C[2, 2])
    print(
        f"tp={tp:6.3f}: σ(ra)={sigma_ra:.3f} deg, σ(dec)={sigma_dec:.3f} deg, σ(dist)={sigma_dist:.3f} kpc"
    )

tp= 0.000: σ(ra)=0.281 deg, σ(dec)=0.139 deg, σ(dist)=0.043 kpc
tp= 0.303: σ(ra)=0.291 deg, σ(dec)=0.287 deg, σ(dist)=0.044 kpc
tp= 0.606: σ(ra)=0.317 deg, σ(dec)=0.247 deg, σ(dist)=0.046 kpc
tp= 0.909: σ(ra)=0.286 deg, σ(dec)=0.132 deg, σ(dist)=0.049 kpc
tp= 1.212: σ(ra)=0.220 deg, σ(dec)=0.038 deg, σ(dist)=0.055 kpc

Stream-frame plots with custom_sky_transform

For visualizing thin streams, it’s often useful to rotate to a frame aligned with the stream itself — a (φ1, φ2) system in which the stream lies near φ2 = 0. Pass a 3×3 rotation matrix as custom_sky_transform= and the track exposes .phi1, .phi2, .pmphi1, .pmphi2 accessors and cov(basis='customsky'). The helper Orbit.align_to_orbit() builds such a matrix from the progenitor’s angular momentum:

Calling Orbit.align_to_orbit() also stashes the rotation matrix on the orbit itself, so afterwards prog.phi1(t), prog.phi2(t), prog.pmphi1(t), prog.pmphi2(t) (and prog.plot(d1='phi1', d2='phi2')) work without re-passing T=. Below: the same plot but also overlaying the progenitor’s orbit in the aligned frame — it traces phi2 = 0 by construction (the rotation maps the orbital plane onto the equator).

prog = prog()  # Need to reset because we integrated it above
T = prog.align_to_orbit()
track_cs = spdf_both.streamTrack(n=3000, tail="both", custom_sky_transform=T)

# To overlay the progenitor's *orbit* in the aligned frame we have to
# actually integrate it; `prog` is just initial conditions. The
# progenitor sits in the middle of the stream so we integrate both
# forward (to the leading-arm tip) and backward (to the trailing-arm
# tip) — the second integrate() call extends the first rather than
# overwriting it.
tp_lo = float(track_cs.trailing.tp_grid()[0])
tp_hi = float(track_cs.leading.tp_grid()[-1])
prog.integrate(numpy.linspace(0, tp_hi, 201), lp)
prog.integrate(numpy.linspace(0, tp_lo, 201), lp)

fig, ax = plt.subplots(figsize=(7, 3))
# `prog.align_to_orbit()` stashed the rotation matrix on `prog`, so
# `prog.phi1`/`prog.phi2` work without re-passing T=.
prog.plot(
    d1="phi1", d2="phi2", color="k", lw=1.0, ls="--", label="progenitor orbit", gcf=True
)
track_cs.leading.plot(d1="phi1", d2="phi2", color="C0", label="leading")
track_cs.trailing.plot(d1="phi1", d2="phi2", color="C3", label="trailing")
ax.scatter([180.0], [0.0], color="k", s=30, zorder=5, label="progenitor (today)")
ax.set_xlabel(r"$\phi_1$ [deg]")
ax.set_ylabel(r"$\phi_2$ [deg]")
ax.legend();

The custom_sky_transform is also exposed as a settable property on the StreamTrack (and on StreamTrackPair, which broadcasts the assignment to both arms), so you don’t have to rebuild the track to look at it in a different frame — or to drop the rotation altogether (set to None):

# Build the track without a rotation, then attach one after the fact —
# the phi1/phi2 accessors and cov(basis='customsky') light up immediately.
# (track_cs.particles is the both-arm sample from the previous cell; the
# leading-only xv_dt from earlier in the tutorial would be wrong here.)
track_bare = spdf_both.streamTrack(particles=track_cs.particles, tail="both")
print("before:", track_bare.custom_sky_transform)
track_bare.custom_sky_transform = T  # broadcasts to both arms of the pair
print("after :  shape", track_bare.custom_sky_transform.shape)

fig, ax = plt.subplots(figsize=(7, 2.5))
track_bare.leading.plot(d1="phi1", d2="phi2", color="C0", label="leading")
track_bare.trailing.plot(d1="phi1", d2="phi2", color="C3", label="trailing")
ax.set_xlabel(r"$\phi_1$ [deg]")
ax.set_ylabel(r"$\phi_2$ [deg]")
ax.legend();

before: None
after :  shape (3, 3)

Reusing samples for a smoothing sweep

The expensive part of streamTrack is the orbit-integrated sample (and, on first call, the GCV per-spline s selection). The cheap part is the spline fit. To explore smoothing variants without redoing the expensive work, save the particles and pass them back as particles= to subsequent calls — only the bin-and-smooth step re-runs.

There are three reuse handles depending on what object you have in hand:

• track.particles (a single-arm StreamTrack) — the xv array the fit saw. Pass back as spdf.streamTrack(particles=..., tail='leading') (or 'trailing').

• pair.particles (a StreamTrackPair from tail='both') — concatenates the two arms’ xv arrays in the leading-first order that streamTrack(tail='both', particles=...) expects.

• spdf.sample(tail='both', returndt=False, return_orbit=False, integrate=True) — pre-sample without building any track at all; returns the same concatenated xv array.

Two further knobs that take the cheap-refit path:

• track.smoothing_s is the per-spline effective s. Passing it back as smoothing= (with smoothing_factor=1) skips the GCV step entirely and re-fits at exactly that smoothness.

• smoothing_factor= is the simplest dial: rerun with the same particles= and a different multiplier.

Below: the same track_cs sample fit at smoothing_factor \(\in \{0.25, 1, 4\}\) in the aligned frame, both arms. The factor 0.25 track follows individual bin means and the factor 4 track sands out the corners — the default 1.0 is a balanced fit.

# Reuse the both-arms sample from track_cs (built in the previous cell) —
# pair.particles concatenates the two arms in the leading-first order that
# streamTrack(tail='both') expects, so re-fitting only re-runs the splines.
xv_dt = track_cs.particles
tr_cs_rough = spdf_both.streamTrack(
    particles=xv_dt, tail="both", custom_sky_transform=T, smoothing_factor=0.25
)
tr_cs_smooth = spdf_both.streamTrack(
    particles=xv_dt, tail="both", custom_sky_transform=T, smoothing_factor=4.0
)
fig, ax = plt.subplots(figsize=(8, 3.5))
for arm in ("leading", "trailing"):
    getattr(tr_cs_rough, arm).plot(
        d1="phi1",
        d2="phi2",
        color="C3",
        alpha=0.7,
        label="factor=0.25 (rougher)" if arm == "leading" else None,
    )
    getattr(track_cs, arm).plot(
        d1="phi1",
        d2="phi2",
        color="C0",
        label="factor=1.0 (default)" if arm == "leading" else None,
    )
    getattr(tr_cs_smooth, arm).plot(
        d1="phi1",
        d2="phi2",
        color="C2",
        label="factor=4.0 (smoother)" if arm == "leading" else None,
    )
ax.set_xlabel(r"$\phi_1$ [deg]")
ax.set_ylabel(r"$\phi_2$ [deg]")
ax.legend(loc="best", fontsize=9);

Non-uniform stripping times

By default, streamspraydf draws stripping times uniformly between -tdisrupt and 0. Real tidal streams, however, are stripped preferentially near pericenter passages of the progenitor’s orbit, where the tidal field is strongest. The stripping_pdf= keyword lets you supply an arbitrary probability density over the progenitor time axis t \in [-t_\mathrm{disrupt}, 0] (t=0 is the present day, t<0 the past). Internally the PDF is converted to a CDF and inverted; samples are drawn via inverse-CDF transform.

The PDF need not be normalized. It may accept and return floats, numpy arrays, or astropy Quantity (input a time, output a 1/time); the detection is the same as for AnyAxisymmetricRazorThinDiskPotential.

Example A: a custom Gaussian stripping window

Suppose stripping was concentrated around 2.5 Gyr ago, with a width of 0.3 Gyr. We pass a Gaussian PDF in Quantity form and sample the stream.

from galpy.util import conversion

ro_, vo_ = 8.0, 220.0


def gauss_pdf(t):
    return numpy.exp(-0.5 * ((t - (-2.5 * u.Gyr)) / (0.3 * u.Gyr)) ** 2) / u.Gyr


spdf_gauss = chen24spraydf(
    progenitor_mass=2 * 10**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    stripping_pdf=gauss_pdf,
    tail="both",
)
numpy.random.seed(0)
stream_gauss, dt_gauss = spdf_gauss.sample(n=600, returndt=True)
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.hist(dt_gauss * conversion.time_in_Gyr(vo_, ro_), bins=40, color="steelblue")
plt.xlabel(r"$\Delta t$ since stripping [Gyr]")
plt.ylabel("count")
plt.title("Gaussian stripping_pdf")
plt.subplot(1, 2, 2)
plt.plot(stream_gauss.x(), stream_gauss.y(), "k.", ms=1)
plt.plot(prog.x(), prog.y(), "ro", ms=6)
plt.xlabel(r"$x$ [kpc]")
plt.ylabel(r"$y$ [kpc]")
plt.tight_layout()
plt.show()

Example B: pericenter-driven stripping

A common physical model is to bunch stripping into Gaussians centered on each pericenter passage. galpy.df.pericenter_stripping_pdf builds such a PDF from the progenitor and potential, given a width sigma. The returned PDF can be passed straight into stripping_pdf=, and it exposes a .pericenter_times attribute listing the pericenter times in internal units.

from galpy.df import pericenter_stripping_pdf
from galpy.util import conversion

ro_, vo_ = 8.0, 220.0
peri_pdf = pericenter_stripping_pdf(
    prog, lp, 4.5 * u.Gyr, sigma=50 * u.Myr, ro=ro_, vo=vo_
)
spdf_peri = chen24spraydf(
    progenitor_mass=2 * 10**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    stripping_pdf=peri_pdf,
    tail="both",
)
spdf_uniform = chen24spraydf(
    progenitor_mass=2 * 10**4.0 * u.Msun,
    progenitor=prog,
    pot=lp,
    tdisrupt=4.5 * u.Gyr,
    tail="both",
)
numpy.random.seed(1)
stream_peri, dt_peri = spdf_peri.sample(n=600, returndt=True)
numpy.random.seed(1)
stream_uniform, dt_uniform = spdf_uniform.sample(n=600, returndt=True)
peri_gyr = peri_pdf.pericenter_times * conversion.time_in_Gyr(vo_, ro_)
plt.figure(figsize=(10, 4))
plt.subplot(1, 2, 1)
plt.hist(
    dt_uniform * conversion.time_in_Gyr(vo_, ro_),
    bins=40,
    color="0.6",
    alpha=0.6,
    label="uniform",
)
plt.hist(
    dt_peri * conversion.time_in_Gyr(vo_, ro_),
    bins=40,
    color="C3",
    alpha=0.6,
    label="pericenter",
)
for tp in peri_gyr:
    plt.axvline(-tp, color="k", lw=0.5, alpha=0.6)
plt.xlabel(r"$\Delta t$ since stripping [Gyr]")
plt.ylabel("count")
plt.legend()
plt.subplot(1, 2, 2)
plt.plot(stream_uniform.x(), stream_uniform.y(), "k.", ms=1, alpha=0.4, label="uniform")
plt.plot(stream_peri.x(), stream_peri.y(), "C3.", ms=1, alpha=0.6, label="pericenter")
plt.plot(prog.x(), prog.y(), "yo", ms=6)
plt.xlabel(r"$x$ [kpc]")
plt.ylabel(r"$y$ [kpc]")
plt.legend()
plt.tight_layout()
plt.show()