Bessel Beam Simulation from an annulus of light about lightpipes HOT 4 OPEN

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)

from lightpipes.

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

from lightpipes.

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

from lightpipes.

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

from lightpipes.