from visual import * from visual.graph import * from random import random entropygraph = gdisplay(x=400, y=400, xtitle='Time', ytitle='Response (click and drag mouse to see coordinates)') entropyplot = gcurve(color=color.green) NSites=20 NOccupied=NSites/2 leftprob=0.25 rightprob=1-leftprob Sites = [] statearray=[] checkarray=[] leftsum=0 rightsum=0 # Set up the initial array with all occupied sites on the left edge for i in range(NSites): if i<=(NOccupied-1): initstate=1 initcolor=color.red else: initstate=0 initcolor=color.white if i<(NSites/2): leftsum=leftsum+initstate else: rightsum=rightsum+initstate print(i,leftsum) Sites=Sites+[sphere(radius=0.01*NSites, pos=vector((-NSites/2)+i,0,0), color=initcolor, state=initstate)] statearray=statearray+[initstate] checkarray=checkarray+[random()] n=0 #calculate the initial entropy leftentropy=log(factorial(NSites/2)/(factorial((NSites/2)-leftsum)*factorial(leftsum))) rightentropy=log(factorial(NSites/2)/(factorial((NSites/2)-rightsum)*factorial(rightsum))) entropyplot.plot( pos=(n, leftentropy+rightentropy) ) #print(statearray) print(n,leftsum,rightsum,leftentropy+rightentropy) while n<1000: rate(100) for i in range(1,NSites): #print(i) if Sites[i].state==0 and Sites[i-1].state==1: if checkarray[i-1]>rightprob: # If current site is empty and site to left is full # have "atom" jump right with some probability Sites[i].state=1 Sites[i].color=color.red Sites[i-1].state=0 Sites[i-1].color=color.white else: if Sites[i].state==1 and Sites[i-1].state==0: if checkarray[i]