Bessel Beam Simulation from an annulus of light

[DerKamin]

Hi,

I have some difficulties with the simulation of the Bessel Beam simulation. It would be glad if you could help me out or give a few tips.

Here is what I try to achieve: My goal is to investigate the propagation of a bessel beam resulting from the focusing of an annulus of light. Especially the width of the central lobe and the propagation distance dependent on the ring width and diameter is very interesting for me. Meanwhile I create a collimated annulus of light with some axicons lenses / and an SLM in the lab and will see if the results are comparable…

Here is my Code so far:

from LightPipes import* 
import matplotlib.pyplot as plt
import numpy as np
from tifffile import imsave

GridSize= 10*mm
GridDimension=1500
BesselGridDim=200
lambda_=800*nm
innerR=350
outerR=380
center_x=GridDimension//2
center_y=GridDimension//2
f=20*cm
dz=0.05*mm
nsteps=400

#Create an annulus of light
Ring=np.zeros([GridDimension,GridDimension])
for a in range(0,GridDimension):
    for b in range(0,GridDimension):
        distance=np.sqrt((a-center_x)**2+(b-center_y)**2)

        if innerR <= distance <= outerR:
            Ring[a,b]=1

#Create and propagate field 180mm     
Field=Begin(GridSize,lambda_,GridDimension)        
Field = SubIntensity(Field,Ring)
Field = Lens(f,0,0,Field)
Field=Forward(Field, 190*mm, 0.2*mm, BesselGridDim)

#Go through Bessel region and safe Intensity profile
FocusInten=np.zeros([nsteps,BesselGridDim,BesselGridDim])
for a in range(0,nsteps):
    Field=Steps(dz, Field)
    I=Intensity(0,Field)
    FocusInten[a,:,:]=I

#Normalizing and safe
Focusnormalized=np.array([nsteps,BesselGridDim,BesselGridDim],dtype=type(FocusInten))
Focusnormalized=(FocusInten - np.min(FocusInten))/(np.max(FocusInten)-np.min(FocusInten))
FocusInten=(Focusnormalized*65535).astype(np.uint16)
imsave('BesselRegion.tiff',FocusInten)

1. What are the best ways to propagate the light? The problem is that I have a large annulus of light (several millimeters) and very small focal region (several µm)
2. Is there any faster implementation?
3. At the moment, some errors occur since the Bessel beam seams to disappear in focus region. Any suggestion why are welcome.

Thank you very much in advance.

Greetings from Hannover😊

Edited by ldoyle: used code markup

[ldoyle]

Hi DerKamin,

I had a look at your code and I think I managed to get an acceptable solution. Some quick comments first:

• The generation of the annular beam can be done conveniently using LP functions CircScreen and CircAperture. In my example below I therefore omitted your loops and changed the radii definitions to meters instead of pixels.
• I simplified the final intensity normalization assuming the typical minimum should be 0.

As you noticed, using the Forward propagator is extremely slow (I actually did not test your code to more than BesselGridDim=50 because I was impatient 😉 ). Since you are using a lens, your example fits well with the spherical coordinates found in LightPipes [1]. This will allow you to change the grid size while using the faster propagators Fresnel or Forvard (as LensFresnel or LensForvard in this case).

A second problem I noticed is that independent of the method I tested, your BesselGridSize=0.2*mm is too small. There is significant intensity at the borders of the resulting grid, which create problems in further propagation using Steps (but also any other method, intensity should always drop off to almost 0 at the edges).

So I increased the BesselGridSize=1*mm and now it looks alright:

from LightPipes import * 
import matplotlib.pyplot as plt
import numpy as np
from tifffile import imsave
from tqdm import tqdm # fancy progress bar, recomment installing if you don't have it yet

GridSize= 10*mm
GridDimension=200
BesselGridSize = 1*mm #0.2*mm
lambda_=800*nm
innerR_m = 2.333*mm
outerR_m = 2.5333*mm
f=20*cm
dz=0.05*mm
nsteps=400

#Create field
Field=Begin(GridSize,lambda_,GridDimension)
#Create an annulus of light
Field = CircAperture(Field, outerR_m)
Field = CircScreen(Field, innerR_m)

plt.figure()
plt.title('Field before lens')
plt.imshow(Intensity(Field))
plt.show()

#propagate field some distance before focus
dist_z = 190*mm
f1 = GridSize/(GridSize - BesselGridSize)*dist_z
f2 = f1*f/(f1-f) # since 1/f = 1/f1 + 1/f2 (two thin lenses)
Field = Lens(Field,f2) # new syntax, easier and shorter
Field = LensFresnel(Field, f1, dist_z)
Field = Convert(Field) # need to convert back from spherical coords

