Skip to content

Create PenguinMotion3D_final #10

New issue

Have a question about this project? Sign up for a free GitHub account to open an issue and contact its maintainers and the community.

Sign up for GitHub

By clicking “Sign up for GitHub”, you agree to our terms of service and privacy statement. We’ll occasionally send you account related emails.

Already on GitHub? Sign in to your account

Open
slnhermans wants to merge 1 commit into
base: master
Choose a base branch
from
Open

Create PenguinMotion3D_final #10

Changes from all commits
Commits
Show all changes
1 commit
Select commit
File filter

Filter by extension

Filter by extension
Conversations
Failed to load comments.
Loading
Jump to
Jump to file
Failed to load files.
Loading
Diff view
Diff view
253 changes: 253 additions & 0 deletions PenguinMotion3D_final
View file
Open in desktop
Original file line number Diff line number Diff line change
@@ -0,0 +1,253 @@
import numpy as np
import scipy as sp
import random
import matplotlib.pyplot as plt
from mpl_toolkits.mplot3d import Axes3D
from pylab import *
#from matplotlib.backends.backend_pdf import PdfPages
#pp = PdfPages('16_30frames.pdf')

global a, dev_a, N_xy, N_z, critical_radius, neighbour_radius,neighbour_radius_far

# Computational Parameters
a = 1. # Average size of the penguins
dev_a = a/10. # Deviation of the average size

N_xy = 6 # Number of penguins in the x and y direction
N_z = 11 # Number of penguins in the z direction

dev_i = 0.01
dev_j = 0.01
dev_or = 0.25*np.pi

critical_radius = a
neighbour_radius = 2.2 * a
neighbour_radius_far = 6. * a

t_final = 37
start_drawing = 0

# Physical Constants
F_self = 0.01
F_in = 1.
k = 1.

T_in = 3.
T_n = 0.03
T_align = 3.

class Penguin(object):

def __init__(self, radius, position, alignment, boundary):
"""
Arguments:
radius := float representing penguin radius
position := numpy array representing penguin position (x, y, z)
alignment := numpy array representing alignment unit vector
boundary := boolean value, True if boundary penguin, False if bulk penguin
"""

self.radius = radius
self.position = position
self.alignment = alignment
self.boundary = boundary

def get_distance(self, penguin2):

return self.position - penguin2.position
#return np.dot((self.position - penguin2.position),(self.position - penguin2.position))

def net_force(self, F_self, F_in, k, penguin_list, a):
"""
Arguments:
F_self := multiplicative strength factor of penguin self-propulsion
F_in := multiplicative strength factor of the penguin boundary force
k := spring constant relating to the repulsive force of overlapping penguins
penguin_list := python list containing all penguins in the system
a := average radius of all penguins in the system
"""
abs_align = np.dot(self.alignment,self.alignment)
F_selfPropulsion = F_self * self.alignment/abs_align

F_repulsion = 0

if (self.boundary == True):

F_boundary = F_in * self.alignment/abs_align

else:

F_boundary = 0

above_plane = []
r_nb = np.zeros(3,dtype=float)
neigh_nr = 0
for i in range(len(penguin_list)):

if (penguin_list[i] != self):

r = self.get_distance(penguin_list[i])
r_mag = np.linalg.norm(r)
if (r_mag < neighbour_radius_far):
r_nb += r
if (r_mag < neighbour_radius):
neigh_nr += 1
#Repulstion force if neighbours are too close
if (r_mag < critical_radius):
F_repulsion -= k * r
#Boundary determination
#alignment [a b c], position [x_0, y_0, z_0], neighbour at [x, y, z]
#Plane defined by: a(x-x_0)+b(y-y_0)+c(z-z_0)=0
plane = self.alignment[0]*(penguin_list[i].position[0] - self.position[0]) + self.alignment[1]*(penguin_list[i].position[1] - self.position[1]) + self.alignment[2]*(penguin_list[i].position[2] - self.position[2])
#print(i, plane)
if plane > 0:
above_plane.append(0)
else:
above_plane.append(1)

#neigh_nr = len(above_plane)
#if (all(j == 0 for j in above_plane)) or (all (j == 1 for j in above_plane)):
if neigh_nr < 3:
self.boundary = True
F_boundary = F_in * -r_nb
else:
self.boundary = False
return F_boundary + F_repulsion + F_selfPropulsion

def net_torque(self, T_in, T_n, T_align, penguin_list, a):
"""
Arguments:
T_in := multiplicative strength factor of the boundary torque term
T_n := multiplicative strength factor of the random torque term
T_align := multiplicative strength factor of the alignment torque term
penguin_list := python list containing all penguins in the system
a := average radius of all penguins in the system
"""

