<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 In Class Exercise - Rutherford Scattering of a Gaussian Beam of Particles
When a positively charged particle passes close to an atom, its path will be deflected (or scatter) by an angle $\theta$ which obeys the relation
Step1:  Notice something about $b$?
Step2:  In Class Exercise
<ASSISTANT_TASK:>
Python Code:
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 
print(const.e.value)
print(const.e)
#Atomic Number of Gold
Z = 72
e = const.e.value
E = 7.7e6*e
eps0 = const.eps0.value
sigma = const.a0.value/100.
#print(Z,e,E,eps0,sigma)
N = 1000000 #Start small, and increase to 1 million when you're sure the code runs correctly.
#Function to generate two sets of random Gaussian numbers. 
def gaussian():
    r = np.sqrt(-2*sigma*sigma*np.log(1-np.random.random()))
    theta=2*np.pi*np.random.random()
    x=r*np.cos(theta)
    y=r*np.sin(theta)
    return x,y
#Main Programme
count = 0 #Initate count of particles bounced back
for i in range(N):
    x,y=gaussian()
    b=np.sqrt(x*x+y*y)
#If this is true the particle is bounced back
    if b<Z*e*e/(2*np.pi*eps0*E): 
        count +=1
        
print(count, "particles were reflected out of ", N, "incident")
print("this is a bounce fraction of {0:.5f} +/- {1:.5f}".format(count/N,np.sqrt(count)/N))
Explanation: In Class Exercise - Rutherford Scattering of a Gaussian Beam of Particles
When a positively charged particle passes close to an atom, its path will be deflected (or scatter) by an angle $\theta$ which obeys the relation: 
$$\tan (\theta/2) = \frac{Z e^2}{2 \pi \epsilon_0 E b} $$
where $Z$ is the atomic number, $e$ is the electric charge, $\epsilon_0$ is the permittivity of free space, $E$ is the kinetic energy of the incident particle, and $b$ is the impact parameter (see diagram). 
This process is called "Rutherford Scattering" after Ernest Rutherford, who was among the first physicists to explore this process. 
We will model a beam of 1 million $\alpha$ particles with a 2D Gaussian profile incident on a single atom of Gold. We would like to calculate the fraction of these particles which "bounce back" (ie. scatter through angles greater than $\theta = 90^\circ$. 
When bounce back happens, therefore $\tan (\theta/2) \gt 1$, so 
$$ b \lt \frac{Z e^2}{2 \pi \epsilon_0 E} $$
Write a programme which simulates the incident Gaussian beam and calculates the fraction of particles which bounce back. 
Please write your own function to calculate Gaussian random numbers using this format: 
def gaussian():
YOUR ALGORTHIM HERE
return x,y
Make use of these parameters for the incident beam: 
$E = $7.7 MeV
$\sigma = a_0/100$ (where $a_0$ is the Bohr radius). 
Hint: you can make use of the astropy.constants module to import various constants you need for this problem.
End of explanation
#Atomic Number of Gold
Z = 79
e = const.e.value
E = 7.7e6*e
eps0 = const.eps0.value
sigma = const.a0.value/100.
#print(Z,e,E,eps0,sigma)
N = 1000000 #Start small, and increase to 1 million when you're sure the code runs correctly.
   
#Main Programme
count = 0 #Initate count of particles bounced back
for i in range(N):
    b= np.sqrt(-2*sigma*sigma*np.log(1-np.random.random()))
#If this is true the particle is bounced back
    if b<Z*e*e/(2*np.pi*eps0*E): 
        count +=1
        
print(count, "particles were reflected out of ", N, "incident")
print("this is a bounce fraction of {0:.5f} +/- {1:.5f}".format(count/N,np.sqrt(count)/N))
?np.random.normal
#Atomic Number of Gold
Z = 79
e = const.e.value
E = 7.7e6*e
eps0 = const.eps0.value
sigma = const.a0.value/100.
print(Z,e,E,eps0,sigma)
N = 1000 #Start small, and increase to 1 million when you're sure the code runs correctly.
#Main Programme
count = 0 #Initate count of particles bounced back
for i in range(N):
    x=np.random.normal(0,sigma,1) 
    y=np.random.normal(0,sigma,1)
    b=np.sqrt(x*x+y*y)
#If this is true the particle is bounced back
    if b<Z*e*e/(2*np.pi*eps0*E): 
        count +=1
        
print(count, "particles were reflected out of ", N, "incident")
print("this is a bounce fraction of {0:.5f} +/- {1:.5f}".format(count/N,np.sqrt(count)/N))
Explanation: Notice something about $b$?
End of explanation
#Define the function
def f(x):
    fx = (np.sin(1/(x*(2-x))))**2
    return fx
#Integrate the function from x=0-2
#Note that you need to know the maximum value of the function
#over this range (which is y=1), and therefore the area of the box
#from which we draw random number is A=2. 
N=1000000
k=0
for i in range(N):
    x=2*np.random.random()
    y=np.random.random()
    if y<f(x):
        k+=1
A=2.        
I=A*k/N
print("The integral is equal to I = ",I)
Explanation: In Class Exercise: Monte Carlo Integration
Write a programme following the method of Monte Carlo Integration to calculate 
$$ I = \int_0^2 \sin^2 [\frac{1}{x(2-x)}] dx. $$
As you will need to calculate $f(x) = \sin^2 [\frac{1}{x(2-x)}]$ many times please write a user defined function for this part of your programme.
End of explanation
<END_TASK>