import numpy as np import random as R from scipy.misc.common import factorial # Structural explanation # A short summary of what the classes are meant to do: # # - Environment: # This class implements the drift velocity field and the boundary # conditions for the tea cup in part 2 of the homework. # # - Walker: # A simple random walker in dim dimensions which does random steps along # the axes of the coordinate system and obeys the drift and boundary # conditions of the environment. Note that each walker carries its # own position, thereby creating the square grid in part 1. The # class Grid in fact implements a grid of cells. # # - Grid: # A square grid of cells in two dimensions, used to count particles # within the cells in order to calculate the entropy of the system, # based on the number of micro-states. # # - Group: # A group of N random walkers in an environment. All walkers in a group # are initialised at the same position. The group can tell all its # walkers to do a step, and can calculate the entropy of the system. # def logfact(x): if x <= 120: return np.log(factorial(x)) else: return 0.5*np.log(2.0*np.pi*x) + x*np.log(x) - x class Environment: def __init__(self, limited=False, drifting=False, vr=0.02, vphi=10): # Initialise an environment: # - The boundaries are, by construction, circular, with radius # set, by default, on self.radius = 42. The boundaries become # relevant only, if limited is set to "True". # - The drift can be switched on and off by the bool "drifting". # The drift field has two components, a radial one with constant # velocity vr, and a tangential one with a profile given by # Question (2.b). The two drifts should be implemented in # r_profile and phi_profile, respectively. self.limited = limited self.radius = 42. self.drifting = drifting self.vphi = vphi self.vphi_rmax = 20. self.vr = vr # which number of dimensions can we cope with? self.accepts_dim = [2] def r_profile(self,r): # Fill this function, Q(2.b) return 0. def phi_profile(self,r): # Fill this function, Q(2.b) return 0. def drift(self,pos): # Complete this function, use r_profile and phi_profile. if not self.drifting: return 0. # Relevant stuff should go here. return np.array([0.,0.]) def impose_limit(self,pos): # Fill this function, Q(2.a) pass class Walker: def __init__(self, dim, env, pos=None): # - Initialise a random walker in dim dimensions, at position pos # or the origin; # - Fix the stepsize to be 0.5; # - Set the environment variable if pos is None: pos = np.zeros(dim) self.pos = pos.copy() self.Ndim = dim self.stepsize = .5 self.env = env if env is not None and dim not in env.accepts_dim: print "Environment cannot cope with %s dimensions." % dim, print "Acceptable are %s." % env.accepts_dim print "Environment disabled." self.env = None def step(self): # This function is used in parts 1 and 2 of the homework. # # For the environment in part 1, you will have no drift, and # no boundaries. See the "Environment"-class for more detail. # # In part 2, you will have an environment, with both a drift # switched on, and a check for boundaries. Again, check the # documentation of "Environment". self.randStep() if self.env is not None: self.pos += self.env.drift(self.pos) self.env.impose_limit(self.pos) def randStep(self): # Fill this function, Q(1.a) # Implement a step of size self.stepsize in dim dimensions # # You may use the python function randint(N,M), which # yields an integer random number between N and M # You can find this function in module "random". pass class Grid: def __init__(self,d=5): # Initialise a square grid of infinite size with cell size self.d, # used to calculate the number of micro-states for determining the # entropy. # - self.totalN is the total number of particles in the grid; # - self.cells is a dictionary having the cell coordinates as keys # and the number of particles in the cell as value. self.d = d self.totalN = 0 self.cells = {} def addPoint(self,pt): # Fill this function, Q(1.b) # This function expects the position of a walker in 2 dimensions as # argument. It increases the counter for the total number of # particles in the grid, as well as the counter of the cell the # particle is in. Use the walker's coordinates and the grid # spacing self.d to calculate the cell coordinates which you can # then use as key for the self.cells dictionary. pass def counts(self): # Return a list of the number of particles in each cell that has at # least one particle. return self.cells.values() def total(self): # Return the total number of particles return self.totalN class Group: def __init__(self, N, env, pos=np.zeros(2), style='o'): # Initialise a group of N 2D walkers at position pos with an # environment env. self.style is the plotting style for the # walkers. self.env = env self.walkers = [ Walker(2,env,pos) for i in range(N) ] self.style = style def step(self): # Let each walker in the group make a step, and, depending on the # environment, experience drift and boundaries. for w in self.walkers: w.step() def pos(self): # Return an array with the positions of all the walkers in the # group. This is used for plotting. return np.array([ w.pos for w in self.walkers ]) def xs(self): # Return an array of x-positions of all the walkers in the group. return self.pos()[:,0] def ys(self): # Return an array of y-positions of all the walkers in the group. return self.pos()[:,1] def entropy(self): # Fill this function, Q(1.b) # # Initialise and use the "Grid" class to calculate the total # number of particles and the number of micro-states. # Then use the formulae given in the lecture to calculate # S=k_B*ln(Omega) and return S. You may want to use the function # logfact provided to you. return 100. ##### show histo in x-projection! if __name__ == "__main__": import pylab as pl choice = None while choice not in ['e','m']: choice = raw_input("Type 'e' for Entropy or 'm' for Milk: ") entropy = ( choice == 'e' ) if entropy: env = None groups = [ Group(200, env, np.array([0.,0.]), 'b.') ] steps = 5000 else: drift, limits = None, None while drift not in ['y','n']: drift = raw_input("Drift field (y/n): ") while limits not in ['y','n']: limits = raw_input("Limits (y/n): ") env = Environment( drifting = (drift == 'y'), limited = (limits == 'y') ) groups = [ Group(50, env, np.array([15.,0.]), 'bo'), Group(50, env, np.array([25.,0.]), 'ro'), Group(50, env, np.array([35.,0.]), 'yo'), Group(50, env, np.array([0.,0.]), 'ko') ] steps = 500 pl.ion() if entropy: scatter = pl.subplot('211') entro = pl.subplot('212') else: scatter = pl.subplot('111') plots = {} for g in groups: plots[g] = scatter.plot(g.xs(), g.ys(), g.style)[0] eplots = {} if entropy: ent_xs = [0.] ent_ys = {} for g in groups: ent_ys[g] = [g.entropy()] # print ent_ys.items() eplots[g] = entro.plot(ent_xs,ent_ys[g],'r-')[0] #################### ### Draw a circle #################### if env is not None and env.limited: circ = np.linspace(0,2*np.pi,100) scatter.plot( 42.*np.cos(circ), 42.*np.sin(circ), 'k' ) ######################## ### Draw velocity field ######################## if env is not None and env.drifting: xs = np.linspace(-70,70,50) ys = np.zeros(len(xs)) drifts = [ env.drift(np.array([x,0])) for x in xs ] vx, vy = zip(*drifts) scatter.quiver(xs,ys,vx,vy,width=2e-3,color='grey') scatter.set_xlim(-70,70) scatter.set_ylim(-50,50) if entropy: entro.set_xlim(0,steps) entro.set_ylim(0,logfact(200)) entro.set_ylabel('Entropy [k_B]') entro.set_xlabel('Iteration steps') scatter.set_aspect('equal') if entropy: update = 10 # only plot every nth iteration else: update = 1 for i in range(steps): for g in groups: g.step() if i % update == 0: plots[g].set_data(g.xs(), g.ys()) if entropy: ent_ys[g].append(g.entropy()) ent_xs.append(float(i)) eplots[g].set_data(ent_xs,ent_ys[g]) if i % update == 0: scatter.set_title( '%s of %s steps' % (i,steps) ) pl.draw() scatter.set_title( '%s steps' % steps ) pl.draw() pl.ioff() pl.show()