<jupyter_start><jupyter_text>Parameters:	*data (dict)* – This is a dict of numpy arrays or (numpy array, unit name) tuples, where the keys are the field names. Particles positions must be named “particle_position_x”, “particle_position_y”, and “particle_position_z”.*length_unit (float)* – Conversion factor from simulation length units to centimetersmass_unit (float) – Conversion factor from simulation mass units to gramstime_unit (float) – Conversion factor from simulation time units to secondsvelocity_unit (float) – Conversion factor from simulation velocity units to cm/smagnetic_unit (float) – Conversion factor from simulation magnetic units to gaussbbox (array_like (xdim:zdim, LE:RE), optional) – Size of computational domain in units of the length_unitsim_time (float, optional) – The simulation time in secondsperiodicity (tuple of booleans) – Determines whether the data will be treated as periodic along each axisn_ref (int) – The number of particles that result in refining an oct used for indexing the particles.<jupyter_code>from yt.units import kiloparsec, Msun 

length_unit=3.08567758149137e21 #kpc to cm  
mass_unit=1.98847e33 #solar mass to g 
time_unit=3.15576e16 #Gyr to s, 
velocity_unit=1e5 #km to cm   
unit_system="galactic" #'cgs'
sim_time = part.snapshot['time']*time_unit
field = 'io'
data = {(field, 'particle_position_x'): part['gas']['position'][:,0],
        (field, 'particle_position_y'): part['gas']['position'][:,1],
        (field, 'particle_position_z'): part['gas']['position'][:,2],
        (field, 'particle_velocity_x'): part['gas']['velocity'][:,0],
        (field, 'particle_velocity_y'): part['gas']['velocity'][:,1],
        (field, 'particle_velocity_z'): part['gas']['velocity'][:,2],
        (field, 'density'):    part['gas']['density'],
        (field, 'hydrogen.neutral.fraction'): part['gas']['hydrogen.neutral.fraction'],
        (field, 'massfraction'): part['gas']['massfraction'],
        (field, 'metallicity.total'): part['gas']['metallicity.total'],
        (field, 'smoothing_length'): part['gas']['smooth.length']*u.kpc,
        (field, 'particle_mass'): part['gas']['mass']
        }
half_box = part.info['box.length']/2.*length_unit
bbox = np.array([[-half_box, half_box], [-half_box, half_box], [-half_box, half_box]])

ds = yt.load_particles(data, length_unit=length_unit, mass_unit=mass_unit, velocity_unit=velocity_unit, 
                       time_unit=time_unit, n_ref=64, bbox=bbox, unit_system=unit_system, sim_time=sim_time)
ad = ds.all_data()
xarray = np.linspace(-1, 1, 16)
yarray = np.linspace(-1, 1, 16)
zarray = np.linspace(-1, 1, 16)
ad[(field, 'particle_position_x')]
from yt.units import kpc
left_edge = [0, 0, 0]*kpc
right_edge = [.99, .99, .99]*kpc
ag = ds.arbitrary_grid(left_edge, right_edge, dims=[128]*3)
ag_density = ag[(field, 'density')]
arbitrary_grid = ds.arbitrary_grid([0.0, 0.0, 0.0], [0.99, 0.99, 0.99] ,
                              dims=[128, 128, 128])

print('arbitrary grid shape: ', ag.shape)
print('density grid shape: ', ag_density.shape)

plt.hist(ad[('gas', 'particle_velocity_x')].in_units('km/s')/yt.YTArray(part['gas']['velocity'][:,0], input_units='km/s') - 1.)
dx = ad[('gas', 'particle_position_x')].in_units('kpc') - yt.YTArray(part['gas']['position'][:,0], input_units='kpc')
dy = ad[('gas', 'particle_position_y')].in_units('kpc') - yt.YTArray(part['gas']['position'][:,1], input_units='kpc')
dz = ad[('gas', 'particle_position_z')].in_units('kpc') - yt.YTArray(part['gas']['position'][:,2], input_units='kpc')
xbins = np.linspace(-1, 1, 100)
ybins = np.linspace(-1, 1, 100)
zbins = np.linspace(-1, 1, 100)

fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for xx, yy, xl, yl, ax in zip([dx, dx, dy],[dy, dz, dz], ['dx', 'dx', 'dy'], ['dy', 'dy', 'dz'], axes):
    H, xe, ye = np.histogram2d(xx, yy, bins=(xbins, ybins))
    ax.pcolormesh(xe, ye, H.T,
      norm=mpl.colors.LogNorm(), cmap=plt.get_cmap('Blues'))
    ax.set_xlabel(xl)
    ax.set_ylabel(yl)