plt.figure()
plt.title(f'Field after distance {dist_z/mm:.1f}mm')
plt.imshow(Intensity(Field))
plt.show()

#Go through Bessel region and safe Intensity profile
FocusInten=np.zeros([nsteps,Field.N,Field.N])
for a in tqdm(range(0,nsteps)):
    Field=Steps(dz, Field)
    I=Intensity(Field) # flag 0 is default anyway
    FocusInten[a,:,:]=I

plt.figure()
plt.title(f'Field after final step')
plt.imshow(Intensity(Field))
plt.show()

#Normalizing and save
Focusnormalized = FocusInten / FocusInten.max()
FocusInten=(Focusnormalized*65535).astype(np.uint16)
imsave('BesselRegion.tiff',FocusInten)

The calculation of f1 and f2 might seem obscure and maybe does not become clear from the manual [1]. The basic trick is to use 2 lenses, where one is applied as a phase term and the other with the geometrical coordinate transform. If the distance propagated with LensFresnel matched the focal length f, the geometric grid would be infinitely small. With the formula in the code above, the geometric lens will be chosen to have a focal length f1 given by the intercept theorem ("Strahlensatz") and the phase will be applied with a focal length f2. Together both lenses have the focusing strength f.

Hope this helps, greetings from Munich, Lenny

[1] https://opticspy.github.io/lightpipes/manual.html#spherical-coordinates

[ldoyle]

Here is the resulting beam after many steps for 1mm grid size:

And here it is for only 0.2mm grid size after much fewer steps:

PS: My personal experience with the Steps command is that I trust it the least of all propagators. For example, the results change when doing more smaller steps (compared to the same distance in larger steps). So for e.g. measuring exact focal depths, I would cross check with a different method (e.g. using LensFresnel over and over, with a different dist_z every time)

[DerKamin]

Hi Lenny,

thanks a lot for your answer, it really helped me alot! :) Also having the progress bar feels so much better.

I understand the importance of the gridsize much better now. First I had some troubles understanding the calculation for the first of the two lenses. After some time I think i worked it out and combining the two lenses to avoid the problem of grid with the size zero makes sense. I did two quick sketches to illustrate the grid size problem while focussing an annulus of light (First, for anyone who stumbles upon the same problem and second if I am wrong with something i someone can correct me.

Still i have got some questions for this simulation:

1. Would it be better if i avoided the loop which propagates the field with the step command. Instead of the current loop would it be possible / advised to include almost the hole code snippet in a loop and just alter the dist_z parameter within each iteration. If the simulation is done this way would the different calculations of f1 and f2 in each iteration alter the result? @ldoyle Is this what you ment with your last comment?

Why is the the second lens called weak phase mask?

Thanks and greetings, Hannes

[FredvanGoor]

Hi Lenny and Hannes, Nice project! I made a script using MatPlotLib's interactive mode. Now you can see the range over which there is a Bessel beam. You can pause and resume the process by pressing the F10 key (I discovered that on my windows computer, maybe a (secret?) feature of plt.ion() ?)

Hannes, do you like to make a contribution to our documentation-examples when finished? See: https://opticspy.github.io/lightpipes/examples.html#non-diffractive-bessel-beam

Here is my code:

from LightPipes import * 
import matplotlib.pyplot as plt
from tqdm import tqdm # fancy progress bar, recomment installing if you don't have it yet

GridSize= 10*mm
GridDimension=200
BesselGridSize = 1*mm #0.2*mm
lambda_=800*nm
innerR_m = 2.333*mm
outerR_m = 2.5333*mm
f=20*cm
dz=0.15*mm #0.05*mm
nsteps=400

#Create field
Field=Begin(GridSize,lambda_,GridDimension)

#Create an annulus of light
Field = CircAperture(Field, outerR_m)
Field = CircScreen(Field, innerR_m)

#propagate field some distance before focus
dist_z = 190*mm
f1 = GridSize/(GridSize - BesselGridSize)*dist_z
f2 = f1*f/(f1-f) # since 1/f = 1/f1 + 1/f2 (two thin lenses)
Field = Lens(Field,f2) # new syntax, easier and shorter
Field = LensFresnel(Field, f1, dist_z)
Field = Convert(Field) # need to convert back from spherical coords

#enable matplotlib interactive mode
plt.ion()
fig, ax = plt.subplots()

aximg=ax.imshow(Intensity(Field))

#Go through Bessel region, press F10 to pause/resume
for a in tqdm(range(0,nsteps)):
    Field=Steps(dz, Field)
    I=Intensity(Field) # flag 0 is default anyway
    fig.suptitle(f"distance = {a*dz/mm:2.2f}mm")
    aximg.set_data(I)
    fig.canvas.draw()
    fig.canvas.flush_events()

Greetings from, Fred van goor