SHARE
TWEET

Untitled

a guest Apr 25th, 2019 65 Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
  1. nbins = 33
  2. sp = ds.sphere(center, radius)
  3. nz_center = sp['spherical_r'].nonzero()
  4. dr_min = sp['spherical_r'][nz_center].min()
  5. extrema = dict(x=(0.99*dr_min, outer_factor * radii[_cid]),
  6.                y=(-0.001*np.pi, 1.001*np.pi),
  7.                z=(-1.001*np.pi, 1.001*np.pi))
  8. p3d = yt.Profile3D(sp, 'spherical_r', nbins, extrema['x'][0], extrema['x'][1], True,
  9.                    'spherical_theta', nbins, extrema['y'][0], extrema['y'][1], False,
  10.                    'spherical_phi', nbins, extrema['z'][0], extrema['z'][1],False,
  11.                    weight_field='cell_mass')
  12.  
  13. print("Binning data into spherical grid")
  14. p3d.add_fields(['density', 'velocity_x', 'velocity_y', 'velocity_z', 'ones'])
  15.  
  16. # Calculate the bin numbers for each data point
  17.  
  18. logr = np.log10(sp['spherical_r'].to('cm'))
  19. r0 = np.log10(extrema['x'][0].to('cm'))
  20. r1 = np.log10(extrema['x'][1].to('cm'))
  21. dr = (r1 - r0) / (nbins)
  22. ir = ((logr - r0) / dr).astype('int')
  23. itheta = ((sp['spherical_theta'] - extrema['y'][0]) / (p3d.y[1] - p3d.y[0])).astype('int')
  24. iphi = ((sp['spherical_phi'] - extrema['z'][0]) / (p3d.z[1] - p3d.z[0])).astype('int')
  25. ir = np.minimum(np.maximum(ir, 0), nbins-1)
  26.  
  27. # Calculate rms velocity in each spherical grid cell
  28.  
  29. print("Calculating fluid velocities (ncell=%d)" % (ncells))
  30. ncells = sp['ones'].size
  31. deltav = ds.arr(np.zeros((ncells,3)), 'cm/s')
  32. for idim, dim in enumerate('xyz'):
  33.     f = 'velocity_%c' % (dim)
  34.     deltav[:,idim] = sp[f] - p3d[f][ir,itheta,iphi]
  35. vrms = (deltav**2).sum(1)**0.5
RAW Paste Data
We use cookies for various purposes including analytics. By continuing to use Pastebin, you agree to our use of cookies as described in the Cookies Policy. OK, I Understand
 
Top