dr = np.sqrt(dx**2 + dy**2 + dz**2)
r = np.sqrt(np.sum(part['gas']['position']**2., axis=1))
plt.plot(r[np.argsort(r)], dr[np.argsort(r)])
plt.xlabel('distance from center of simulation [kpc]')
plt.ylabel('dr [kpc]')
plt.title('Distance between "true" positions and yt ds positions')
plt.savefig('particlePositionDifferences.png')
from scipy.stats import binned_statistic_dd 

xbins = np.linspace(-10, 10, 10)
ybins = np.linspace(-10, 10, 10)
zbins = np.linspace(-10, 10, 10)

sample = part['gas']['position']
values = dr
H, be, bn = binned_statistic_dd(sample, values, statistic='mean', bins=(xbins, ybins, zbins))

r = np.sqrt(np.sum(part['gas']['position']**2., axis=1))
H.shape
from mpl_toolkits.axes_grid1 import make_axes_locatable

def colorbar(mappable, ax):
    """ sol'n from https://joseph-long.com/writing/colorbars/"""
    #ax = mappable.axes
    fig = ax.figure
    divider = make_axes_locatable(ax)
    cax = divider.append_axes("right", size="5%", pad=0.05)
    return fig.colorbar(mappable, cax=cax)

fig, axes = plt.subplots(1, 5, figsize=(15, 2.5))
for index, ax in zip([0, 3, 4, 5, 8], axes):
    im = H[:,:,index].T
    im[np.isnan(im)] = 0.0
    im = ax.pcolormesh(be[0], be[1], im, vmin=1e-4, vmax=3e-3, norm=mpl.colors.LogNorm(), cmap=plt.get_cmap('Blues'))
    colorbar(im, ax)
    ax.set_xlabel('x')
    ax.set_ylabel('y')
    ax.axis('equal')
    ax.set_title('z = {0:.2f}'.format(0.5*(zbins[index+1] + zbins[index])))
plt.tight_layout()
fig.suptitle('Average distance between "true" particle positions and yt ds positions', y=1.05)
fig.savefig('differenceInPositions.pdf', rasterized=True)
fig, axes = plt.subplots(1, 3, figsize=(15, 5))

for xx, yy, xl, yl, ax in zip([dx, dx, dy],[dy, dz, dz], ['dx', 'dx', 'dy'], ['dy', 'dy', 'dz'], axes):
    H, be, bn = 
    H, xe, ye = np.histogram2d(xx, yy, bins=(xbins, ybins))
    ax.pcolormesh(xe, ye, H.T,
      norm=mpl.colors.LogNorm(), cmap=plt.get_cmap('Blues'))
    ax.set_xlabel(xl)
    ax.set_ylabel(yl)

xe
ad[('gas')]
ad[('gas', 'smoothing_length')].in_units('kpc')
part['gas']['smooth.length']
part
filename = '/mnt/ceph/users/firesims/ananke/Latte/' + simname + '/snapdir_600/snapshot_600.0.hdf5'
#filename = simname + '/snapdir_600/snapshot_600.0.hdf5'
ds = yt.load(filename)
ad = ds.all_data()
rotate_gizmo.read.Read.read_galaxy_center_coordinates(ad, directory='.')
ad.principal_axes_vectors
rotate_gizmo.read.Read.read_lsr_coordinates(ad, directory='.', lsr_index=2)
ad.set_field_parameter('center', ad.center_position*u.kpc)
ad.set_field_parameter('bulk_velocity', ad.center_velocity*u.km/u.s)
ad.set_field_parameter('rotation_vectors', ad.principal_axes_vectors)
ad.principal_axes_vectors[0]
ad.get_field_parameter("rotation_vectors")
ad.ds.arr(np.eye(3))[0]
part = gizmo_read.read.Read.read_snapshot(species='all', 
                                          properties= 'all', 
                                          directory=directory)
for key in part.keys():
    part[key]['position'] = part[key]['position'] - part.lsr_position
    part[key]['velocity'] = part[key]['velocity'] - part.lsr_velocity

print(part.center_position, part.center_velocity, part.principal_axes_vectors)
ad.set_field_parameter?
ad[('all', 'rotCoordinates')]
ds.unit_registry.keys()
ad[('PartType0', 'Coordinates')]
part['gas']['position']
gizmo_read.read.Read.read_galaxy_center_coordinates?
ds.add_field?
positions = ad['all', 'Coordinates']  - part.lsr_position
velocities = ad['all', 'Velocities'] - part.lsr_velocity
#density ==> rho, smooth.length ==> smooth, mass ==> mass 
sim.gas['pos'] = pynbody.array.SimArray(positions)
sim.gas['pos'].units = 'kpc'

