I am trying to re-purpose a Matlab code re-written in python for the diffraction on a circular aperture shifted with respect to the center of the image. I have looked at earlier posts and did not find a similar one. I am also looking at Computational Fourier Optics as a reference. I expect to get a maximum dislocated from the center in the same direction as the displacement of the initial mask. However, I get the diffraction image that is in the center no matter what the location of the circular aperture is. Is it me needing to go and read undergrad Phys or something is wrong with the code? Thanks in advance! The code is below.
from numba import jit
import scipy
import numpy as np
import time
import matplotlib
import matplotlib.pyplot as plt
import matplotlib.colors as colors
import math
from math import sqrt
padding = 0
complex_type = np.complex128
float_type = np.float64
lambda_0 = 0.5e-6
f = 10e-2
(4*lambda_0*f**3/np.pi)**(1/4)
def propFF(u1,L1,lambda_0,z):
# Fraunhofer Propagation kernel
u1 = u1.astype(complex_type)
(M,N)=np.shape(u1)
L1 = float_type(L1)
lambda_0 = float_type(lambda_0)
z = float_type(z)
M = float_type(M)
N = float_type(N)
dx1 = float_type(L1/M)
k = float_type(2*np.pi/lambda_0)
L2 = float_type(lambda_0*z/dx1)
dx2 = float_type(lambda_0*z/L1)
x2 = np.linspace(-L1/float_type(2),L1/float_type(2),int(M+2 * padding))
(X2,Y2) = np.meshgrid(x2,x2)
X2 = X2.astype(float_type)
Y2 = Y2.astype(float_type)
c = float_type(1)/(complex_type(1j)*lambda_0*z)\
*np.exp(complex_type(1j)*k/(float_type(2)*z)*(X2**2+Y2**2))
u2=c*np.fft.ifftshift(np.fft.fft2(np.fft.fftshift(u1),s=[int(M)+2*padding, int(M)+2*padding]))*(dx1**2)
return u2
def rect(x, y, x_0, x_1, y_0, y_1):
X = np.zeros(np.shape(x)).astype(float)
X[np.abs(x-x_0) < 0.5] = 1.0
X[np.abs(x-x_1) == 0.5] = 0.5
Y = np.zeros(np.shape(y)).astype(float)
Y[np.abs(y-y_0) < 0.5] = 1.0
Y[np.abs(y-y_1) == 0.5] = 0.5
return(X * Y)
def circle(x, y, x_0, x_1, y_0, y_1):
R = np.zeros(np.shape(x)).astype(float)
R[(np.abs(x-x_0)**2+np.abs(y-y_0)**2)**.5 < 0.5] = 1.0
R[(np.abs(x-x_1)**2+np.abs(y-y_1)**2)**.5 == 0.5] = 0.5
return(R)
L1 = float(0.5)
M=float(160)
lambda_0=0.5e-6
k=2*np.pi/lambda_0
w=0.051
z=2000
dx1=L1/M
x1=np.linspace(-L1/2,L1/2,np.int(L1/dx1))
y1=x1
(X1,Y1)=np.meshgrid(x1,y1)
u1 = circle(X1/(2*w), Y1/(2*w), 1., 0, 0, 0)
plt.figure(figsize=(4.5,4))
plt.pcolormesh(x1,y1,np.abs(u1**2))
u2_tf = propFF(u1,L1,lambda_0,z)
I2_tf = np.real(u2_tf*np.conj(u2_tf))
x1 = np.linspace(-L1/2, L1/2,np.int(M + padding * 2))
y1 = x1
plt.figure(figsize=(8,7))
plt.pcolormesh(x1,y1,I2_tf**(1.0/3.0))
plt.colorbar()
plt.title('Fresnel Transfer Function')
plt.show()