import math
import numpy
import random

def densite(a,N):
    ninf = N/math.pow(a,3)
    p = int(math.floor(a/2))
    g = numpy.zeros(p)
    h = numpy.zeros(p)
    n = numpy.zeros(p)
    random.seed()
    for i in range(N):
        x = random.uniform(-a/2,a/2)
        y = random.uniform(-a/2,a/2)
        z = random.uniform(-a/2,a/2)
        r = math.sqrt(x*x+y*y+z*z)
        k = int(math.floor(r))
        if k<p:
            g[k] += 1
    for k in range(1,p):
        h[k] = h[k-1]+g[k]
    for k in range(1,p):
        n[k] = h[k]/(4.0/3*math.pi*math.pow(k,3))
    return[ninf,g,h,n]
             

a=1000.0
N = int(1e5)
[ninf,g,h,n] = densite(a,N)
             

from matplotlib.pyplot import *             

figure(figsize=(8,6))
plot(g)
xlabel("r")
ylabel("g")
grid()
             

figure(figsize=(8,6))
plot(h)
xlabel("r")
ylabel("h")
grid()
             

figure(figsize=(8,6))
plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

figure(figsize=(8,6))
for i in range(5):
    [ninf,g,h,n] = densite(a,N)
    plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

a=1000.0
N = int(1e6)
[ninf,g,h,n] = densite(a,N)
figure(figsize=(8,6))
for i in range(5):
    [ninf,g,h,n] = densite(a,N)
    plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()
             

a=1000.0
N = int(1e3)
[ninf,g,h,n] = densite(a,N)
figure(figsize=(8,6))
for i in range(5):
    [ninf,g,h,n] = densite(a,N)
    plot(n)
xlabel("r")
ylabel("n")
axis([0,500,0,ninf*5])
grid()