sim.gas['vel'] = pynbody.array.SimArray(velocities)
sim.gas['vel'].units = 'km s^-1'

sim.gas['eps'] = sim.gas['smooth']
sim.gas['eps'].units = 'kpc'


sim.gas['massfraction.all'] = part['gas']['massfraction'][:,0]
sim.gas['massfraction.all'].units = '1'
sim.gas['massfraction.he']  = part['gas']['massfraction'][:,1]
sim.gas['massfraction.he'].units = '1'

sim.gas['massfraction.all.relativesolar'] = 10.**part['gas']['metallicity.total']
sim.gas['massfraction.all.relativesolar'].units = '1'

for key in ['hydrogen.neutral.fraction']:
    sim.gas[key] = pynbody.array.SimArray(part['gas'][key])
    sim.gas[key].units = '1'

sim.gas['smooth']
solarMetallicity = 0.2
fractionHydrogen = 1 - sim.gas['massfraction.he'] - sim.gas['massfraction.all']
neutralHydrogenMass = sim.gas['mass']*(fractionHydrogen)*sim.gas['hydrogen.neutral.fraction']
metalWeightedNeutralHydrogenMass = neutralHydrogenMass*sim.gas['massfraction.all.relativesolar']
metalWeightedNeutralHydrogenMass
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
ax[0].hist(np.log10(fractionHydrogen*sim.gas['hydrogen.neutral.fraction']*sim.gas['massfraction.all.relativesolar']), 
           range = [-10, 1], bins=100, log=True, histtype='step')
ax[0].set_title('Metal Weighted Neutral Hydorgen Fraction')
ax[0].set_xlabel('Metal Weighted Neutral Hydorgen Fraction')

ax[1].hist(np.log10(neutralHydrogenMass), bins=100, log=True, range = [-5, 1], histtype='step')
ax[1].set_title('Neutral Hydrogen Mass')
ax[1].set_xlabel('Neutral Hydrogen Mass [Msol]')

ax[2].hist(np.log10(metalWeightedNeutralHydrogenMass), bins=100, log=True, range = [-5, 1], histtype='step')
ax[2].set_title('Metal Weighted Neutral Hydrogen Mass')
ax[2].set_xlabel('Metal Weighted Neutral Hydrogen Mass [Msol]')

sim.gas['metalWeightedNeutralHydrogenMass'] = pynbody.array.SimArray(metalWeightedNeutralHydrogenMass.in_units('Msol'))
sim.gas['metalWeightedNeutralHydrogenMass'].units = 'Msol'

sim.gas['neutralHydrogenMass'] = pynbody.array.SimArray(neutralHydrogenMass.in_units('Msol'))
sim.gas['neutralHydrogenMass'].units = 'Msol'

sim.gas['metalWeightedNeutralHydrogenNumber'] = pynbody.array.SimArray(metalWeightedNeutralHydrogenMass.in_units('kg'), 
                                                                       dtype='float64')/m_p.to('kg')
sim.gas['metalWeightedNeutralHydrogenNumber'].units = '1'
sim.gas['metalWeightedNeutralHydrogenNumber']
sim.gas['metalWeightedNeutralHydrogenMass']
sim.gas.keys()
for k in sim.gas.keys(): print(k, sim.gas[k].dtype)
for key in ['smooth', 'mass', 'metalWeightedNeutralHydrogenMass', 'eps', 'rho']:
    arr = np.asanyarray(sim.gas[key], dtype=sim.gas['pos'].dtype)
    del sim.gas[key]
    sim.gas[key]= arr
x1 = "-750 pc"
x2 = "750 pc"
y1 = "-750 pc"
y2 = "750 pc"
z1 = "-750 pc"
z2 = "750 pc"
box = pynbody.filt.Cuboid(x1, y1=y1, z1=z1, x2=x2, y2=y2, z2=z2)

im = pynbody.sph.to_3d_grid(sim[box].gas, qty='rho', nx=15, ny=15, nz=15, threaded=False, approximate_fast=False)

scale = '0.5 kpc'
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
pynbody.plot.image(sim.gas, qty='rho', width = '100 kpc', 
                   subplot=ax[0], show_cbar=False, title='Gas Density')
pynbody.plot.image(sim.gas, qty='metalWeightedNeutralHydrogenMass', width = '100 kpc', 
                   subplot=ax[1], show_cbar=False, title='Metal Weighted')
