Skip to content

Update Penguin3_withboundary #6

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

Update Penguin3_withboundary #6

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
248 changes: 150 additions & 98 deletions Penguin3_withboundary
View file
Open in desktop
Original file line number Diff line number Diff line change
Expand Up @@ -2,152 +2,204 @@ import numpy as np
import scipy as sp
import random

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

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

N_xy = 10 # Number of penguins in the x and y direction
N_z = 25 # Number of penguins in the z direction
N_z = 12 # 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 = 1.3 * a

t_final = 3

# Physical Constants
F_self = 0.1
F_in = 0.3
k = 1.

T_in = 3.
T_n = 0.03
T_align = 0.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 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
"""

critical_radius = a
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
"""

critical_radius = a
neighbour_radius = 1.3 * a
F_selfPropulsion = F_self * self.alignment
F_selfPropulsion = F_self * self.alignment

F_repulsion = 0
F_repulsion = 0

if (self.boundary == True):
if (self.boundary == True):

F_boundary = F_in * self.alignment
F_boundary = F_in * self.alignment

else:
else:

F_boundary = 0
F_boundary = 0
above_plane = []

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 < neighbour_radius):
#Repulstion force if neighbours are too close
if (r < critical_radius):
F_repulsion += -k * r
#Boundary determination
#alignment [a b c], position [x_0, y_0, z_0], neighbour at [x, y, z]
for i in range(len(penguin_list)):
r_nb = [0,0,0]
if (penguin_list[i] != self):

r = self.get_distance(penguin_list[i])
r_mag = np.linalg.norm(r)
if (r_mag < neighbour_radius):
r_nb += r
#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])
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)
#print(above_plane)
nrofn.append(len(above_plane))
if (all(j == 0 for j in above_plane)) or (all (j == 1 for j in above_plane)):

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 < 6:
self.boundary = True
F_boundary = F_in * -r #/neigh_nr
else:
self.boundary = False
return nrofn
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
T_alignment = np.zeros(3)

if (self.boundary == True):

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):
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


distance_sum += self.get_distance(penguin_list[i])
# Torque Alignment

if (penguin_list[i] != self) and (r_mag < critical_radius):
align = np.zeros(3)

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_alignment = T_align * align

if (self.boundary == True):

## 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

T_boundary = T_in #* delta_theta

T_alignment += T_align * (self.alignment - penguin_list[i].alignment)
else:

exterior_bisector = distance_sum / np.linalg.norm(distance_sum)
T_boundary = 0

delta_theta = self.alignment - exterior_bisector
eta = np.random.uniform(-1.0,1.0,3)

T_boundary = T_in * delta_theta
T_noise = eta * T_n

else:
return T_boundary + T_noise + T_alignment

T_boundary = 0
##############################################################################
########################PROPEGATION###########################################
##############################################################################

eta = random.uniform(-1.0,1.0)
##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)

T_noise = eta * T_n

return T_boundary + T_noise + T_alignment
##############################################################################
########################INITIALIZING##########################################
########################MAIN RUN##############################################
##############################################################################

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

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))
cnt = 0
for penguin in Penguin_list:
penguin.net_force(1,1,1,Penguin_list,1)
if penguin.boundary:
cnt += 1

print(cnt)
print(nrofn)


# Main RUN
for l in range(t_final):
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)

##for penguin in Penguin_list:
## print('Position', penguin.position)
## print('Alignment', penguin.alignment)
## print('Boundary',penguin.boundary)