Fraunhofer Diffraction -- possibly shifting of fft2 issue

Viewed 212

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()
0 Answers
Related