plt.tight_layout()
plt.savefig('m12flsr2faceon.pdf', rasterized=True)
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
pynbody.plot.image(sim.gas, qty='rho', width = scale, 
                   subplot=ax[0], show_cbar=False, title='gas density')
pynbody.plot.image(sim.gas, qty='metalWeightedNeutralHydrogenMass', width = scale, 
                   subplot=ax[1], show_cbar=False, title='Metal Weighted')
plt.tight_layout()
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
sim.gas.rotate_x(90)
pynbody.plot.image(sim.gas, qty='rho', width = '100 kpc', 
                   subplot=ax[0], show_cbar=False, title='Gas Density')
pynbody.plot.image(sim.gas, qty='metalWeightedNeutralHydrogenMass', width = '100 kpc', 
                   subplot=ax[1], show_cbar=False, title='Metal Weighted')
sim.gas.rotate_x(-90)
plt.tight_layout()
plt.savefig('m12flsr2sideone.pdf', rasterized=True)
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
sim.gas.rotate_x(90)
pynbody.plot.image(sim.gas, qty='rho', width = scale, 
                   subplot=ax[0], show_cbar=False, title='gas density')
pynbody.plot.image(sim.gas, qty='metalWeightedNeutralHydrogenMass', width = scale, 
                   subplot=ax[1], show_cbar=False, title='Metal Weighted')
sim.gas.rotate_x(-90)
plt.tight_layout()
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
sim.gas.rotate_y(90)
pynbody.plot.image(sim.gas, qty='rho', width = '100 kpc', 
                   subplot=ax[0], show_cbar=False, title='Gas Density')
pynbody.plot.image(sim.gas, qty='metalWeightedNeutralHydrogenMass', width = '100 kpc', 
                   subplot=ax[1], show_cbar=False, title='Metal Weighted')
sim.gas.rotate_y(-90)
plt.tight_layout()
plt.savefig('m12flsr2sideon2.pdf', rasterized=True)
fig, ax = plt.subplots(1, 2, figsize=(10, 5))
sim.gas.rotate_y(90)
pynbody.plot.image(sim.gas, qty='rho', width = scale, 
                   subplot=ax[0], show_cbar=False, title='Neutral Hydrogen')
pynbody.plot.image(sim.gas, qty='metalWeightedNeutralHydrogenMass', width = scale, 
                   subplot=ax[1], show_cbar=False, title='Metal Weighted')
sim.gas.rotate_y(-90)
plt.tight_layout()
np.linspace(-400, 400.0001, 5)
x1 = "-10 kpc"
x2 = "10 kpc"
y1 = "-10 kpc"
y2 = "10 kpc"
zup = np.linspace(-1000, 1000.0001, 10)[1:]
zdown = np.linspace(-1000, 1000.0001, 10)[:-1]

for i, (zu, zd) in enumerate(zip(zup, zdown)):
        
    box = pynbody.filt.Cuboid(x1, y1=y1, z1=str(zd) + ' pc', x2=x2, y2=y2, z2=str(zu) + ' pc')
    fig, ax = plt.subplots(1, 2, figsize=(10, 5))
    im = pynbody.plot.image(sim[box].gas, qty='rho', width = scale, 
                   subplot=ax[0], show_cbar=False, title='rho z {0:.1f} - {1:.1f} pc'.format(zu, zd))

    im = pynbody.plot.image(sim[box].gas, qty='metalWeightedNeutralHydrogenMass', width =  scale, 
                   subplot=ax[1], show_cbar=True, title='Metal Weighted', vmax=1e4, vmin=1e2)

    plt.tight_layout()
    fig.savefig('anankedust{0}.pdf'.format(i), rasterized=True)
radius = "2 kpc" # or might be a float if you already know your units
centre = (0,0,0) # Take the origin for now
sphere = pynbody.filt.Sphere(radius, centre)
sim[sphere].gas
pynbody.sph.to_3d_grid?
pynbody.__version__
x1 = "-750 pc"
x2 = "750 pc"
y1 = "-750 pc"
y2 = "750 pc"
z1 = "-750 pc"
z2 = "750 pc"
box = pynbody.filt.Cuboid(x1, y1=y1, z1=z1, x2=x2, y2=y2, z2=z2)

im = pynbody.sph.to_3d_grid(sim[box].gas, qty='mass', nx=15, ny=15, nz=15, threaded=False, approximate_fast=False)

