[Python] Issues with boundary layers in step propagation

I am trying to create a function which receives the transversal profile of a beam and outputs the transversal profile propagated by some distance. To exemplify the task, here is the profile of a Gaussian beam at z=z_0:



After propagating this profile a distance d=z_0, it is expected that the propagated profile matches the profile at z'=z+d=2 z_0:



My function works while the profile fits in the grid, but produces errors once the propagated beam is larger than the grid. That is an expected error. To avoid it, this article recommends that I implement a propagation in steps, with absorbing boundary layers that prevent the intermediate profiles from having significant intensity next to the borders while also not introducing frequencies in the Fourier transform of the profiles that may disturb propagation.

Following the article, I am trying to implement step propagation (with fixed step size, for now) that propagates a gaussian beam correctly. Here is my code (apologies in advance for the names and comments in Portuguese):

import numpy as np
import scipy
from matplotlib import pyplot as plt
from numpy.fft import fft, ifft, fft2, ifft2, fftfreq, fftshift, ifftshift
from scipy.special import genlaguerre

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 feixe:
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

#funções que dependem de z
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): #fase de Gouy
return arctg(z, self.z_0)

class FeixeGaussiano(feixe):
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) #garantindo que z é um array

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

def AmplitudeTransversal(fx, z: float = 0.0):
W_0=fx.W_0
gr=1000
X = np.linspace(-4*W_0,4*W_0,gr)
Y = np.linspace(-4*W_0,4*W_0,gr)
x,y=np.meshgrid(X,Y)
rho=sqrt(x**2+y**2)
phi=arctg(y,x)

u=fx.Amplitude(z,rho,phi) #Amplitude

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()

def plota(u,X,Y): #plota uma matriz qualquer. Tipicamente, queremos plotar perfis transversais propagados
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()


compOnda=800*(10**-9) #nm
z_0=pi*(10**(-3))/2 #m -> valor tal que W_0=20.10⁻6

def passo(U,dz,kz):
#amplitude complexa radial
Fourier_U = fft2(U) #espectro angular
if dz>0.0:
Propag_FU = Fourier_U * exp(1j * kz * dz).conj() #propaga pra frente
elif dz<0.0:
Propag_FU = Fourier_U * exp(1j * kz * (-dz)) #propaga pra trás
U_propag = ifft2(Propag_FU)
return U_propag


def propaga(U,z,X,Y,k,filtro):
#adicionamos células de fronteira ao perfil
#nossa grade é 1250x1250, com 125 células de fronteira em cada borda.
drho=0.2 #pois a fronteira equivale a 20% da grade. implementar drho dependente do perfil depois

#Coordenadas de frequência
grx,gry = U.shape
dx=(X[-1]-X[0])/(grx-1)
dy=(Y[-1]-Y[0])/(gry-1)
fx = fftfreq(grx, d=dx)
fy = fftfreq(gry, d=dy)
fx, fy = np.meshgrid(fx, fy)

#constantes de propagação:
kx = 2*pi*fx
ky = 2*pi*fy
kz = sqrt(k**2 - kx**2 - ky**2 + 0j)
krho=sqrt(kx**2+ky**2)

U_passo=U

#calculamos a distância máxima de propagação em z que não ultrapassa as bordas:
while z!=0.0:
U_=U_passo*filtro
#dz=max((drho*kz/krho).flatten()) #menor dz é o máximo que dá pra propagar sem ter erros #maximo, minimo...
dz=0.00005
dz=min(dz,z) #se z for menor que a distância máxima de propagação sem erro, propagamos z. (por enquanto só dá pra propagar pra frente)
U_passo=passo(U_,dz,kz)
z-=dz

return U_passo

#bordas gaussianas
fx=FeixeGaussiano(z_0,compOnda)
W_0=fx.W_0
k=fx.k
I_0=fx.I_0
b=500
gr=1000+2*b
X = np.linspace(-8*W_0,8*W_0,gr)
Y = np.linspace(-8*W_0,8*W_0,gr)
x,y=np.meshgrid(X,Y)
rho=sqrt(x**2+y**2)

phi=arctg(y,x)

def gaussiana(x,b):
sigma = 1
x-=4*W_0
x/=W_0
gaussiana = exp(-(x**2) / (2 * sigma**2))
return gaussiana

filtro=np.where(rho>4*W_0, gaussiana(rho,b),1)

x = np.linspace(-8*W_0,8*W_0,2000)
y = np.linspace(-8*W_0,8*W_0,2000)
X,Y = np.meshgrid(x,y)
rho=sqrt(X**2+Y**2)
phi=arctg(Y,X)
U=fx.Amplitude(z_0,rho,phi)

U_propag=propaga(U,5*z_0,x,y,k,filtro)
x_ = np.linspace(-4*W_0,4*W_0,1000)
y_ = np.linspace(-4*W_0,4*W_0,1000)
plota(U_propag[-1500:1500,-1500:1500],x_,y_)
AmplitudeTransversal(fx,6*z_0)


The code will output the transversal profile from z_0 propagated 5*z_0 and the analytical transversal profile in 6*z_0, respectively. It is expected that both are the same, but they aren't. Other large propagation distances also give incorrect outputs, but small ones (such as z_0) don't.

Since wrong outputs are only observed once the beam reaches the borders, I have analysed my filter in hopes of finding what is wrong, to no success. Below are some graphs of it.





As far as I could check, the filter should work according to the article I'm following. Changing the step size does not work, and neither does increasing the size of the boundary layer. What am I doing wrong?

Continue reading...