1import matplotlib.pyplot as plt
2
3import yt
4
5ds = yt.load("GasSloshing/sloshing_nomag2_hdf5_plt_cnt_0150")
6
7# Get the first sphere
8sp0 = ds.sphere(ds.domain_center, (500.0, "kpc"))
9
10# Compute the bulk velocity from the cells in this sphere
11bulk_vel = sp0.quantities.bulk_velocity()
12
13
14# Get the second sphere
15sp1 = ds.sphere(ds.domain_center, (500.0, "kpc"))
16
17# Set the bulk velocity field parameter
18sp1.set_field_parameter("bulk_velocity", bulk_vel)
19
20# Radial profile without correction
21
22rp0 = yt.create_profile(
23    sp0,
24    ("index", "radius"),
25    ("gas", "radial_velocity"),
26    units={("index", "radius"): "kpc"},
27    logs={("index", "radius"): False},
28)
29
30# Radial profile with correction for bulk velocity
31
32rp1 = yt.create_profile(
33    sp1,
34    ("index", "radius"),
35    ("gas", "radial_velocity"),
36    units={("index", "radius"): "kpc"},
37    logs={("index", "radius"): False},
38)
39
40# Make a plot using matplotlib
41
42fig = plt.figure()
43ax = fig.add_subplot(111)
44
45ax.plot(
46    rp0.x.value,
47    rp0[("gas", "radial_velocity")].in_units("km/s").value,
48    rp1.x.value,
49    rp1[("gas", "radial_velocity")].in_units("km/s").value,
50)
51
52ax.set_xlabel(r"$\mathrm{r\ (kpc)}$")
53ax.set_ylabel(r"$\mathrm{v_r\ (km/s)}$")
54ax.legend(["Without Correction", "With Correction"])
55
56fig.savefig(f"{ds}_profiles.png")
57