Page MenuHomec4science
Files F105184529

InterferenceCode.py
No OneTemporary
Actions

File Metadata

Created
Sat, Mar 15, 06:43

InterferenceCode.py
View Options

from __future__ import print_function
import matplotlib.pyplot as plt
import numpy as np
from matplotlib.widgets import Slider
from scipy.integrate import quad
import math
from scipy.integrate import simps
from ipywidgets import interact, interactive, fixed, interact_manual
import ipywidgets as widgets
import numpy as np
import matplotlib.pyplot as plt
from scipy.fftpack import *
def get_r(x, y, x0=0, y0=0):
return np.sqrt((x-x0)**2 + (y-y0)**2)
def point_source(x, y, x0=0, y0=0,lambda_=1.):
k = 2*np.pi / (lambda_)
return np.cos(k * get_r(x,y, x0, y0))
xmin = -1e-1
xmax = 1e-1
ymin = -1e-1
ymax = 1e-1
grid_spacing = 1 #mm
x = np.arange(-1e2*grid_spacing, 1e2*grid_spacing, grid_spacing) #m
y = x #m
X, Y = np.meshgrid(x, y)
def n_slits(n, xmin=xmin, xmax=xmax, ymin=ymin, ymax=ymax, y0=0, lambda_=1.):
"""Compute the interference pattern for n equidistant sources from xmin
to xmax in units of the grid_spacing. The wavelength is also given in
grid_spacing units."""
x0 = np.linspace(xmin*grid_spacing, xmax*grid_spacing, n)
y0*=grid_spacing
final = np.zeros_like(X)
for i in x0:
final += point_source(X, Y, x0=i, y0=y0, lambda_=lambda_*grid_spacing)
return np.abs(final)**2
def plot_point_sources(n_src, lambda_exp, y_slice, x_slice, d_sources):
l = 10**(lambda_exp) * grid_spacing #mm
lambda_ = 10 ** lambda_exp ## grid_spacings
# input d_sources is in lambdas
d_sources = float(d_sources * lambda_) # grid_spacings
d_sources_lambda = d_sources / lambda_ #lambdas
xrange = 0.5 * (n_src - 1) * d_sources #grid_spacings
sources = np.arange(-xrange, xrange+1, d_sources)*grid_spacing
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
intensity = n_slits(n_src, lambda_=lambda_, xmin=-xrange, xmax=xrange) # Accepts quantities in grid spacings
ax[0].pcolorfast(x , y, intensity)
ax[0].set_xlabel('$x/\lambda$')
ax[0].set_ylabel('$y/\lambda$')
ax[0].axvline(x[x_slice],
c='r',
lw=3,
label='x = %.2f mm' % (x[x_slice]))
ax[0].axhline(y[y_slice],
c='magenta',
lw=3,
label='y = %.2f mm' % (y[y_slice]))
ax[0].scatter(sources*grid_spacing,
np.zeros_like(sources),
color='white',
s=50,
marker='*',
label='d = %.2f$\lambda$' %
(d_sources_lambda))
ax[1].plot(x , intensity[y_slice, :], c='magenta')
ax[1].set_xlabel('$x/\lambda$')
ax[2].plot(y , intensity[:, x_slice], c='r')
ax[2].set_xlabel('$y/\lambda$')
ax[1].set_ylabel('I [a.u.]')
ax[2].set_ylabel('I [a.u.]')
ax[1].set_ylim(0, n_src**2)
ax[2].set_ylim(0, n_src**2)
ax[0].set_xlim(x[0], x[-1])
ax[0].legend(loc=(-0.7, 0.8))
for a in ax.ravel():
labels = a.get_xticks()
a.set_xticklabels([float('%.2f'%(label/l)) for label in labels])
labels = ax[0].get_yticks()
ax[0].set_yticklabels([float('%.2f'%(label/l)) for label in labels])
ax[0].set_title('Zoom = %.2fX' % (np.log10(lambda_)))
def slit_pattern(k, x, y, a=0.1, b=0.1):
# return (np.sin(k*x*a)*np.sin(k*y*b)/(k*x*a*k*y*b))**2
return (np.sinc(k*x*a)*np.sinc(k*y*b))**2
x = np.linspace(-50, 50, 100)
y = x
X,Y = np.meshgrid(x,y)
def plot_sinc_2d(lambda_exp, a, b):
"""Plots the 2D sinc function corresponding to the interference pattern
from a slit of dimensions 2a*2b when hit by light of wavelength 10**lambda_exp."""
lambda_ = 10**lambda_exp
k = 2 * np.pi / lambda_
fig, ax = plt.subplots(1, 2)
intensity = slit_pattern(k, X, Y, a=a, b=b)
intensity/=np.max(intensity)
ax[0].pcolorfast(x, y, intensity)
ax[1].plot(x, intensity[int(0.5 * intensity.shape[1]), :])
fig.tight_layout()
#Pupil function
def draw_rectangle_pupil(aperture_fraction_x,
aperture_fraction_y,
canvas,
x_pos=0.5,
y_pos=0.5):
HEIGHT, WIDTH = canvas.shape
aperture_indices_x = WIDTH * np.array([
x_pos - 0.5 * aperture_fraction_x, x_pos + 0.5 * aperture_fraction_x
])
aperture_indices_y = HEIGHT * np.array([
y_pos - 0.5 * aperture_fraction_y, y_pos + 0.5 * aperture_fraction_y
])
masked = np.copy(canvas)
masked[int(aperture_indices_y[0]):int(aperture_indices_y[1]),
int(aperture_indices_x[0]):int(aperture_indices_x[1])] = 1
return masked
def draw_circular_pupil(radius, canvas, x_pos = 0.5, y_pos = 0.5):
HEIGHT, WIDTH = canvas.shape
x = np.arange(HEIGHT)
y = np.arange(WIDTH)
X, Y = np.meshgrid(x, y)
x_pos_index = int(WIDTH * x_pos)
y_pos_index = int(WIDTH * y_pos)
distances = get_r(X, Y, x_pos_index, y_pos_index)
mask = np.copy(canvas)
mask[distances < radius*WIDTH] = 1
return mask
def interference_from_pupil(pupil):
HEIGHT, WIDTH = pupil.shape
# Intensity
pupil_fft = np.abs(fft2(pupil))**2
# Shift frequencies
pupil_fft_shift = fftshift(pupil_fft)
# Plot
fig, ax = plt.subplots(1, 3, figsize=(15, 5))
ax[0].imshow(pupil, cmap="Greys_r")
ax[1].imshow(np.log((pupil_fft_shift + 1))) # Log stretch
ax[1].set_title("Log stretched interference pattern.")
ax[2].plot(pupil_fft_shift[int(0.5 * HEIGHT)]/np.max(pupil_fft_shift))
ax[2].set_title('$I/I_0$')
[a.set_xlabel('x [px]') for a in ax]
ax[0].set_ylabel('y [px]')
ax[1].set_ylabel('y [px]')
plt.tight_layout()
def plot_interact_circular(r):
"""Plot the interference pattern for a circular aperture of radius r in
units of WIDTH"""
HEIGHT = 1000
WIDTH = HEIGHT
canvas = np.zeros([HEIGHT, WIDTH])
interference_from_pupil(draw_circular_pupil(r, canvas))
def plot_interact_rectangular(h_x, h_y):
"""Plot the interference pattern for a rectangular aperture of width h_x
and height h_y in units of WIDTH"""
HEIGHT = 1000
WIDTH = HEIGHT
canvas = np.zeros([HEIGHT, WIDTH])
interference_from_pupil(draw_rectangle_pupil(h_x, h_y, canvas))
def plot_interact_slits(n_slits, d_slits, h_x, h_y):
"""Plot the interference pattern for n_slits slits of width h_x
and height h_y, separated by a distance d_slits, all in units of WIDTH"""
HEIGHT = 1001
WIDTH = HEIGHT
canvas = np.zeros([HEIGHT, WIDTH])
xrange = 0.5 * (n_slits - 1) * d_slits
x_sources = np.linspace(0.5 - xrange, 0.5 + xrange, n_slits)
for x in x_sources:
canvas = draw_rectangle_pupil(h_x, h_y, canvas, x_pos=x, y_pos=0.5)
interference_from_pupil(canvas)
def plot_interact_holes(n_holes, d_holes, r):
"""Plot the interference pattern for n_holes circular holes of radius r,
separated by a distance d_holes, all in units of WIDTH"""
HEIGHT = 1000
WIDTH = HEIGHT
canvas = np.zeros([HEIGHT, WIDTH])
total_interval = n_holes*d_holes
xrange = 0.5 * (n_holes - 1) * d_holes
x_sources = np.linspace(0.5 - xrange, 0.5 + xrange, n_holes)
for x in x_sources:
canvas = draw_circular_pupil(r, canvas, x_pos=x, y_pos=0.5)
interference_from_pupil(canvas)
def diffraction_intensity (slit_width, wavelength, screen_distance, distance_between_slits, number_of_slits, X):
te1 = np.sin(np.pi*X*slit_width*10**(-6)/(wavelength*(10**-9)*screen_distance*10**(-2)))/(np.pi*X*slit_width*10**(-6)/(wavelength*(10**-9)*screen_distance*10**(-2)))
te2 = (np.sin(number_of_slits*np.pi*distance_between_slits*10**(-3)*X/(wavelength*(10**-9)*screen_distance*10**(-2))))/(number_of_slits*np.sin((np.pi*distance_between_slits*10**(-3)*X)/(wavelength*(10**-9)*screen_distance*10**(-2))))
intensity = (te1**2)*(te2**2)
plt.plot(X, intensity)
plt.xlabel("y (m)")
plt.ylabel("Relative intensity")
plt.show()

Event Timeline

c4science · Help