Add code comments
This commit is contained in:
parent
ed96bc6c0d
commit
02e81cf7bb
GPG key ID:
9814758436430045
2 changed files with 50 additions and 57 deletions
|
|
@ -1,9 +1,9 @@
|
|||
import matplotlib.pyplot as plt
|
||||
import numpy as np
|
||||
from scipy import integrate
|
||||
from tqdm import tqdm
|
||||
from tqdm import tqdm #Import all needed modules
|
||||
|
||||
def Fresnel2dreal(yp, xp, y, x, k, z):
|
||||
def Fresnel2dreal(yp, xp, y, x, k, z): #Define functions for integral kernels
|
||||
kernel = np.cos((k/(2*z))*((x-xp)**2+(y-yp)**2))
|
||||
return kernel
|
||||
|
||||
|
|
@ -11,85 +11,80 @@ def Fresnel2dimag(yp, xp, y, x, k, z):
|
|||
kernel = np.sin((k/(2*z))*((x-xp)**2+(y-yp)**2))
|
||||
return kernel
|
||||
|
||||
c = 3e8
|
||||
c = 3e8 #Universal constants
|
||||
e0 = 8.85e-12
|
||||
|
||||
def plot1D(aperture, z, k, screen_range, resolution):
|
||||
def genData(aperture, z, k, screen_range, resolution):
|
||||
y = 0
|
||||
def plot1D(aperture, z, k, screen_range, resolution): #Function for part 1
|
||||
def genData(aperture, z, k, screen_range, resolution): #Function to generate data
|
||||
y = 0 #As in 1D
|
||||
|
||||
xp1=yp1=-aperture/2
|
||||
xp1=yp1=-aperture/2 #Set range over which we are integrating
|
||||
xp2=yp2=aperture/2
|
||||
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to calculate for
|
||||
intensities = []
|
||||
|
||||
constant = k/2*np.pi*z
|
||||
completion = 0
|
||||
constant = k/2*np.pi*z #Relative intensity constant
|
||||
|
||||
for x in tqdm(xs):
|
||||
realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)
|
||||
for x in tqdm(xs): #Loop through all values of x
|
||||
realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10) #Calculate both real and imaginary parts
|
||||
imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)
|
||||
|
||||
I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)
|
||||
intensities.append(I)
|
||||
completion = completion + 100/resolution
|
||||
print(completion)
|
||||
I = c*e0*((realpart*constant)**2+(imagpart*constant)**2) #Combine parts and constants to get intensity
|
||||
intensities.append(I) #Add calculated intensity to list
|
||||
return xs, intensities
|
||||
|
||||
ax = plt.axes()
|
||||
xs, intensities = genData(aperture, z, k, screen_range, resolution)
|
||||
xs, intensities = genData(aperture, z, k, screen_range, resolution) #Plot intensity against distance
|
||||
ax.plot(xs, intensities)
|
||||
# xs, intensities = genData(2e-5, 0.05, 8.377e6, 0.015)
|
||||
# ax.plot(xs, intensities)
|
||||
plt.show()
|
||||
|
||||
def plot2Drectangular(aperture, z, k, screen_range, resolution):
|
||||
def genData(aperture, z, k, screen_range, resolution):
|
||||
xp1=yp1=-aperture/2
|
||||
def plot2Drectangular(aperture, z, k, screen_range, resolution): #Function for part 2
|
||||
def genData(aperture, z, k, screen_range, resolution): #Function to generate data
|
||||
xp1=yp1=-aperture/2 #Set integration limits
|
||||
xp2=yp2=aperture/2
|
||||
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to integrate for
|
||||
ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
|
||||
intensities = []
|
||||
|
||||
constant = k/2*np.pi*z
|
||||
constant = k/2*np.pi*z #Relative intensity constant
|
||||
completion = 0
|
||||
|
||||
for y in tqdm(ys):
|
||||
xIntensities = []
|
||||
for x in xs:
|
||||
realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)
|
||||
realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)#Calculate both parts
|
||||
imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1, yp2, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)
|
||||
|
||||
I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)
|
||||
I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)#Combine both parts and constants
|
||||
xIntensities.append(I)
|
||||
completion = completion + 100/resolution**2
|
||||
print(completion)
|
||||
intensities.append(xIntensities)
|
||||
intensities = np.array(intensities)
|
||||
return intensities
|
||||
|
||||
intensity = genData(aperture, z, k, screen_range, resolution)
|
||||
extents = (-screen_range/2,screen_range/2,-screen_range/2,screen_range/2)
|
||||
intensity = genData(aperture, z, k, screen_range, resolution) #Generate the data
|
||||
extents = (-screen_range/2,screen_range/2,-screen_range/2,screen_range/2) #Set limits of the plot
|
||||
|
||||
plt.imshow(intensity,vmin=0.0,vmax=1.0*intensity.max(),extent=extents,origin="lower",cmap="nipy_spectral_r")
|
||||
plt.imshow(intensity,vmin=0.0,vmax=1.0*intensity.max(),extent=extents,origin="lower",cmap="nipy_spectral_r") #Plot the graphs
|
||||
plt.colorbar()
|
||||
plt.show()
|
||||
|
||||
def plot2Dcircular(aperture, z, k, screen_range, resolution):
|
||||
def plot2Dcircular(aperture, z, k, screen_range, resolution): #Function for part 3
|
||||
def genData(aperture, z, k, screen_range, resolution):
|
||||
xp1=-aperture/2
|
||||
xp1=-aperture/2 #Set integration limits
|
||||
xp2=aperture/2
|
||||
|
||||
def yp1func(xp):
|
||||
return -np.sqrt((aperture/2)**2-(xp**2))
|
||||
return -np.sqrt((aperture/2)**2-(xp**2)) #Define y limits in terms of x
|
||||
|
||||
def yp2func(xp):
|
||||
return np.sqrt((aperture/2)**2-(xp**2))
|
||||
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to integrate for
|
||||
ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
|
||||
intensities = []
|
||||
|
|
@ -99,8 +94,8 @@ def plot2Dcircular(aperture, z, k, screen_range, resolution):
|
|||
for y in tqdm(ys):
|
||||
xIntensities = []
|
||||
for x in xs:
|
||||
realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1func, yp2func, args=(y, x, k, z))
|
||||
imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1func, yp2func, args=(y, x, k, z))
|
||||
realpart, realerror = integrate.dblquad(Fresnel2dreal, xp1, xp2, yp1func, yp2func, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10) #Calculate functions using circular limits
|
||||
imagpart, imagerror = integrate.dblquad(Fresnel2dimag, xp1, xp2, yp1func, yp2func, args=(y, x, k, z), epsabs=1e-10, epsrel=1e-10)
|
||||
|
||||
I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)
|
||||
xIntensities.append(I)
|
||||
|
|
@ -115,32 +110,32 @@ def plot2Dcircular(aperture, z, k, screen_range, resolution):
|
|||
plt.colorbar()
|
||||
plt.show()
|
||||
|
||||
def monte(aperture, z, k, screen_range, resolution, samples):
|
||||
N = samples
|
||||
def monte(aperture, z, k, screen_range, resolution, samples): #Function for part 4
|
||||
N = samples #Number of samples
|
||||
|
||||
def doubleInteg(x, y, xp, yp, z, k, aperture):
|
||||
def doubleInteg(x, y, xp, yp, z, k, aperture): #Function to return calculated values for each sample
|
||||
values = []
|
||||
for i in range(len(xp)):
|
||||
if (xp[i]**2+yp[i]**2) > (aperture/2)**2:
|
||||
if (xp[i]**2+yp[i]**2) > (aperture/2)**2: #Check if point is in the aperture
|
||||
values.append(0)
|
||||
else:
|
||||
value = np.exp(((1j*k)/(2*z))*((x-xp[i])**2+(y-yp[i])**2))
|
||||
value = np.exp(((1j*k)/(2*z))*((x-xp[i])**2+(y-yp[i])**2)) #Calculate value if in aperture
|
||||
values.append(value.imag)
|
||||
return np.array(values)
|
||||
|
||||
def monteCarlo(x, y, z, k, aperture):
|
||||
xp = np.random.uniform(low=(-aperture/2), high=aperture/2, size=N)
|
||||
def monteCarlo(x, y, z, k, aperture): #Perform onte carlo method
|
||||
xp = np.random.uniform(low=(-aperture/2), high=aperture/2, size=N) #Generate random samples
|
||||
yp = np.random.uniform(low=(-aperture/2), high=aperture/2, size=N)
|
||||
values = doubleInteg(x, y, xp, yp, z, k , aperture)
|
||||
values = doubleInteg(x, y, xp, yp, z, k , aperture) #Calculate value for each pair of samples
|
||||
mean = values.sum()/N
|
||||
meansq = (values*values).sum()/N
|
||||
integral = aperture*mean
|
||||
integral = aperture*mean #Calculate the integral end error
|
||||
error = aperture*np.sqrt((meansq-mean*mean)/N)
|
||||
return integral, error
|
||||
|
||||
def genData(aperture, z, k, resolution, screen_range):
|
||||
def genData(aperture, z, k, resolution, screen_range): ~Function to generate the data
|
||||
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution) #Generate values to integrate for
|
||||
ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
||||
|
||||
|
||||
|
|
@ -152,32 +147,28 @@ def monte(aperture, z, k, screen_range, resolution, samples):
|
|||
for y in tqdm(ys):
|
||||
xIntensities = []
|
||||
for x in tqdm(xs):
|
||||
integral, error = monteCarlo(x, y, z, k, aperture)
|
||||
integral, error = monteCarlo(x, y, z, k, aperture) #Find integral vaule using monte carlo method
|
||||
I = c*e0*constant*integral
|
||||
if I < 0.005:
|
||||
I = 0
|
||||
# completion = completion + 100/resolution**2
|
||||
# print(completion)
|
||||
xIntensities.append(I)
|
||||
intensities.append(xIntensities)
|
||||
intensities = np.array(intensities)
|
||||
return intensities
|
||||
|
||||
intensity = genData(aperture, z, k, resolution, screen_range)
|
||||
intensity = genData(aperture, z, k, resolution, screen_range) #Generate data
|
||||
extents = (-screen_range/2,screen_range/2,-screen_range/2,screen_range/2)
|
||||
|
||||
plt.imshow(intensity,vmin=1.0*intensity.min(),vmax=1.0*intensity.max(),extent=extents,origin="lower",cmap="nipy_spectral_r")
|
||||
plt.colorbar()
|
||||
plt.show()
|
||||
|
||||
MyInput = '0'
|
||||
MyInput = '0' #Selection menu
|
||||
while MyInput != 'q':
|
||||
MyInput = input('Enter a choice, "1", "2", "3", "4" or "q" to quit: ')
|
||||
print('You entered the choice: ',MyInput)
|
||||
if MyInput == '1':
|
||||
print('You have chosen part (1): 1D rectangular diffraction')
|
||||
aperture = input("Please input the desired aperture (m): ")
|
||||
z = input("Please enter the desired distnce from the screen (m): ")
|
||||
z = input("Please enter the desired distance from the screen (m): ")
|
||||
wl = input("Please enter the desired wavelength of light (m): ")
|
||||
k = (2*np.pi)/float(wl)
|
||||
screen_range = input("Please enter the diameter of the screen (m): ")
|
||||
|
|
@ -186,7 +177,7 @@ while MyInput != 'q':
|
|||
elif MyInput == '2':
|
||||
print('You have chosen part (2): 2D rectangular diffraction')
|
||||
aperture = input("Please input the desired aperture (m): ")
|
||||
z = input("Please enter the desired distnce from the screen (m): ")
|
||||
z = input("Please enter the desired distance from the screen (m): ")
|
||||
wl = input("Please enter the desired wavelength of light (m): ")
|
||||
k = (2*np.pi)/float(wl)
|
||||
screen_range = input("Please enter the diameter of the screen (m): ")
|
||||
|
|
@ -195,7 +186,7 @@ while MyInput != 'q':
|
|||
elif MyInput == '3':
|
||||
print('You have chosen part (3): 2D circular diffraction')
|
||||
aperture = input("Please input the desired aperture (m): ")
|
||||
z = input("Please enter the desired distnce from the screen (m): ")
|
||||
z = input("Please enter the desired distance from the screen (m): ")
|
||||
wl = input("Please enter the desired wavelength of light (m): ")
|
||||
k = (2*np.pi)/float(wl)
|
||||
screen_range = input("Please enter the diameter of the screen (m): ")
|
||||
|
|
@ -204,7 +195,7 @@ while MyInput != 'q':
|
|||
elif MyInput == '4':
|
||||
print('You have chosen part (3): 2D circular diffraction usinf Monte Carlo')
|
||||
aperture = input("Please input the desired aperture (m): ")
|
||||
z = input("Please enter the desired distnce from the screen (m): ")
|
||||
z = input("Please enter the desired distance from the screen (m): ")
|
||||
wl = input("Please enter the desired wavelength of light (m): ")
|
||||
k = (2*np.pi)/float(wl)
|
||||
screen_range = input("Please enter the diameter of the screen (m): ")
|
||||
|
|
|
|||
|
|
@ -85,6 +85,8 @@ To allow for this to be calculated computtionally, we will simplify it, integrat
|
|||
|
||||
\section{Explanation of Code}
|
||||
|
||||
|
||||
|
||||
\section{Results and Discussion}
|
||||
|
||||
\section{Conclusion}
|
||||
|
|
|
|||
Loading…
Add table
Add a link
Reference in a new issue