RandWalk.py

import random
import numpy as np
import matplotlib.pyplot as plt
import matplotlib.mlab as mlab
from math import sqrt

def magnitude(x, y, z):
    return sqrt(x**2 + y**2 + z**2)

for k in xrange(1,3):
    Rx, Ry, Rz = 0, 0, 0
    RxSqrS, RySqrS, RzSqrS = 0, 0, 0
    RxT, RyT, RzT = 0, 0, 0
    RSqrHist = []
    
    for n in xrange(100):
        RxS, RyS, RzS = 0, 0, 0
        
        for N in xrange(10* 10**(k - 1)):
            x = random.uniform(-1, 1)
            y = random.uniform(-1, 1)
            z = random.uniform(-1, 1)
            mag = magnitude(x, y, z)
            ux, uy, uz = x / mag, y / mag, z / mag
            
            # Sanity check
            assert round(magnitude(ux, uy, uz), 5) == 1.0
            
            # Add unit vectors to R#S
            RxS += ux
            RyS += uy
            RzS += uz
            
            # Calculate R#SqrS
            RxSqrS = RxS**2
            RySqrS = RyS**2
            RzSqrS = RzS**2
            
            RSqrHist.append(magnitude(RxSqrS, RySqrS, RzSqrS))
        
        RxT += RxS
        RyT += RyS
        RzT += RzS

# Plot the histogram(s)        
fig = plt.figure()
ax = fig.add_subplot(111)
n, bins, patches = ax.hist(RSqrHist, 50)
ax.autoscale()
ax.grid(True)

plt.show()