critical_radius = a


# Torque Alignment

align = np.zeros(3, dtype=float)
T_b = np.zeros(3, dtype=float)

for i in range(len(penguin_list)):
if (penguin_list[i] != self):

r = self.get_distance(penguin_list[i])
r_mag = np.linalg.norm(r)

if r_mag < neighbour_radius:
align += penguin.alignment
T_b -= r

T_alignment = T_align * align

if (self.boundary == True):
T_boundary = T_in * T_b


## distance_sum = np.zeros(3)
##
## for i in range(len(penguin_list)):
##
## r = self.get_distance(penguin_list[i])
## r_mag = np.linalg.norm(r)
##
## if (penguin_list[i] != self) and (penguin_list[i].boundary == True):
##
## distance_sum += np.linalg.norm(self.get_distance(penguin_list[i]))
##
## if (penguin_list[i] != self) and (r_mag < critical_radius):
##
## T_alignment += T_align * (self.alignment - penguin_list[i].alignment)
##
## exterior_bisector = distance_sum / np.linalg.norm(distance_sum)
##
## delta_theta = self.alignment - exterior_bisector

#* delta_theta

else:

T_boundary = 0

eta = np.random.uniform(-1.0,1.0,3)

T_noise = eta * T_n

return T_boundary + T_noise + T_alignment

##############################################################################
########################PROPEGATION###########################################
##############################################################################

##def propegation(Penguin_list):
## for penguin in Penguin_list:
## penguin.position += penguin.net_force(F_self,F_in,k,Penguin_list,a)
## penguin.alignment += penguin.net_torque(T_in,T_n,T_align,Penguin_list,a)




##############################################################################
########################MAIN RUN##############################################
##############################################################################

# Initializing
# Placing all the penguins on a cuboid grid, with small deviations on the grid points
Penguin_list = []

for i in range(N_xy):
for j in range(N_xy):
for k in range(N_z):
Penguin_list.append(Penguin(np.random.normal(a,dev_a),
np.array([np.random.normal(i,dev_i),np.random.normal(j,dev_j),k]),
np.array([np.random.normal(0,1),np.random.normal(0,1),1]), False))


# Main RUN
x, y, z, size = np.zeros(N_xy**2*N_z, dtype=float), np.zeros(N_xy**2*N_z, dtype=float), np.zeros(N_xy**2*N_z, dtype=float), np.zeros(N_xy**2*N_z, dtype=float)
u, v, w = np.zeros(N_xy**2*N_z, dtype=float), np.zeros(N_xy**2*N_z, dtype=float), np.zeros(N_xy**2*N_z, dtype=float)
fig = plt.figure()
ax = fig.add_subplot(111, projection='3d')

#plt.show()

for l in range(t_final):
print(l)
cnt = 0
color = []
for penguin in Penguin_list:
penguin.position += penguin.net_force(F_self,F_in,k,Penguin_list,a)
penguin.alignment += penguin.net_torque(T_in,T_n,T_align,Penguin_list,a)
abs_alignment = (np.dot(penguin.alignment,penguin.alignment))**(0.5)
penguin.alignment = penguin.alignment/abs_alignment
if l > start_drawing:
x[cnt] = penguin.position[0]; y[cnt] = penguin.position[1]; z[cnt] = penguin.position[2]
u[cnt] = penguin.alignment[0]; v[cnt] = penguin.alignment[1]; w[cnt] = penguin.alignment[2]
size[cnt]= 100.*penguin.radius; #color.append(cm.spectral(1.3*(penguin.radius-0.6*a)))
color.append((1-int(penguin.boundary),int(penguin.boundary),0.))
cnt += 1
if l > start_drawing:
ax.clear()
ax.set_xlim3d(-40,40)
ax.set_ylim3d(-40,40)
ax.set_zlim3d(-40,40)
plt.ion()
ax.scatter(x,y,z,s=size,facecolors = 'none', edgecolors = color)
ax.quiver(x,y,z,u,v,w)
plt.draw()
plt.ioff()
plt.savefig('image00' + str(l) + '.png')
#pp.close()
##for penguin in Penguin_list:
## print('Position', penguin.position)
## print('Alignment', penguin.alignment)
## print('Boundary',penguin.boundary)

#fig = plt.figure()
#ax = fig.add_subplot(111, projection='3d')
#ax.scatter(x,y,z,s=size,facecolors = 'none', edgecolors = color)
#plt.draw()