In [6]:
%pylab inline
In [7]:
import readhaloHDF5
import snapHDF5 as snap
import readsubfHDF5 as sub
from scipy.interpolate import griddata
In [8]:
rcParams['figure.figsize'] = 10, 10
In [9]:
def enlarge_labels(ax,nsize=25):
ticklabels_x = ax.get_xticklabels()
ticklabels_y = ax.get_yticklabels()
for label_x in ticklabels_x:
label_x.set_fontsize(nsize)
label_x.set_family('serif')
for label_y in ticklabels_y:
label_y.set_fontsize(nsize)
label_y.set_family('serif')
In [10]:
def load_snap_dm(path, snapnum, filenum):
filename="%s/snap_%d.%d"%(path, snapnum, filenum)
position = snap.read_block(filename, "POS ", parttype=0)
velocity = snap.read_block(filename, "VEL ", parttype=0)
return position, velocity
In [11]:
def load_snap(path, snapnum, filenum):
filename="%s/snap_%d.%d"%(path, snapnum, filenum)
position = snap.read_block(filename, "POS ", parttype=0)
velocity = snap.read_block(filename, "VEL ", parttype=0)
rho_gas = snap.read_block(filename, "RHO ", parttype=0)
volume = snap.read_block(filename, "VOL ", parttype=0)
vgrad = snap.read_block(filename, "VGRD", parttype=0)
return position, velocity, rho_gas, volume, vgrad
In [12]:
def make_dens_plot(density, snapnum=1, filenum=1):
fig = figure(1, figsize=(9.5,6.5))
min_val = amin(log10(density))
n = hist(log10(dens),50, normed=1, color='silver')
ax = axes()
enlarge_labels(ax)
ax.set_ylabel("${\mathrm{Histogram}}$", fontsize=25)
ax.set_xlabel(r"$\log_{10} \rho [10^{10}M_{\odot} \mathrm{kpc}^{-3} h^{-2}]$",fontsize=25)
filename = 'dens_hist_snap%d_file%d'%(snapnum,filenum)
texto = "${\mathrm{file\ %d}}$"%(filenum)
print texto
ax.set_xlim([-12.0,-5.0])
ax.set_ylim([0.0,1.0])
ax.text(-11.5,0.8,texto, fontsize=35)
savefig(filename + '.eps',format = 'eps', transparent=True)
savefig(filename + '.pdf',format = 'pdf', transparent=True)
savefig(filename + '.png',format = 'png', transparent=True)
close(fig)
In [17]:
def dm_h5_to_ascii(path, snapnum, filenum):
filename="%s/snap_%d.%d"%(path, snapnum, filenum)
position = snap.read_block(filename, "POS ", parttype=1)
velocity = snap.read_block(filename, "VEL ", parttype=1)
filenameout = "%s/snap_%d.%d.dm_pos_ascii"%(path, snapnum, filenum)
out = open(filenameout,"w")
n_points = size(position[:,0])
for i in range(n_points):
out.write("%f %f %f\n"%(position[i,0], position[i,1], position[i,2]))
out.close()
filenameout = "%s/snap_%d.%d.dm_vel_ascii"%(path, snapnum, filenum)
out = open(filenameout,"w")
n_points = size(position[:,0])
for i in range(n_points):
out.write("%f %f %f\n"%(position[i,0], position[i,1], position[i,2]))
out.close()
return
In [18]:
snapnum=136
filenum=0
data_path= "/Users/forero/work/ArepoWeb/data/L75n455FP/snapdir_136"
dm_h5_to_ascii(data_path, snapnum, filenum)
In [35]:
snapnum=136
filenum=0
data_path= "/Users/forero/work/ArepoWeb/data/L75n455FP/snapdir_136"
pos, vel, dens, vol, v_gradient = load_snap(data_path, snapnum, filenum)
In [36]:
pos_x = pos[:,0]
pos_y = pos[:,1]
pos_z = pos[:,2]
In [37]:
index = where((pos_z>40000.0) & (pos_z<49000.0))
index = index[0]
scatter(pos_x[index], pos_y[index],s=0.1)
Out[37]:
In [38]:
n = hist((abs(v_gradient[index,8])),10)
In [39]:
scatter((v_gradient[index,0]), (v_gradient[index,4]))
xlabel(r"$\frac{\partial v_x}{\partial x}[\mathrm{arepo units}]$", fontsize=25)
ylabel(r"$\frac{\partial v_y}{\partial y}[\mathrm{arepo units}]$", fontsize=25)
Out[39]:
In [44]:
scatter((log10(dens[index])),(v_gradient[index,4] + v_gradient[index,0] +v_gradient[index,8]),s=0.1)
xlabel(r"$\log_{10}\rho_{\rm gas}\ [\mathrm{arepo units}]$", fontsize=25)
ylabel(r"$\nabla\cdot \vec{v}\ [\mathrm{arepo units}]$", fontsize=25)
Out[44]:
In [54]:
snapnum=136
slice_x = empty((0))
slice_y = empty((0))
slice_dens = empty((0))
slice_v_grad = empty((0,0))
data_path= "/Users/forero/work/ArepoWeb/data/L75n455FP/snapdir_136"
for i in range(32):
pos, vel, dens, vol, v_grad = load_snap(data_path, snapnum, i)
index = where((pos[:,2]>35000.0) & (pos[:,2]<35100.0))
index = index[0]
slice_x = append(slice_x, pos[index,0])
slice_y = append(slice_y, pos[index,1])
slice_dens = append(slice_dens, dens[index])
slice_v_grad = append(slice_v_grad , v_grad[index,0] + v_grad[index,4]+v_grad[index,8])
print size(index)
print size(slice_x), size(slice_dens)
In [55]:
min_y = amin(slice_y)
max_y = amax(slice_y)
min_x = amin(slice_x)
max_x = amax(slice_x)
grid_x, grid_y = mgrid[min_x:max_x:500j, min_y:max_y:500j]
print median(slice_dens), min(slice_dens), max(slice_dens)
In [56]:
n_points = size(slice_x)
points = ones((n_points,2))
points[:,0] = slice_x
points[:,1] = slice_y
grid_d = griddata(points, log10(slice_dens), (grid_x, grid_y), method='nearest')
In [57]:
imshow((grid_d.T), extent=(min_x,max_x,min_y,max_y), aspect='equal',origin='lower')
print amin(grid_d), amax(grid_d)
In [63]:
print (slice_v_grad)
grid_d = griddata(points, log10(+(slice_v_grad)), (grid_x, grid_y), method='nearest')
imshow((grid_d.T), extent=(min_x,max_x,min_y,max_y), aspect='equal',origin='lower')
Out[63]: