In [9]:

    

import numpy as np
import matplotlib.pyplot as plt
import matplotlib.tri as tri
%matplotlib inline

# Yao-Yuan Mao's Gadget snapshot reader
from readGadgetSnapshot import readGadgetSnapshot


    

In [15]:

    

# USER OPTIONS

# snapshot name
snapfile = 'snap16'

# z-bounds of the particle slab
zmin = 10.0
zmax = 20.0

# output file names
svgfile = 'delaunay.svg'
epsfile = 'delaunay.eps'


    

In [16]:

    

# Loading and massaging the data

# get particle locations
header, pos = readGadgetSnapshot(snapfile, read_pos=True)
pos_x = pos.T[0]
pos_y = pos.T[1]
pos_z = pos.T[2]

# Add a boundary buffer to make the 
# Delaunay tessellation appear periodic
L = header.BoxSize
bfac = 0.1
# x-periodicity
lo_inds = np.where(pos_x < bfac*L)
hi_inds = np.where(pos_x >= bfac*L)
pos_x = np.append(pos_x, pos_x[lo_inds] + L)
pos_x = np.append(pos_x, pos_x[hi_inds] - L)
pos_y = np.append(pos_y, pos_y[lo_inds])
pos_y = np.append(pos_y, pos_y[hi_inds])
pos_z = np.append(pos_z, pos_z[lo_inds])
pos_z = np.append(pos_z, pos_z[hi_inds])
# y-periodicity
lo_inds = np.where(pos_y < bfac*L)
hi_inds = np.where(pos_y >= bfac*L)
pos_x = np.append(pos_x, pos_x[lo_inds])
pos_x = np.append(pos_x, pos_x[hi_inds])
pos_y = np.append(pos_y, pos_y[lo_inds] + L)
pos_y = np.append(pos_y, pos_y[hi_inds] - L)
pos_z = np.append(pos_z, pos_z[lo_inds])
pos_z = np.append(pos_z, pos_z[hi_inds])

# pick out a particle slab in the z-coordinate
slab_inds = np.where((pos_z >= zmin) & (pos_z < zmax))
pos_x = pos_x[slab_inds]
pos_y = pos_y[slab_inds]


    

In [17]:

    

# Make the delaunay triangulation
triang = tri.Triangulation(pos_x, pos_y)


    

In [18]:

    

# plot it up!
fig, ax1 = plt.subplots(1, 1, figsize=(15,15))
ax1.set_xlim(0,L)
ax1.set_ylim(0,L)
ax1.axis('off')
ax1.triplot(triang, '-')

# save the SVG figure (can import SVG and do postprocessing directly in LaTex)
fig.savefig(svgfile, format='svg',transparent=True, bbox_inches='tight', pad_inches=0)
fig.savefig(epsfile, format='eps')


    

In [ ]: