# Particle Game of Life

Apr 15th, 2021
475
Never
Not a member of Pastebin yet? Sign Up, it unlocks many cool features!
1. import numpy as np
2. import matplotlib.pyplot as plt
3. import matplotlib.animation as anim
4. import colorsys as colsys
5.
6. np.set_printoptions(precision=3)
7. plt.rcParams['axes.facecolor'] = '#000000'
8.
9. def block_matrix(p0,p_std,pos=False):
10.     """
11.     generates block matrix (used for assigning
12.     parameters to groups of particles)
13.     """
14.     m = np.random.normal(p0,p_std,(groups,groups))
15.     m = m.repeat(particles_per_group,axis=1).repeat(particles_per_group,axis=0)
16.
17.     if pos:
18.         return np.abs(m)
19.     else:
20.         return m
21.
22. def bound(x,y):
23.     """
24.     periodic boundaries
25.     """
26.     x[x>=x_width*lim] -= 2*x_width*lim
27.     x[x<-x_width*lim] += 2*x_width*lim
28.     y[y>=lim] -= 2*lim
29.     y[y<-lim] += 2*lim
30.
31.     return x,y
32.
33.
34. def attract(x,y,vx,vy):
35.     """
36.     calculate forces between particles
37.     """
38.
39.     dx = x.reshape(total_particles,1) - x
40.     dy = y.reshape(total_particles,1) - y
41.     dx,dy = bound(dx,dy)
42.
43.     r = np.sqrt(dx*dx + dy*dy)
44.     f = np.zeros((total_particles,total_particles))
45.     i0 = np.where(r < separation)
46.     i1 = np.where((r > separation) & (r < (separation+force_radius/2)))
48.
49.     f[i0] = rep_force[i0]*(separation[i0]-r[i0])
50.     f[i1] = forces[i1] * (r[i1]-separation[i1])
51.     f[i2] = -forces[i2] * (r[i2]-(separation[i2]+force_radius[i2]))
52.
53.     np.fill_diagonal(f,0) # no self-force
54.
55.     # increase force at long range
56.     f *= close_range_factor + (dist_range_factor-close_range_factor) * (r/lim)
57.
58.     a = np.arctan2(dy,dx)
59.
60.     # add contributions from all particles
61.     fx = np.sum(f * np.cos(a), axis=0)
62.     fy = np.sum(f * np.sin(a), axis=0)
63.
64.     return fx,fy
65.
66. def friction(vx,vy):
67.     """
68.     friction and thermal forces
69.     """
70.
71.     v = np.sqrt(vx*vx + vy*vy)
72.     a = np.arctan2(vy,vx)
73.     fx = friction_coefficient * v * np.cos(a) + np.random.normal(0,temperature,total_particles)
74.     fy = friction_coefficient * v * np.sin(a) + np.random.normal(0,temperature,total_particles)
75.
76.     return fx,fy
77.
78. def init():
79.     particles.set_offsets([])
80.     particles_area.set_offsets([])
81.     return particles, particles_area,
82.
83. def animate(i):
84.
85.     global x,y,vx,vy
86.
87.     x,y = bound(x,y)
88.     fx,fy = attract(x,y,vx,vy)
89.     ffx,ffy = friction(vx,vy)
90.
91.     x += vx * dt
92.     y += vy * dt
93.     vx += (-fx-ffx) * dt
94.     vy += (-fy-ffy) * dt
95.
96.     xy = np.vstack((x,y)).T
97.
98.     particles.set_offsets(xy)
99.     particles_area.set_offsets(xy)
100.
101.     return particles, particles_area,
102.
103. # global settings
104. dt = 0.01
105. lim = 60
106. x_width = 1
107. groups = 3
108. particles_per_group = 200
109. total_particles = groups * particles_per_group
110.
111. # parameters
112. repelling_force = 60 # force active if r < separation radius
113. temperature = 5 # controls random fluctuations in particle velocities
114. friction_coefficient = 90
117. force_strength = 10 # inter-particle force strength
118. close_range_factor = 1 # force strength multiplier at r=0
119. dist_range_factor = 5 # force strength multiplier at r=lim
120. deviation = 0.1 # spread in group parameters
121. seed_range = 0.9 # initial position spread
122.
123.
124. # parameter matrices
127. forces = block_matrix(0,force_strength)
128. rep_force = block_matrix(repelling_force,repelling_force*deviation,True)
129.
130.
131.
132. # group properties
133. cs,ss,sa = [],[],[]
134. ppg = particles_per_group
135. for g in range(groups):
136.
138.     fs = np.abs(forces[0,g*ppg])/np.max(np.abs(forces))
139.     sep = np.abs(separation[0,g*ppg])/np.max(separation)
140.     rep = np.abs(rep_force[0,g*ppg])/np.max(rep_force)
141.     fd = np.abs(forces[g*ppg,g*ppg])/np.max(np.abs(forces))
142.
143.     size = 20 + 50 * fs
144.
145.     fc = (forces[0,g*ppg]-np.min(forces))/(np.max(forces)-np.min(forces))
146.
147.     color = colsys.hsv_to_rgb(fc,0.3+0.7*fd,0.5*(1+fs))
148.     for p in range(particles_per_group):
149.         cs.append(color)
150.         ss.append(size)
151.         sa.append(2*fr*15*size)
152.
153. # initialization
154. x = np.random.uniform(-lim*seed_range,lim*seed_range,total_particles)
155. y = np.random.uniform(-lim*seed_range,lim*seed_range,total_particles)
156. vx = np.random.uniform(-1,1,total_particles)
157. vy = np.random.uniform(-1,1,total_particles)
158.
159. # parameter settings
160. # titlestring = ''
161. # titlestring += f'dt {dt:.2f} | '
162. # titlestring += f'Rep. F {repelling_force:.1f} | '
163. # titlestring += f'Fric. {friction_coefficient:.1f} | '
164. # titlestring += f'Sep. {separation_radius:.1f} | '
165. # titlestring += f'Int. R {interaction_radius:.1f} | '
166. # titlestring += f'F Str. {force_strength:.1f} | '
167. # titlestring += f'Temp. {temperature:.1f} | '
168.
169.
170.
171. fig,ax = plt.subplots(figsize=(14,7),facecolor='k')
172. # plt.title(titlestring)
173. ax.set_aspect('equal')
174. plt.xlim(-x_width*lim,x_width*lim)
175. plt.ylim(-lim,lim)
176. particles = plt.scatter(x, y, marker='.', c=cs, s=ss)
177. particles_area = plt.scatter(x, y, marker='.', c=cs, s=sa, alpha=0.1, edgecolors=None, linewidths=0)
178.
179. animation = anim.FuncAnimation(fig, animate, init_func=init, interval=0, blit=True)
180.
181. plt.axis('off')
184. # figManager = plt.get_current_fig_manager()
185. # figManager.window.showMaximized()
186. plt.show()
RAW Paste Data