In [1]:

    

from astropy import constants as const
import numpy as np
import matplotlib.pyplot as plt
#This just needed for the Notebook to show plots inline. 
%matplotlib inline


    

In [22]:

    

#Number of galaxies
N = 1000
#Velocity dispersion (km/s)
sigma = 600
#Hubble constant km/s/Mpc
H = 74.0
#distance
dist=[]
vr=[]

#Draw 100 random distances out to 9Mpc distribution with probability ~ d cubed.
#This used the rejection method
i=0
while i<N:
    d = 9*np.random.random()
    f = d*d*d
    f2 = 9*9*9*np.random.random()
    if (f>f2): 
        dist.append(d)
        vr.append(H*d+np.random.normal(0,600))
        i+=1
                
n, bins, patches = plt.hist(dist, 50, normed=1, facecolor='m', alpha=0.75)


    

In [23]:

    

plt.plot(dist,vr,'.')
plt.plot([0,9],[0,666])
#Title
plt.title('Hubble plot with sigma=600 km/s')
plt.ylim(-2000,2000)
plt.xlabel("D (Mpc)")
plt.ylabel( "V_r (km/s)" )
plt.show()


    

In [21]:

    

#Number of galaxies
N = 1000
#Velocity dispersion (km/s)
sigma = 300
#Hubble constant km/s/Mpc
H = 74.0
#distance
dist=[]
vr=[]

#Draw 100 random distances out to 9Mpc distribution with probability ~ d cubed.
#This used the rejection method
i=0
while i<N:
    d = 9*np.random.random()
    f = d*d*d
    f2 = 9*9*9*np.random.random()
    if (f>f2): 
        dist.append(d)
        vr.append(H*d+np.random.normal(0,600))
        i+=1
                
plt.plot(dist,vr,'.')
plt.plot([0,9],[0,666])

#Title
plt.title('Hubble plot with sigma=300 km/s')
plt.ylim(-2000,2000)
plt.xlabel("D (Mpc)")
plt.ylabel( "V_r (km/s)" )
plt.show()


    

In [20]:

    

#Number of galaxies
N = 1000
#Velocity dispersion (km/s)
sigma = 1200
#Hubble constant km/s/Mpc
H = 74.0
#distance
dist=[]
vr=[]

#Draw 100 random distances out to 9Mpc distribution with probability ~ d cubed.
#This used the rejection method
i=0
while i<N:
    d = 9*np.random.random()
    f = d*d*d
    f2 = 9*9*9*np.random.random()
    if (f>f2): 
        dist.append(d)
        vr.append(H*d+np.random.normal(0,600))
        i+=1
                
plt.plot(dist,vr,'.')
plt.plot([0,9],[0,666])
#Title
plt.title('Hubble plot with sigma=300 km/s')
plt.ylim(-2000,2000)
plt.xlabel("D (Mpc)")
plt.ylabel( "V_r (km/s)" )
plt.show()


    

In [ ]: