PenguinMotion/Initializing at master · meesteph/PenguinMotion · GitHub

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
import scipy as sp
import numpy as np

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

N_xy = 2        # Number of penguins in the x and  y direction
N_z = 3         # Number of penguins in the z direction

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

def heaviside(number):

	if number < 0:
		return 0
	else:
		return 1

# def find_theta(vector1, vector2):

# 	v1_mag = np.linalg.norm(vector1)
# 	v2_mag = np.linalg.norm(vector2)

# 	return np.arccos(np.dot(vector1, vector2) / (v1_mag * v2_mag))

class Penguin(object):

	def __init__(self, radius, position, alignment):

		self.radius = radius
		self.position = position	# Use numpy array for position
		self.alignment = alignment

	def __str__(self):

		return "Penguin located at {}".format(self.position)

	def __repr__(self):

		return self.__str__()

	def get_distance(self, penguin2):

		return self.position - penguin2.position

	def net_force(self, F_self, F_in, k, penguin_list):

		crit_radius = 1.0

		F_sp = F_self * self.alignment

		theta_list = []
		net_repulse_distance = np.array([0,0,0])

		F_boundary = np.array([0,0,0])

		for i in range(len(penguin_list)):

			if (penguin_list[i] != self):

				r = self.get_distance(penguin_list[i])
				net_repulse_distance += r

				if np.linalg.norm(r) < crit_radius:

					theta = np.arctan2(r[1], r[0])
					theta_list.append(theta)

		theta_list.sort()

		for j in range(len(theta_list)):

			try:

				difference = theta_list[j+1] - theta_list[j]

			except IndexError:

				difference = (np.pi - theta_list[j]) + (theta_list[0] + np.pi)

			print(difference * 180 / np.pi)
			F_boundary += (difference - np.pi) * F_in * heaviside(difference - np.pi) * self.alignment

			print(F_boundary)

		F_repulsion = -k * net_repulse_distance

		return F_sp + F_boundary + F_repulsion

##particle1 = Penguin(0.4, np.array([0,0,0]), np.array([1,1,1]))
##particle2 = Penguin(0.4, np.array([0.25,0.25,0]), np.array([1,1,1]))
##particle3 = Penguin(0.4, np.array([-0.25,0.25,0]), np.array([1,1,1]))
##
##penguin_list = [particle1, particle2, particle3]
##
##print(particle1.net_force(1,1,0,penguin_list))

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