I am trying to write a code which propagates a beam a certain distance through Fourier propagation. Here is the relevant excerpt of the code:
import numpy as np
import scipy
from matplotlib import pyplot as plt
from numpy.fft import fft, ifft, fft2, ifft2, fftfreq, fftshift, ifftshift
exp=np.exp
pi=np.pi
sqrt=np.sqrt
tg=np.tan
arctg=np.arctan2
angle=np.angle
array=np.array
plt.rcParams["figure.autolayout"] = True
class beam:
def __init__(self, z_0: float, compOnda: float):
self.z_0 = z_0 #distância de Rayleigh
self.compOnda = compOnda #comprimento de onda
self.k = 2*pi / self.compOnda #número de onda
self.W_0 = sqrt(self.compOnda * self.z_0 / pi)
self.teta_0 = self.W_0 / self.z_0
self.A_0 = sqrt(2 / (pi * self.W_0)) #amplitude em 0; definimos ela de forma a normalizar o pulso (gaussiano - preciso verificar se normaliza o laguerre gaussiano também)
self.I_0 = self.A_0 * self.A_0.conjugate() #intensidade em 0
def W(self,z: np.array):
return self.W_0 * sqrt(1 + (z / self.z_0)**2)
def R(self,z: np.array):
z_=np.where(z==0.0,1e-20,z) #evitando divisão por zero
R=z_ * (1 + (self.z_0 / z_)**2)
R=np.where(z==0.0,np.inf,R)
return R
def zeta(self,z: np.array):
return arctg(z, self.z_0)
class GaussianBeam(beam):
def __init__(self, z_0: float, compOnda: float):
super().__init__(z_0,compOnda)
def Amplitude(self,z: array = 0.0, rho: array = 0.0, phi: array = 0.0) -> array:
z=np.asarray(z,dtype=float)
z_0=self.z_0
A_0=self.A_0
k=self.k
W_0=self.W_0
W=self.W(z)
R=self.R(z)
zeta=self.zeta(z)
U=A_0*(W_0/W)*exp(-rho**2/W**2 + (-k*z-k*(rho**2)/(2*R)+zeta)*1.0j)
return U
wavelength=800*(10**-9) #nm
z_0=pi*(10**(-3))/2 #m ->approx. 0.0015. W_0=20.10⁻6
dz=0.006 #propagation distance (m)
def Propag(U,dz,gr,d):
#frequency coordinates
fx = fftfreq(gr, d=d)
fy = fftfreq(gr, d=d)
fx, fy = np.meshgrid(fx, fy)
kx = 2*pi*fx
ky = 2*pi*fy
kz = sqrt(k**2 - kx**2 - ky**2 + 0j) #propagation constant in z
#radial complex amplitude
Fourier_U = fft2(U) #angular spectrum
Propag_FU = Fourier_U * exp(1j * kz * dz)
U_propag = ifft2(Propag_FU)
return U_propag
fx=GaussianBeam(z_0,wavelength)
W_0=fx.W_0
gr=1000
k=fx.k
x = np.linspace(-4*W_0,4*W_0,gr) #grid size
y = np.linspace(-4*W_0,4*W_0,gr)
d=8*W_0/(gr-1)
X,Y = np.meshgrid(x,y)
rho=sqrt(X**2+Y**2)
phi=arctg(Y,X)
U=fx.Amplitude(z_0,rho,phi) #this is a bidimensional Gaussian
U_propag=Propag(U,dz,gr,d)
def plot(u):
modulo = abs(u)
fase = angle(u)
# Gráfico
plt.figure(figsize=(12, 5))
plt.subplot(1, 2, 1)
plt.pcolor(X, Y, modulo, shading='auto', cmap='viridis')
plt.title('Módulo')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
plt.subplot(1, 2, 2)
plt.pcolor(X, Y, fase, shading='auto', cmap='twilight')
plt.title('Fase')
plt.xlabel('x')
plt.ylabel('y')
plt.colorbar()
plt.tight_layout()
plt.show()
plot(U_propag)
The code works, but it propagates the beam backwards - which can be seen by the shrunkage of the profile when dz=z_0, meaning that the beam is propagating towards its waist. I want to fix this and make the beam propagate forward.
I read that the Fourier transform in optics has a positive sign, but swapping fft and ifft produced the same results. I then tried to give exp(1j * kz * dz) a negative sign, but that gave me an error:
/tmp/ipython-input-17-2768534464.py:14: RuntimeWarning: overflow encountered in exp
Propag_FU = Fourier_U * exp(-1j * kz * dz)
/usr/local/lib/python3.11/dist-packages/numpy/fft/_pocketfft.py:94: RuntimeWarning: invalid value encountered in fft
return ufunc(a, fct, axes=[(axis,), (), (axis,)], out=out)
Why is my beam being propagated backwards, and how can I propagate it forward?
exp(1j * kz * dz)to exp(-1j * kz * dz). Otherwise no error occurs except for the beam propagating backwards.