One D spheres.py

From Werner KRAUTH

Here is the program discussed in Lecture 05, for computing the probability to have a particle at position $x$

import pylab

def binomialCoeff(n, k):
    result = 1
    for i in range(1, k+1):
        result = result * (n-i+1) / i
    return result

def Z(N, L, sigma):
    freespace = L - 2.0 * N * sigma
    if freespace > 0.0:
        result = freespace ** N
    else:
        result = 0.0
    return result

def pi(x, N, L, sigma):
    tot = 0.0
    for k in range(0, N):
        Z1 = Z(k, x - sigma, sigma)
        Z2 = Z(N - k - 1, L - x - sigma, sigma)
        tot += binomialCoeff( N - 1, k) * Z1 * Z2
    Ztotal = Z(N, L, sigma)
    return tot / Ztotal

L = 105.0
N = 100
sigma = .50
xr = pylab.linspace(0.0, L, 2001)
yr = [pi(x, N, L, sigma)  for x in xr]
pylab.plot(xr, yr, 'red', linewidth=2.0)
pylab.title('One-dimensional hard-sphere density, $L$ = ' + str(L) + ' $N$ = '
+ str(N) )
pylab.xlabel('$x$ (position)')
pylab.ylabel('$\pi(x)$ (probability for a particle to be at $x$)' )
pylab.show()