diff --git a/Exercise 2/exercise2.py b/Exercise 2/exercise2.py index 1c4817a..b6c5fc8 100644 --- a/Exercise 2/exercise2.py +++ b/Exercise 2/exercise2.py @@ -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): ") diff --git a/Exercise 2/main.tex b/Exercise 2/main.tex index b8c5b1f..f1f1bd7 100644 --- a/Exercise 2/main.tex +++ b/Exercise 2/main.tex @@ -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}