Datasets used in this example are a system of hard hexagons, simulated in the NVT thermodynamic ensemble in HOOMD-Blue, for a dense fluid (phi065) and a solid (phi075)
from bokeh.io import output_notebook
output_notebook()
from bokeh.models import Legend
from bokeh.plotting import figure, output_file, show
from bokeh.layouts import gridplot
import numpy as np
import time
from freud import parallel
parallel.setNumThreads(4)
# define vertices for hexagons
verts = [[0.537284965911771, 0.31020161970069976],
[3.7988742065678664e-17, 0.6204032394013997],
[-0.5372849659117709, 0.31020161970070004],
[-0.5372849659117711, -0.31020161970069976],
[-1.1396622619703597e-16, -0.6204032394013997],
[0.5372849659117711, -0.3102016197006997]]
verts = np.array(verts)
# define colors for our system
c_list = ["#30A2DA", "#FC4F30", "#E5AE38", "#6D904F", "#9757DB",
"#188487", "#FF7F00", "#9A2C66", "#626DDA", "#8B8B8B"]
c_dict = dict()
c_dict[6] = c_list[0]
c_dict[5] = c_list[1]
c_dict[4] = c_list[2]
c_dict[3] = c_list[7]
c_dict[7] = c_list[4]
def default_bokeh(p):
"""
wrapper which takes the default bokeh outputs and changes them to more sensible values
"""
p.title.text_font_size = "18pt"
p.title.align = "center"
p.xaxis.axis_label_text_font_size = "14pt"
p.yaxis.axis_label_text_font_size = "14pt"
p.xaxis.major_tick_in = 10
p.xaxis.major_tick_out = 0
p.xaxis.minor_tick_in = 5
p.xaxis.minor_tick_out = 0
p.yaxis.major_tick_in = 10
p.yaxis.major_tick_out = 0
p.yaxis.minor_tick_in = 5
p.yaxis.minor_tick_out = 0
p.xaxis.major_label_text_font_size = "12pt"
p.yaxis.major_label_text_font_size = "12pt"
def cubeellipse(theta, lam=0.5, gamma=1., s=4.0, r=1., h=1.):
"""Create an RGB colormap from an input angle theta. Takes lam (a list of
intensity values, from 0 to 1), gamma (a nonlinear weighting power),
s (starting angle), r (number of revolutions around the circle), and
h (a hue factor)."""
import numpy
lam = lam**gamma
a = h*lam*(1 - lam)*.5
v = numpy.array([[-.14861, 1.78277], [-.29227, -.90649], [1.97294, 0.]], dtype=numpy.float32)
ctarray = numpy.array([numpy.cos(theta*r + s), numpy.sin(theta*r + s)], dtype=numpy.float32)
# convert to 255 rgb
ctarray = (lam + a*v.dot(ctarray)).T
ctarray *= 255
ctarray = ctarray.astype(dtype=np.int32)
return ctarray
def local_to_global(verts, positions, orientations):
"""
Take a list of vertices, positions, and orientations and create
a list of vertices in the "global coordinate system" for plotting
in bokeh
"""
num_particles = len(positions)
num_verts = len(verts)
# create list of vertices in the "local reference frame" i.e.
# centered at (0,0)
l_verts = np.zeros(shape=(num_particles, num_verts, 2), dtype=np.float32)
l_verts[:] = verts
# create array of rotation matrices
rot_mat = np.zeros(shape=(num_particles, 2, 2), dtype=np.float32)
for i, theta in enumerate(orientations):
rot_mat[i] = [[np.cos(theta), -np.sin(theta)], [np.sin(theta), np.cos(theta)]]
# rotate; uses einsum for speed; please see numpy documentation
# for more information
r_verts = np.einsum("lij,lkj->lki", rot_mat, l_verts)
# now translate to global coordinates
# need to create a position array with same shape as vertex array
l_pos = np.zeros(shape=(num_particles, num_verts, 2), dtype=np.float32)
for i in range(num_particles):
for j in range(len(verts)):
l_pos[i,j] = positions[i]
# translate
output_array = np.add(r_verts, l_pos)
return output_array
def clamp(x):
"""
limit values between 0 and 255
http://stackoverflow.com/questions/3380726/converting-a-rgb-color-tuple-to-a-six-digit-code-in-python
"""
return max(0, min(x, 255))
One of the basic computations required for higher-level computations (such as the [hexatic order parameter](AIChE Hexatic.ipynb)) is finding the nearest neighbors of a particle. This tutorial will show you how to compute the nearest neighbors and visualize that data.
The algorithm is straightforward:
for each particle i:
for each particle j in neighbor_cells(i):
r_ij = position[j] - position[i]
r = sqrt(dot(r_ij, r_ij))
l_r_array.append(r)
l_n_array.append(j)
# sort by distance
sort(n_array, r_array)
neighbor_array[i] = n_array[:k]
# import the box and locality modules;
# Nearest Neighbors calculation is provided by the locality module
from freud import box, locality
# load the data
data_path = "ex_data/phi065"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]
Before proceeding, we should probably view our system first. Freud does not make any assumptions about your data and is not specifically designed for any one visualization package. Here we use bokeh to render our system. Bokeh is not appropriate for real-time interaction with your simulation data, nor is it appropriate for 3D data, but is perfectly fine for rendering individual simulation frames, so we will use it here
# grab data from last frame
l_box = box_data[-1]
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:,3]), np.copy(l_quat[:,0]))
# create box
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
side_length = max(fbox.Lx, fbox.Ly)
l_min = -side_length / 2.0
l_min *= 1.1
l_max = -l_min
# take local vertices and rotate, translate into system coordinates
patches = local_to_global(verts, l_pos[:,0:2], l_ang)
# plot
p = figure(title="System Visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),
fill_color=(42,126,187), line_color="black", line_width=1.5, legend="hexagons")
# box display
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
By eye, we can see regions where the hexagons appear close-packed, as well as regions where there are vacancies. We will be using the Nearest Neighbor object to investigate this in our system
This module will give the indices of the $k$ particles which are nearest to another particle. Freud provides two different modes by which to compute the nearest neighbors, selected by the strict_cut
variable:
strict_cut=False
(default): The value for rmax
is expanded until every particle has $k$ nearest neighborsstrict_cut=True
: the rmax
value is not expanded, so that any "vacancies" in the number of neighbors found are filled with UINTMAX
strict_cut=False
¶First we show how to use the strict_cut=False
mode to find the neighbors of a specific particle
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 6
# create freud nearest neighbors object
nn = locality.NearestNeighbors(rmax=1.5,n_neigh=n_neigh,strict_cut=False)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = nn.getNeighborList()
# get the number of particles
num_particles = nn.getNRef()
# get the neighbors for particle 0
pidx = 0
n_idxs = n_list[pidx]
# get position, orientation for the central particle
center_pos = np.zeros(shape=(1, 3),dtype=np.float32)
center_ang = np.zeros(shape=(1),dtype=np.float32)
center_pos[:] = l_pos[pidx]
center_ang[:] = l_ang[pidx]
# get the positions, orientations for the neighbor particles
neigh_pos = np.zeros(shape=(n_neigh, 3),dtype=np.float32)
neigh_ang = np.zeros(shape=(n_neigh),dtype=np.float32)
neigh_pos[:] = l_pos[n_idxs]
neigh_ang[:] = l_ang[n_idxs]
# render in bokeh
# create array of transformed positions
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)
# turn into list of colors
# bokeh (as of this version) requires hex colors, so convert rgb to hex
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
# plot
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="centers")
# box display
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.023049116134643555
Notice that nearest neighbors properly handles periodic boundary conditions
We do the same thing below, but for a particle/neighbors not spanning the box
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 6
# create freud nearest neighbors object
nn = locality.NearestNeighbors(rmax=1.5,n_neigh=n_neigh,strict_cut=False)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = nn.getNeighborList()
# get the number of particles
num_particles = nn.getNRef()
# get the neighbors for particle 1000
pidx = 1000
n_idxs = n_list[pidx]
# get position, orientation for the central particle
center_pos = np.zeros(shape=(1, 3),dtype=np.float32)
center_ang = np.zeros(shape=(1),dtype=np.float32)
center_pos[:] = l_pos[pidx]
center_ang[:] = l_ang[pidx]
# get positions, orientations for the neighbor particles
neigh_pos = np.zeros(shape=(n_neigh, 3),dtype=np.float32)
neigh_ang = np.zeros(shape=(n_neigh),dtype=np.float32)
neigh_pos[:] = l_pos[n_idxs]
neigh_ang[:] = l_ang[n_idxs]
non_neigh_pos = np.zeros(shape=(1, 3),dtype=np.float32)
non_neigh_ang = np.zeros(shape=(1),dtype=np.float32)
non_neigh_pos[:] = neigh_pos[-1]
non_neigh_ang[:] = neigh_ang[-1]
# render in bokeh
# create array of transformed positions
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)
non_n_patches = local_to_global(verts, non_neigh_pos[:,0:2], non_neigh_ang)
# turn into list of colors
# bokeh (as of this version) requires hex colors, so convert rgb to hex
center_color = np.array([c_list[0] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
non_neigh_color = np.array([c_list[-2] for _ in range(non_neigh_pos.shape[0])])
# color the last particle differently
# plot
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
p.patches(xs=non_n_patches[:,:,0].tolist(), ys=non_n_patches[:,:,1].tolist(),
fill_color=non_neigh_color.tolist(), line_color="black", legend="non-neighbor")
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="centers")
# box display
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.02012920379638672
Notice that while Freud found the 6 nearest neighbors, one of the particles isn't really in a neighbor position (which we have colored purple). How do we go about finding particles with a deficit or surplus of neighbors?
strict_cut=True
¶Now for strict_cut=True
. This mode allow you to find particles which have fewer than the specified number of particles. For this system, we'll search for 8 neighbors, so that we can display particles with both a deficit and a surplus of neighbors.
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 8
# create freud nearest neighbors object
# set rmax to some value
rmax=1.7
nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = np.copy(nn.getNeighborList())
# get the number of particles
num_particles = nn.getNRef()
# now for array manipulation magic
# create an integer array of the same shape as the neighbor list array
int_arr = np.ones(shape=n_list.shape, dtype=np.int32)
# "search" for non-indexed particles (missing neighbors)
# while it would be most accurate to use the UINTMAX value
# provided by nn.getUINTMAX(), but this works just as well
int_arr[n_list > (num_particles-1)] = 0
# sum along particle index axis to determine the number of neighbors per particle
n_neighbors = np.sum(int_arr, axis=1)
# find the complement (if desired) to find number of missing neighbors per particle
n_deficits = n_neigh - n_neighbors
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
for k in np.unique(n_neighbors):
# find particles with k neighbors
c_idxs = np.copy(np.where(n_neighbors==k)[0])
center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)
center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.004876852035522461
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 8
# create freud nearest neighbors object
# set rmax to some value
rmax=1.7
nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = np.copy(nn.getNeighborList())
# get the number of particles
num_particles = nn.getNRef()
# now for array manipulation magic
# create an integer array of the same shape as the neighbor list array
int_arr = np.ones(shape=n_list.shape, dtype=np.int32)
# "search" for non-indexed particles (missing neighbors)
# while it would be most accurate to use the UINTMAX value
# provided by nn.getUINTMAX(), but this works just as well
int_arr[n_list > (num_particles-1)] = 0
# sum along particle index axis to determine the number of neighbors per particle
n_neighbors = np.sum(int_arr, axis=1)
# find the complement (if desired) to find number of missing neighbors per particle
n_deficits = n_neigh - n_neighbors
for k in np.unique(n_neighbors):
p = figure(title="Nearest Neighbors: k={}".format(k),
x_range=(l_min, l_max), y_range=(l_min, l_max))
# find particles with k neighbors
c_idxs = np.copy(np.where(n_neighbors==k)[0])
center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)
center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
# take local vertices and rotate, translate into system coordinates
patches = local_to_global(verts, l_pos[:,0:2], l_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),
fill_color=(128,128,128,0.1), line_color=(0,0,0,0.1), legend="other")
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.0025110244750976562
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 8
# create freud nearest neighbors object
# set rmax to some value
rmax=1.7
nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = np.copy(nn.getNeighborList())
# get the number of particles
num_particles = nn.getNRef()
# now for array manipulation magic
# create an integer array of the same shape as the neighbor list array
int_arr = np.ones(shape=n_list.shape, dtype=np.int32)
# "search" for non-indexed particles (missing neighbors)
# while it would be most accurate to use the UINTMAX value
# provided by nn.getUINTMAX(), but this works just as well
int_arr[n_list > (num_particles-1)] = 0
# sum along particle index axis to determine the number of neighbors per particle
n_neighbors = np.sum(int_arr, axis=1)
# find the complement (if desired) to find number of missing neighbors per particle
n_deficits = n_neigh - n_neighbors
for k in np.unique(n_neighbors):
p = figure(title="Nearest Neighbors: k={}".format(k),
x_range=(l_min, l_max), y_range=(l_min, l_max))
# find particles with k neighbors
c_idxs = np.copy(np.where(n_neighbors==k)[0])
center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)
center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
neigh_pos = np.zeros(shape=(k*len(c_idxs), 3),dtype=np.float32)
neigh_ang = np.zeros(shape=(k*len(c_idxs)),dtype=np.float32)
for i, pidx in enumerate(c_idxs):
# create a list of positions, angles to draw
n_idxs = n_list[pidx,:k]
for j, nidx in enumerate(n_idxs):
neigh_pos[k*i+j] = l_pos[nidx]
neigh_ang[k*i+j] = l_ang[nidx]
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
n_patches = local_to_global(verts, neigh_pos[:,0:2], neigh_ang)
patches = local_to_global(verts, l_pos[:,0:2], l_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
neigh_color = np.array([c_list[-1] for _ in range(neigh_pos.shape[0])])
p.patches(xs=patches[:,:,0].tolist(), ys=patches[:,:,1].tolist(),
fill_color=(128,128,128,0.1), line_color=(0,0,0,0.1), legend="other")
p.patches(xs=n_patches[:,:,0].tolist(), ys=n_patches[:,:,1].tolist(),
fill_color=neigh_color.tolist(), line_color="black", legend="neighbors")
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.002343893051147461
Visualize in the same way, but with a solid system
data_path = "ex_data/phi075"
box_data = np.load("{}/box_data.npy".format(data_path))
pos_data = np.load("{}/pos_data.npy".format(data_path))
quat_data = np.load("{}/quat_data.npy".format(data_path))
n_frames = pos_data.shape[0]
# grab data from last frame
l_box = box_data[-1]
l_pos = pos_data[-1]
l_quat = quat_data[-1]
l_ang = 2*np.arctan2(np.copy(l_quat[:,3]), np.copy(l_quat[:,0]))
# create box
fbox = box.Box(Lx=l_box["Lx"], Ly=l_box["Ly"], is2D=True)
side_length = max(fbox.Lx, fbox.Ly)
l_min = -side_length / 2.0
l_min *= 1.1
l_max = -l_min
# create freud nearest neighbor object
# set number of neighbors
n_neigh = 8
# create freud nearest neighbors object
# set rmax to some value
rmax=1.65
nn = locality.NearestNeighbors(rmax=rmax,n_neigh=n_neigh,strict_cut=True)
# compute nearest neighbors for 6 nearest neighbors
start_time = time.time()
nn.compute(fbox, l_pos, l_pos)
stop_time = time.time()
print("time to calc 1 frame = {}".format(stop_time-start_time))
# get the neighborlist
n_list = np.copy(nn.getNeighborList())
# get the number of particles
num_particles = nn.getNRef()
# now for array manipulation magic
# create an integer array of the same shape as the neighbor list array
int_arr = np.ones(shape=n_list.shape, dtype=np.int32)
# "search" for non-indexed particles (missing neighbors)
# while it would be most accurate to use the UINTMAX value
# provided by nn.getUINTMAX(), but this works just as well
int_arr[n_list > (num_particles-1)] = 0
# sum along particle index axis to determine the number of neighbors per particle
n_neighbors = np.sum(int_arr, axis=1)
# find the complement (if desired) to find number of missing neighbors per particle
n_deficits = n_neigh - n_neighbors
p = figure(title="Nearest Neighbors visualization",x_range=(l_min, l_max), y_range=(l_min, l_max))
for k in np.unique(n_neighbors):
# find particles with k neighbors
c_idxs = np.copy(np.where(n_neighbors==k)[0])
center_pos = np.zeros(shape=(len(c_idxs), 3),dtype=np.float32)
center_ang = np.zeros(shape=(len(c_idxs)),dtype=np.float32)
center_pos = l_pos[c_idxs]
center_ang = l_ang[c_idxs]
c_patches = local_to_global(verts, center_pos[:,0:2], center_ang)
center_color = np.array([c_dict[k] for _ in range(center_pos.shape[0])])
p.patches(xs=c_patches[:,:,0].tolist(), ys=c_patches[:,:,1].tolist(),
fill_color=center_color.tolist(), line_color="black", legend="k={}".format(k))
p.patches(xs=[[-fbox.Lx/2, fbox.Lx/2, fbox.Lx/2, -fbox.Lx/2]],
ys=[[-fbox.Ly/2, -fbox.Ly/2, fbox.Ly/2, fbox.Ly/2]],
fill_color=(0,0,0,0), line_color="black", line_width=2)
p.legend.location='bottom_center'
p.legend.orientation='horizontal'
default_bokeh(p)
show(p)
time to calc 1 frame = 0.0023381710052490234