# Propagation of a Gaussian beam using FFT

I am trying to simulate the propagation of a gaussian beam through a lens using an FFT approach.

I tried to implement the approach described by Couairon in this paper at page 43: https://link.springer.com/article/10.1140/epjst/e2011-01503-3

Here is a description of the algorithm: .

However, when I propagate the beam from a lens to the focal spot, the intensity I get is much less than if the theoretical intensity that I can compute using the theoretical formulas(https://en.wikipedia.org/wiki/Gaussian_beam):

$$I(r,z)=I_0 (\frac{w_0}{w(z)})^2 e^{(-2r^2/w^2(z)}$$

where $$I_0$$ is the intensity at the center of the beam at its waist.

Here is the code I wrote:

Wavelength=532e-9
n0=1.0 #Freespace index of refraction
f=100e-3 # Lens focal length
r_max=10e-3 # Limits of the grid
w_lens=1e-3 #Waist size at the lens
N_r=2**12 # Number of points
r_grid = np.linspace(-r_max, r_max, N_r)  # Spatial grid
k0 = (n0) * 2 * np.pi / Wavelength  # Freespace wavevector
dr = abs(r_grid[1] - r_grid[0])

E_initial= 1* np.exp(-(r_grid ** 2) / (w_lens ** 2) - 1j * k0 * r_grid ** 2 / (2 * f))
dz=f  # Distance of propagation (propagating to the focal spot)
delta = dz / (4 * k0 * dr ** 2)
# Steps described in last item of step 2
kjdr1 = [2 * np.pi * j / N_r for j in np.arange(0, N_r / 2)]
kjdr2 = [2 * np.pi * (-1 + j / N_r) for j in np.arange(N_r / 2, N_r)]
kjdr = np.concatenate([kjdr1, kjdr2])

A = np.exp(-2 * 1j * delta * (kjdr) ** 2)

E_F = np.fft.fft(np.fft.fftshift(E_initial),norm="ortho")
E_F = (E_F * A) # Convolution
E_Final = np.fft.ifftshift(np.fft.ifft(E_F,norm="ortho"))


If you use these parameters, the ratio $$\frac{I_{focus}}{I_{lens}}$$ is 56.64 . However, if you use the theoretical relation, you get 3488.20.

Do you see anything wrong in my implementation of the algorithm? Anything you see that could explain why the ratio $$\frac{I_{focus}}{I_{lens}}$$ seems to be way off?

PS:My code uses field E but I convert E to I using $$I=2*n_0*\epsilon_0*c*|E|^2$$

• Maybe because in theory it is a 2-D problem but here you implemented it in 1-D? So I guess there might be some relationship like a square somewhere. BTW, I am doing similar things here, and thank you so much for sharing the link of the paper, it really helps a lot! – Wei Li Jul 3 '20 at 9:53