%matplotlib inline
from matplotlib import pyplot as plt
import yt
import numpy as np
ds = yt.load('IsolatedGalaxy/galaxy0030/galaxy0030')
prj_plot = yt.ProjectionPlot(ds, 2, 'density', width=(20, 'kpc'))
prj_plot.set_figure_size(4)
# Construct a cylindrical radius image
surface_density_image = prj_plot.frb['density']
shape = surface_density_image.shape
x = np.linspace(-10, 10, shape[0])
x = np.tile(x, shape[1])
x.shape = shape
y = np.flipud(x.transpose())
radius_image = (np.sqrt(x**2 + y**2)*yt.units.kpc.in_cgs()).d
plt.imshow(radius_image)
<matplotlib.image.AxesImage at 0x1148c9a90>
import numpy as np
radii_bins = (np.linspace(0, 15, 100)*yt.units.kpc.in_cgs()).d
# here's the tricky bit - digitize counts the number of cells in each radius bin
digitized_bins = np.digitize(radius_image.flat, radii_bins)
# a weighted bincount over an unweighted count is the average surface density in each bin
profile = np.bincount(digitized_bins, weights=surface_density_image.flat, minlength=99) / np.bincount(digitized_bins, minlength=99)
-c:9: RuntimeWarning: invalid value encountered in divide
radii = (radii_bins[1:]+radii_bins[:1])/2.0
plt.plot(radii, profile)
plt.xlabel('Radius (cm)')
plt.ylabel('Surface Density ($g/cm^2$)')
<matplotlib.text.Text at 0x1148f6210>