In [1]:
%pylab inline
In [2]:
import readhaloHDF5
import snapHDF5 as snap
from scipy.interpolate import griddata
from struct import *
import h5py
In [3]:
rcParams['figure.figsize'] = 8, 8
In [4]:
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 [5]:
def load_snap(path, snapnum, filenum):
filename="%s/snap_%d.%d"%(path, snapnum, filenum)
position = snap.read_block(filename, "POS ", parttype=0) #position
rho = snap.read_block(filename, "RHO ", parttype=0) # density
volume = snap.read_block(filename, "VOL ", parttype=0) # volume
return position, rho, volume
In [6]:
def load_gas_web(path, snapnum, filenum):
filename="%s/gasweb_snap_%d.%d.hdf5"%(path, snapnum, filenum)
fin = h5py.File(filename, "r")
eigenvals = fin["ShearTensor"]["Eigenvalues"]
eigenvecs = fin["ShearTensor"]["Eigenvectors"]
helicity = fin["ShearTensor"]["Helicity"]
vorticity = fin["ShearTensor"]["Vorticity"]
return eigenvals, eigenvecs, helicity, vorticity
In [7]:
def read_CIC_scalar(filename):
f = open(filename, "rb")
dumb = f.read(38)
dumb = f.read(4)
n_x = f.read(4)
n_y = f.read(4)
n_z = f.read(4)
nodes = f.read(4)
x0 = f.read(4)
y0 = f.read(4)
z0 = f.read(4)
dx = f.read(4)
dy = f.read(4)
dz = f.read(4)
dumb = f.read(4)
n_x = (unpack('i', n_x))[0]
n_y = (unpack('i', n_y))[0]
n_z = (unpack('i', n_z))[0]
nodes = (unpack('i', nodes))[0]
dx = (unpack('f', dx))[0]
dy = (unpack('f', dy))[0]
dz = (unpack('f', dz))[0]
x0 = (unpack('f', x0))[0]
y0 = (unpack('f', y0))[0]
z0 = (unpack('f', z0))[0]
print n_x, n_y, n_z, nodes, dx, dy, dz
total_nodes = n_x * n_y *n_z
dumb = f.read(4)
array_data = f.read(total_nodes*4)
dumb = f.read(4)
format_s = str(total_nodes)+'f'
array_data = unpack(format_s, array_data)
f.close()
array_data = np.array(array_data)
array_data.resize(n_z,n_y,n_x)
array_data = array_data.transpose()
return array_data
In [8]:
snapnum=136
snap_data_path= "/Users/forero/work/ArepoWeb/data/L75n455FP/snapdir_%d/"%(snapnum)
web_data_path= "/Users/forero/work/ArepoWeb/data/L75n455FP/gasweb_%d/"%(snapnum)
pos_gas, rho_gas, vol_gas = load_snap(snap_data_path, snapnum, 0)
eigenvalues_gas, eigenvectors_gas, helicity_gas, vorticity_gas = load_gas_web(web_data_path, snapnum, 0)
center_slab = 30000.0
delta_slab = 200.0
index = where((pos_gas[:,2]>(center_slab - delta_slab*0.5))&(pos_gas[:,2]<(center_slab+delta_slab*0.5)))
index = index[0]
pos_gas = pos_gas[index,:]
rho_gas = rho_gas[index,:]
vol_gas = vol_gas[index,:]
eigenvalues_gas = eigenvalues_gas[index,:]
helicity_gas = helicity_gas[index,:]
vorticity_gas = vorticity_gas[index,:]
print size(index)
for i in range(1,32):
print i
pos_gas_tmp, rho_gas_tmp, vol_gas_tmp = load_snap(snap_data_path, snapnum, i)
eigenvalues_gas_tmp, eigenvectors_gas_tmp, helicity_gas_tmp, vorticity_gas_tmp = \
load_gas_web(web_data_path, snapnum, i)
index = where((pos_gas_tmp[:,2]>(center_slab - delta_slab*0.5))&
(pos_gas_tmp[:,2]<(center_slab+delta_slab*0.5)))
index = index[0]
pos_gas = append(pos_gas, pos_gas_tmp[index,:], axis=0)
rho_gas = append(rho_gas, rho_gas_tmp[index,:], axis=0)
vol_gas = append(vol_gas, vol_gas_tmp[index,:], axis=0)
eigenvalues_gas = append(eigenvalues_gas, eigenvalues_gas_tmp[index,:], axis=0)
eigenvectors_gas = append(eigenvectors_gas, eigenvectors_gas_tmp[index,:], axis=0)
helicity_gas = append(helicity_gas, helicity_gas_tmp[index,:], axis=0)
vorticity_gas = append(vorticity_gas, vorticity_gas_tmp[index,:], axis=0)
print shape(eigenvalues_gas)
In [9]:
def interpolate_grid(pos, to_interpolate):
n_points = size(pos[:,0])
points = ones((n_points,2))
points[:,0] = pos[:,0]
points[:,1] = pos[:,1]
#to_interpolate = log10(abs(helicity_gas[:,0]))
min_y = amin(points[:,0])
max_y = amax(points[:,0])
min_x = amin(points[:,1])
max_x = amax(points[:,1])
grid_x, grid_y = mgrid[min_x:max_x:1000j, min_y:max_y:1000j]
grid_d = griddata(points, to_interpolate, (grid_x, grid_y), method='nearest')
return grid_d
In [10]:
grid_helicity = interpolate_grid(pos_gas, log10(abs(helicity_gas[:,0])))
In [11]:
imshow((grid_helicity.T), aspect='equal',origin='lower')
Out[11]:
In [12]:
in_grid = log10(sqrt(vorticity_gas[:,0]**2 + vorticity_gas[:,1]**2 + vorticity_gas[:,2]**2))
grid_vorticity = interpolate_grid(pos_gas, in_grid)
In [13]:
imshow((grid_vorticity[:300,:300].T), aspect='equal',origin='lower')
Out[13]:
In [14]:
in_grid = log10(rho_gas)
grid_density = interpolate_grid(pos_gas, in_grid)
In [15]:
imshow((grid_density.T), aspect='equal',origin='lower')
imshow((grid_vorticity.T), aspect='equal',origin='lower' ,alpha=0.4)
Out[15]:
In [16]:
in_grid = log10((eigenvalues_gas[:,0] + eigenvalues_gas[:,1] + eigenvalues_gas[:,2]))
grid_gradient = interpolate_grid(pos_gas, in_grid)
In [17]:
imshow((grid_gradient.T), aspect='equal',origin='lower')
Out[17]:
In [60]:
divergence_to_vorticity = abs(eigenvalues_gas[:,0] + eigenvalues_gas[:,1] + eigenvalues_gas[:,2])/\
sqrt(vorticity_gas[:,0]**2 + vorticity_gas[:,1]**2 + vorticity_gas[:,2]**2)
in_grid = log10(divergence_to_vorticity-2)
grid_div_to_vor = interpolate_grid(pos_gas, in_grid)
In [61]:
imshow((grid_div_to_vor.T), aspect='equal',origin='lower')
Out[61]:
In [47]:
print shape(in_grid)
index = where(~(isnan(in_grid)))
in_grid = in_grid[index]
print amin(in_grid), amax(in_grid)
n = hist(in_grid,bins=200)
In [62]:
environment_grid = eigenvalues_gas[:,0].copy()
PEAK = 3
FILAMENT = 2
SHEET = 1
VOID =0
lambda_th = 0.02 #the interesting values is 0.02
peak_id = where(eigenvalues_gas[:,2]>lambda_th)
peak_id = peak_id[0]
filament_id = where((eigenvalues_gas[:,1]>lambda_th) & (eigenvalues_gas[:,2]<lambda_th))
filament_id = filament_id[0]
sheet_id = where((eigenvalues_gas[:,0]>lambda_th) & (eigenvalues_gas[:,1]<lambda_th))
sheet_id = sheet_id[0]
void_id = where(eigenvalues_gas[:,0]<lambda_th)
void_id = void_id[0]
environment_grid[peak_id] = PEAK
environment_grid[filament_id] = FILAMENT
environment_grid[sheet_id] = SHEET
environment_grid[void_id] = VOID
print size(sheet_id), size(peak_id), size(filament_id), size(void_id)
print size(sheet_id) + size(peak_id) + size(filament_id) +size(void_id)
In [63]:
grid_env = interpolate_grid(pos_gas, environment_grid)
shape(grid_env)
Out[63]:
In [64]:
imshow((grid_env.T), aspect='equal',origin='lower')
Out[64]:
In [65]:
n_size = 128
lambda_1 = read_CIC_scalar("../data/L75n455FP/Tweb/%d/snap_136.s1.00.eigen_1"%(n_size))
lambda_2 = read_CIC_scalar("../data/L75n455FP/Tweb/%d/snap_136.s1.00.eigen_2"%(n_size))
lambda_3 = read_CIC_scalar("../data/L75n455FP/Tweb/%d/snap_136.s1.00.eigen_3"%(n_size))
delta = lambda_1 + lambda_2 + lambda_3
In [66]:
cut_z = int(n_size*center_slab/70000.0 - 0.5)
In [67]:
imshow((grid_vorticity.T), aspect='equal',origin='lower' ,alpha=1.0, cmap='Greys')
imshow(log10(delta[:,:,cut_z]+1.56).T, extent=(0,1000,0,1000),aspect='equal',origin='lower', alpha=0.3)
Out[67]:
In [119]:
VOID = 0
FILAMENT = 1
SHEET = 2
PEAK = 3
lambda_th_DM = 0.2
environment = delta.copy()
peak_id = where(lambda_3>lambda_th_DM)
environment[peak_id] = PEAK
filament_id = where((lambda_2>lambda_th_DM) & (lambda_3<lambda_th_DM))
environment[filament_id] = FILAMENT
sheet_id = where((lambda_1>lambda_th_DM) & (lambda_2<lambda_th_DM))
environment[sheet_id] = SHEET
void_id = where(lambda_1<lambda_th_DM)
environment[void_id] = VOID
In [120]:
imshow((grid_helicity.T), aspect='equal',origin='lower' ,alpha=1.0, cmap='Greys')
imshow(environment[:,:,cut_z].T, extent=(0,1000,0,1000),aspect='equal',origin='lower', alpha=0.3)
Out[120]:
In [121]:
imshow((grid_env.T**0.1), aspect='equal',origin='lower', cmap='Greys')
imshow(environment[:,:,cut_z].T, extent=(0,1000,0,1000),aspect='equal',origin='lower', alpha=0.5)
Out[121]:
In [122]:
index = where(environment==VOID)
density_in_voids = delta[index]
index = where(environment==SHEET)
density_in_sheets = delta[index]
In [123]:
histo_void, edges_void = histogram(density_in_voids)
histo_sheets, edges_sheets = histogram(density_in_sheets)
In [125]:
plot(edges_void[1:], log10(histo_void))
Out[125]: