anasqiblawi

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)))
  47.     i2 = np.where((r > separation+force_radius/2) & (r < separation+force_radius))
  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
  115. separation_radius = 10 #mean separation radius
  116. interaction_radius = 25 # mean interaction radius
  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
  125. force_radius = block_matrix(interaction_radius,interaction_radius*deviation,True)
  126. separation = block_matrix(separation_radius,separation_radius*deviation,True)
  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.  
  137.     fr = np.abs(force_radius[0,g*ppg])/np.max(force_radius)
  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')
  182. plt.tight_layout(pad=0.3)
  183. plt.subplots_adjust(hspace=0.01, wspace=0.01)
  184. # figManager = plt.get_current_fig_manager()
  185. # figManager.window.showMaximized()
  186. plt.show()
RAW Paste Data

Adblocker detected! Please consider disabling it...

We've detected AdBlock Plus or some other adblocking software preventing Pastebin.com from fully loading.

We don't have any obnoxious sound, or popup ads, we actively block these annoying types of ads!

Please add Pastebin.com to your ad blocker whitelist or disable your adblocking software.

×