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()