216 lines
8.8 KiB
Python
216 lines
8.8 KiB
Python
import matplotlib.pyplot as plt
|
|
import numpy as np
|
|
from scipy import integrate
|
|
from tqdm import tqdm
|
|
|
|
def Fresnel2dreal(yp, xp, y, x, k, z):
|
|
kernel = np.cos((k/(2*z))*((x-xp)**2+(y-yp)**2))
|
|
return kernel
|
|
|
|
def Fresnel2dimag(yp, xp, y, x, k, z):
|
|
kernel = np.sin((k/(2*z))*((x-xp)**2+(y-yp)**2))
|
|
return kernel
|
|
|
|
c = 3e8
|
|
e0 = 8.85e-12
|
|
|
|
def plot1D(aperture, z, k, screen_range, resolution):
|
|
def genData(aperture, z, k, screen_range, resolution):
|
|
y = 0
|
|
|
|
xp1=yp1=-aperture/2
|
|
xp2=yp2=aperture/2
|
|
|
|
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
intensities = []
|
|
|
|
constant = k/2*np.pi*z
|
|
completion = 0
|
|
|
|
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)
|
|
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)
|
|
return xs, intensities
|
|
|
|
ax = plt.axes()
|
|
xs, intensities = genData(aperture, z, k, screen_range, resolution)
|
|
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
|
|
xp2=yp2=aperture/2
|
|
|
|
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
|
|
intensities = []
|
|
|
|
constant = k/2*np.pi*z
|
|
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)
|
|
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)
|
|
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)
|
|
|
|
plt.imshow(intensity,vmin=0.0,vmax=1.0*intensity.max(),extent=extents,origin="lower",cmap="nipy_spectral_r")
|
|
plt.colorbar()
|
|
plt.show()
|
|
|
|
def plot2Dcircular(aperture, z, k, screen_range, resolution):
|
|
def genData(aperture, z, k, screen_range, resolution):
|
|
xp1=-aperture/2
|
|
xp2=aperture/2
|
|
|
|
def yp1func(xp):
|
|
return -np.sqrt((aperture/2)**2-(xp**2))
|
|
|
|
def yp2func(xp):
|
|
return np.sqrt((aperture/2)**2-(xp**2))
|
|
|
|
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
|
|
intensities = []
|
|
|
|
constant = k/2*np.pi*z
|
|
|
|
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))
|
|
|
|
I = c*e0*((realpart*constant)**2+(imagpart*constant)**2)
|
|
xIntensities.append(I)
|
|
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)
|
|
|
|
plt.imshow(intensity,vmin=0.0,vmax=1.0*intensity.max(),extent=extents,origin="lower",cmap="nipy_spectral_r")
|
|
plt.colorbar()
|
|
plt.show()
|
|
|
|
def monte(aperture, z, k, screen_range, resolution, samples):
|
|
N = samples
|
|
|
|
def doubleInteg(x, y, xp, yp, z, k, aperture):
|
|
values = []
|
|
for i in range(len(xp)):
|
|
if (xp[i]**2+yp[i]**2) > (aperture/2)**2:
|
|
values.append(0)
|
|
else:
|
|
value = np.exp(((1j*k)/(2*z))*((x-xp[i])**2+(y-yp[i])**2))
|
|
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)
|
|
yp = np.random.uniform(low=(-aperture/2), high=aperture/2, size=N)
|
|
values = doubleInteg(x, y, xp, yp, z, k , aperture)
|
|
mean = values.sum()/N
|
|
meansq = (values*values).sum()/N
|
|
integral = aperture*mean
|
|
error = aperture*np.sqrt((meansq-mean*mean)/N)
|
|
return integral, error
|
|
|
|
def genData(aperture, z, k, resolution, screen_range):
|
|
|
|
xs = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
ys = np.linspace(-screen_range/2, screen_range/2, num=resolution)
|
|
|
|
|
|
intensities = []
|
|
completion = 0
|
|
|
|
constant = k/(2*np.pi*z)
|
|
|
|
for y in tqdm(ys):
|
|
xIntensities = []
|
|
for x in tqdm(xs):
|
|
integral, error = monteCarlo(x, y, z, k, aperture)
|
|
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)
|
|
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'
|
|
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): ")
|
|
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): ")
|
|
resolution = input("Please enter the resolution of the plot (pixels): ")
|
|
plot1D(float(aperture), float(z), float(k), float(screen_range), int(resolution))
|
|
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): ")
|
|
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): ")
|
|
resolution = input("Please enter the resolution of the plot (pixels): ")
|
|
plot2Drectangular(float(aperture), float(z), float(k), float(screen_range), int(resolution))
|
|
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): ")
|
|
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): ")
|
|
resolution = input("Please enter the resolution of the plot (pixels): ")
|
|
plot2Dcircular(float(aperture), float(z), float(k), float(screen_range), int(resolution))
|
|
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): ")
|
|
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): ")
|
|
resolution = input("Please enter the resolution of the plot (pixels): ")
|
|
samples = input("Please enter the desired nu,ber of samples for the calculation: ")
|
|
monte(float(aperture), float(z), float(k), float(screen_range), int(resolution), int(samples))
|
|
elif MyInput != 'q':
|
|
print('This is not a valid choice')
|
|
print('You have chosen to finish - goodbye.')
|