#!/usr/bin/env python3 # -*- coding: utf-8 -*- """ Animation of Maxwell's Demon author: Benjamin Preis email: benpreis@mit.edu ----------------------- Modified from Jake Vanderplas's Matplotlib Animation Tutorial of particles in a box website: http://jakevdp.github.com ----------------------- """ import numpy as np import matplotlib.pyplot as plt import matplotlib.animation as animation class Box: def __init__(self, #The overall bounds bounds=[0,8,0,4], #The immovable part bottom_line=[4,4,0,3], #The door that Maxwell's demon opens top_line=[4,4,4,4]): self.bounds=bounds self.bottom_line=bottom_line self.top_line=top_line self.time_elapsed = 0 def step(self, dt): self.time_elapsed +=dt #This is some poorly written code, as it's hugely self-referential #And not easily changable #This is the product of me brute-force teaching myself Python this week #Anyway #Evaluate as true if ... the fast balls (as defined below) are #To the right of the divider demon_r=np.any((balls_right.state[:, 0] > 4) & #But only somewhat to the right (balls_right.state[:, 0] <4.08) & #And they're within the gate area (balls_right.state[:,1] > self.bottom_line[3] + balls_right.size)) #Evaluate as true if ... the slow balls (as defined below) are #To the left of the divider demon_l=np.any((balls_left.state[:, 0] < 4) & #But only somewhat to the left (balls_left.state[:, 0] >3.92) & #And within the gate area (balls_left.state[:,1] > self.bottom_line[3] + balls_left.size)) #If either of the above conditions are true, raise the gate if (demon_r or demon_l): self.top_line[2] = 4 #Otherwise, lower it else: self.top_line[2] = 3 class Particles: """ init_state is an [N x 4] array, where N is the number of particles: [[x1, y1, vx1, vy1], [x2, y2, vx2, vy2], ... ] bounds is the size of the box: [xmin, xmax, ymin, ymax] """ def __init__(self, init_state = [], size = 0, M = 0.05): self.init_state = np.asarray(init_state, dtype=float) self.M = M * np.ones(self.init_state.shape[0]) self.size = size self.state = self.init_state.copy() self.time_elapsed = 0 def step(self, dt): """step once by dt seconds""" self.time_elapsed += dt # update positions self.state[:, :2] += dt * self.state[:, 2:] # check for crossing edge #To the left crossed_x1 = (self.state[:, 0] < box.bounds[0] + self.size) #To the right crossed_x2 = (self.state[:, 0] > box.bounds[1] - self.size) #Below crossed_y1 = (self.state[:, 1] < box.bounds[2] + self.size) #Above crossed_y2 = (self.state[:, 1] > box.bounds[3] - self.size) #Among those particles that crossed an edge, push them back into place self.state[crossed_x1, 0] = box.bounds[0] + self.size self.state[crossed_x2, 0] = box.bounds[1] - self.size self.state[crossed_y1, 1] = box.bounds[2] + self.size self.state[crossed_y2, 1] = box.bounds[3] - self.size #If they crossed an edge, reverse their velocity in that direction self.state[crossed_x1 | crossed_x2, 2] *= -1 self.state[crossed_y1 | crossed_y2, 3] *= -1 #Check for crossing bottom_line, the permanent part #If the particle is ... to the left of the middle middle_x1 = ((self.state[:, 0] < box.bottom_line[0] + self.size) & #And if it was previously to the right (self.state[:,0]+self.size>box.bottom_line[0]) & #And it's moving to the left (self.state[:,2]<0)& #And it's below the opening (self.state[:, 1] < box.bottom_line[3] + self.size)) #If the particle is ... to the right of the middle middle_x2 = ((self.state[:, 0] > box.bottom_line[1] - self.size) & #And was previously to the left (self.state[:,0]-self.size0) & #And is below the opening (self.state[:, 1] < box.bottom_line[3] + self.size)) #Evaluate the above, but with the added condition of velocity (see np.square) middle_demon1 = ((self.state[:, 0] < box.bottom_line[0] + self.size) & #And if it was previously to the right (self.state[:,0]+self.size>box.bottom_line[0]) & #And it's moving to the left (self.state[:,2]<0)& #And it's above the opening (self.state[:, 1] > box.bottom_line[3] + self.size)& #And if it is slow, keep it to the right #32 derives from the initial conditions, given below ((np.square(self.state[:,2])+np.square(self.state[:,3]))<32)) middle_demon2 = ((self.state[:, 0] > box.bottom_line[1] - self.size) & #And was previously to the left (self.state[:,0]-self.size0) & #And is above the opening (self.state[:, 1] > box.bottom_line[3] + self.size) & #And if it is fast, keep it to the left ((np.square(self.state[:,2])+np.square(self.state[:,3]))>32)) #Revert back to its previous side if below the bottom line #Or within the gate area, but on it desired side self.state[middle_x1, 0] = box.bottom_line[0] + self.size self.state[middle_x2, 0] = box.bottom_line[1] - self.size self.state[middle_demon1,0] = box.bottom_line[1] + self.size self.state[middle_demon2,0] = box.bottom_line[1] - self.size #self.state[middle_y1, 1] = box.bottom_line[3] + self.size #If it needs to stay on its side, reverse its x velocity self.state[middle_x1|middle_x2|middle_demon1|middle_demon2, 2] *= -1 #------------------------------------------------------------ # set up initial state np.random.seed(0) init_state1 = np.random.random((50, 4)) init_state1[:, :2] *= 8 init_state1[:,2:]*=4 init_state2 = np.random.random((50,4)) init_state2[:,:2] *= 8 init_state2[:,2:] = 4*init_state2[:,2:]+4 box=Box() balls_left = Particles(init_state1,size=0.04) balls_right = Particles(init_state2,size=0.04) dt = 1./300 # 30fps #------------------------------------------------------------ # set up figure and animation fig = plt.figure() fig.subplots_adjust(left=0, right=1, bottom=0, top=1) ax = fig.add_subplot(111, aspect='equal', autoscale_on=False, xlim=(-.5, 8.5), ylim=(-.5, 4.5)) # particles holds the locations of the particles particlesl, = ax.plot([], [], 'bo', ms=6) particlesr, = ax.plot([],[],'ro',ms=6) # rect is the box edge rect = plt.Rectangle(box.bounds[::2], box.bounds[1] - box.bounds[0], box.bounds[3] - box.bounds[2], ec='none', lw=2, fc='none') ax.add_patch(rect) #Generate the bottom line bottom_line = plt.Line2D((box.bottom_line[0], box.bottom_line[1]), (box.bottom_line[2], box.bottom_line[3]), lw=2.5) ax.add_line(bottom_line) #Create a placeholder for the gate top_line = plt.Line2D([],[], lw=2.5) ax.add_line(top_line) def init(): """initialize animation""" global balls_left, balls_right, box, rect particlesl.set_data([], []) particlesr.set_data([],[]) rect.set_edgecolor('none') top_line.set_data([],[]) return particlesl, particlesr, top_line, rect def animate(i): """perform animation step""" global box, balls_left, balls_right, rect, dt, ax, fig balls_left.step(dt) balls_right.step(dt) box.step(dt) ms = int(fig.dpi * 2 * balls_left.size * fig.get_figwidth() / np.diff(ax.get_xbound())[0]) # update pieces of the animation rect.set_edgecolor('k') particlesl.set_data(balls_left.state[:, 0], balls_left.state[:, 1]) particlesl.set_markersize(ms) particlesr.set_data(balls_right.state[:, 0], balls_right.state[:, 1]) particlesr.set_markersize(ms) top_line.set_data((box.top_line[0],box.top_line[1]),(box.top_line[2],box.top_line[3])) return particlesl, particlesr, top_line, rect ani = animation.FuncAnimation(fig, animate, frames=600, interval=10, blit=True, init_func=init) plt.show()