nbins = 33
sp = ds.sphere(center, radius)
nz_center = sp['spherical_r'].nonzero()
dr_min = sp['spherical_r'][nz_center].min()
extrema = dict(x=(0.99*dr_min, outer_factor * radii[_cid]),
               y=(-0.001*np.pi, 1.001*np.pi),
               z=(-1.001*np.pi, 1.001*np.pi))
p3d = yt.Profile3D(sp, 'spherical_r', nbins, extrema['x'][0], extrema['x'][1], True,
                   'spherical_theta', nbins, extrema['y'][0], extrema['y'][1], False,
                   'spherical_phi', nbins, extrema['z'][0], extrema['z'][1],False,
                   weight_field='cell_mass')

print("Binning data into spherical grid")
p3d.add_fields(['density', 'velocity_x', 'velocity_y', 'velocity_z', 'ones'])

# Calculate the bin numbers for each data point

logr = np.log10(sp['spherical_r'].to('cm'))
r0 = np.log10(extrema['x'][0].to('cm'))
r1 = np.log10(extrema['x'][1].to('cm'))
dr = (r1 - r0) / (nbins)
ir = ((logr - r0) / dr).astype('int')
itheta = ((sp['spherical_theta'] - extrema['y'][0]) / (p3d.y[1] - p3d.y[0])).astype('int')
iphi = ((sp['spherical_phi'] - extrema['z'][0]) / (p3d.z[1] - p3d.z[0])).astype('int')
ir = np.minimum(np.maximum(ir, 0), nbins-1)

# Calculate rms velocity in each spherical grid cell

print("Calculating fluid velocities (ncell=%d)" % (ncells))
ncells = sp['ones'].size
deltav = ds.arr(np.zeros((ncells,3)), 'cm/s')
for idim, dim in enumerate('xyz'):
    f = 'velocity_%c' % (dim)
    deltav[:,idim] = sp[f] - p3d[f][ir,itheta,iphi]
vrms = (deltav**2).sum(1)**0.5