gistfile1.txt

import random
from math import *

MAXSTEPS = 300
NUM_DRUNKARDS = 20000
J = 12
JSQR = J ** 2
DIRECTIONS = 1

random.seed()

class RandomWalk:

    def __init__(self):
        self.directions = [0]
        

    def walk(self):
        """
        Move the particle by one step
        """
        self.directions[0] += random.choice((1, -1))    


    def arrested(self):
        """
        Calculates if the particle (drunkard) has reached the edge of the
        potential well and needs to be stopped (arrested)
        """

        distance = 0
        for i in self.directions:
            distance += i ** 2

        return distance == JSQR
             

def filter(arr):
    """
    Do some post-processing of the data
    """
    out = []
    out2 = []
    for i, a in enumerate(arr):
        if i % 2 == 0:
            # Remove the zero values
            out.append(log(a) if a else 1)
            out2.append(i)
    return [out, out2]


def main():

    bins = [0] * MAXSTEPS
    
    for j in xrange(NUM_DRUNKARDS):
        W = RandomWalk()
        for i in xrange(MAXSTEPS):
            if W.arrested():
                bins[i] += 1
                break
            W.walk()
        if j % 1000 == 0:
            print 'Done 1000 drunkards'

    # The next line is a hack due to a mistake by my supervisor...
    survivors = bins #W.numberSurvived(bins)
    yerr = []    
    x, y = filter(survivors)
    
    for i in y:
        if survivors[i] != 0:
            yerr.append(1.0/sqrt(survivors[i]))
        else:
            yerr.append(0)

    fp = open('1D.dat', 'w')
    for i in range(len(x)):
        fp.write(str(y[i]) + " " + str(x[i]) + " " + str(yerr[i]) + "\n")
    fp.close()


if __name__ == "__main__":
    # Import Psyco if available
    import datetime
    now = datetime.datetime.now()

    try:
        import psyco
        psyco.full()
    except ImportError:
        print 'Cannot import psyco\n'

    main()

    print datetime.datetime.now() - now