H0 = sim.gas['smooth']
m0 = sim.gas['mass']
z0 = part['gas']['metallicity.total']
#z0 = sim.gas['metallicity.total']
yHe = part['gas']['metallicity.he']
#yHe = sim.gas['metallicity.he']
fH = 1 - yHe - z0
fneutral = part['gas']['hydrogen.neutral.fraction']
#fneutral = sim.gas['hydrogen.neutral.fraction']
fH_0 = fH*fneutral
hubble = sim.properties['h']
Hmin = 0.25*0.0131968*hubble
H0 = np.sqrt(H0*H0 + Hmin*Hmin)
z0 = np.sqrt(z0*z0 + 1e-5*1e-5)
wt = fH_0*m0*(z0/0.2)/(H0 * H0 * H0) #*hubble took out, but need to check
fH_0*z0/0.2
numberH = pynbody.array.SimArray(m0.in_units('kg'), dtype='float64')*fH_0*z0/0.2/m_p.value
numberH.sim = sim # this links the array to the simulation

numberH.units = '1'
numberH
pynbody.sph.to_3d_grid(sim.gas, qty=)
radius = np.sqrt(np.sum(part['gas']['position']**2, axis=1))
smoothlength = (part['gas']['mass']/part['gas']['density'])**(1/3)
plt.hist(np.log10(smoothlength[radius < 2]), log=True, bins=100)
plt.xlabel('smoothing lengths [kpc]')
plt.title('smoothing lengths inside 2 kpc')
xbins = np.linspace(0, 100, 100)
ybins = np.linspace(0, 2.5, 100)
H, xe, ye = np.histogram2d(radius, smoothlength, bins=[xbins,ybins])
im = plt.pcolormesh(xe, ye, H.T,
  norm=mpl.colors.LogNorm(), cmap=plt.get_cmap('Blues'))#vmax=75))#, vmax=20) #, vmin=-100, vmax=100)
im.set_rasterized(True)
plt.xlabel('radius [kpc]')
plt.ylabel('smoothing length [kpc]')
plt.colorbar()
plt.savefig('smoothinglengths.pdf')

xbins = np.linspace(-1, 2, 100)
ybins = np.linspace(-3, 0.5, 100)
H, xe, ye = np.histogram2d(np.log10(radius), np.log10(smoothlength), bins=[xbins,ybins])
im = plt.pcolormesh(xe, ye, H.T,
  norm=mpl.colors.LogNorm(), cmap=plt.get_cmap('Blues'))#vmax=75))#, vmax=20) #, vmin=-100, vmax=100)
im.set_rasterized(True)
plt.xlabel('log10 radius [kpc]')
plt.ylabel('log10 smoothing length [kpc]')
plt.colorbar()
plt.savefig('smoothinglengths_logspace.pdf')
part['gas'].keys()
part.info['hubble']

H0 = smoothlength
m0 = part['gas']['mass']
z0 = part['gas']['metallicity.total']
yHe = part['gas']['metallicity.he']
fH = 1 - yHe - z0
fneutral = part['gas']['hydrogen.neutral.fraction']
fH_0 = fH*fneutral
hubble = part.info['hubble']
Hmin = 0.25*0.0131968*hubble
H0 = np.sqrt(H0*H0 + Hmin*Hmin)
z0 = np.sqrt(z0*z0 + 1e-5*1e-5)
wt = fH_0*m0*(z0/0.2)*hubble*2.08854/(H0 * H0 * H0)*6.7702543e-32
wt
plt.hist(smoothlength*2. - sim.gas['smooth'], bins=100, log=True)
plt.title('difference between smoothing lengths')
plt.xlabel('delta smoothing length [kpc]')
smoothlength*2.
sim.gas['smooth']
sim.properties
part.info
pynbody.plot.image(sim.gas, qty='rho', threaded=False)
pynbody.plot.image?
pynbody.plot.sideon_image(sim.gas)
sim.all_keys()
sim['x']<jupyter_output><empty_output>
{'fraction_non_alphanumeric': 0.18664714494875548, 'fraction_numerical': 0.03976573938506588, 'mean_word_length': 4.951899616591147, 'pattern_counts': {'":': 0, '<': 6, '<?xml version=': 0, '>': 8, 'https://': 1, 'lorem ipsum': 0, 'www.': 0, 'xml': 0}, 'pii_count': 0, 'substrings_counts': 0, 'word_list_counts': {'cursed_substrings.json': 7, 'profanity_word_list.json': 0, 'sexual_word_list.json': 0, 'zh_pornsignals.json': 0}}
{'chain_length': '1', 'dir': 'jupyter-structured-clean-dedup', 'id': '359980', 'license': 'MIT', 'n_tokens_mistral': 7260, 'n_tokens_neox': 6613, 'n_words': 1252, 'path': 'domainSimulation_test2.ipynb', 'repo_name': 'andersdot/dust', 'size': '17075'}