EDIT: Now I think it is a matter of convention, I had to use the right wiki page :)
Hey there,
First of all, thank you for putting this great code online!
I've tried to do some sanity test with your example and see the diffraction of a gaussian beam. Zr is known to be k0 w0^2 / 2. Which becomes k0(Dr/2.355)**2 / 2 = 0.0116 (the 2.355 factor is to convert fwhm to waist) in your example. Then I propagate the beam to this distance but unfortunately I do not get the sqrt(2) factor in the waist as the picture below shows:
I tried also different grids, wavelengths and resolutions with no success of getting the right behaviour at Zr. Do you have a hint what do I miss?
This is the code I used, I did it with minimal changes of your example file:
import matplotlib.pyplot as plt
import numpy as np
plt.figure()
plt.plot(r, Irz_norm[:,0])
plt.plot(r, Irz_norm[:,-1])
plt.title('Comparison of Intensities at initial position and at Zr')
plt.legend(['Initial', 'Zr'])
EDIT: Now I think it is a matter of convention, I had to use the right wiki page :)
Hey there, First of all, thank you for putting this great code online! I've tried to do some sanity test with your example and see the diffraction of a gaussian beam. Zr is known to be k0 w0^2 / 2. Which becomes k0(Dr/2.355)**2 / 2 = 0.0116 (the 2.355 factor is to convert fwhm to waist) in your example. Then I propagate the beam to this distance but unfortunately I do not get the sqrt(2) factor in the waist as the picture below shows:
I tried also different grids, wavelengths and resolutions with no success of getting the right behaviour at Zr. Do you have a hint what do I miss?
This is the code I used, I did it with minimal changes of your example file: import matplotlib.pyplot as plt import numpy as np
from pyhank import HankelTransform
def gauss1d(x, x0, fwhm): return np.exp(-2 np.log(2) ((x - x0) / fwhm) ** 2)
nr = 1024 # Number of sample points r_max = 5e-3 # Maximum radius (5mm) r = np.linspace(0, r_max, nr)
Nz = 2 # Number of z positions z_max = 0.0116 # Maximum propagation distance z = np.linspace(0, z_max, Nz) # Propagation axis
Dr = 100e-6 # Beam radius (100um) lambda = 488e-9 # wavelength 488nm k0 = 2 * np.pi / lambda # Vacuum k vector print('w0 is ' + str(k0*(Dr/2.355)**2/2))
H = HankelTransform(order=0, radial_grid=r)
Er = gauss1d(r, 0, Dr) # Initial field ErH = H.to_transform_r(Er) # Resampled field
EkrH = H.qdht(ErH)
Erz = np.zeros((nr, Nz), dtype=complex) kz = np.sqrt(k0 2 - H.kr 2) for i, z_loop in enumerate(z): phi_z = kz z_loop # Propagation phase EkrHz = EkrH np.exp(1j * phi_z) # Apply propagation ErHz = H.iqdht(EkrHz) # iQDHT Erz[:, i] = H.to_original_r(ErHz) # Interpolate output Irz = np.abs(Erz) ** 2
Irz_norm = Irz / Irz[0, :]
plt.figure() plt.plot(r, Irz_norm[:,0]) plt.plot(r, Irz_norm[:,-1]) plt.title('Comparison of Intensities at initial position and at Zr') plt.legend(['Initial', 'Zr'])