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