Propagation of Gaussian Beam through free space and through lens

[IvanOstr]

Hey there, 1) I'm trying to propagate a gaussian beam in free space using the code below. I've tried two method - once with Forvard/Fresnel and once with Steps. In both cases the gaussian beam expands, unlike in a laser system where the gaussian beam stays almost constant because of collimation and a really large Rayleigh length. I want to simulate the case with a really large gaussian beam, and I did not find in your code the full gaussian beam as in Saleh or other references: Should I implement it by myself and put it in the field matrix? Or does the beam expansion is just a result of the FFT approach?

2) I think the both methods does not agree completely. What are the reasons, and is there a guide line when to use each? for example, in steps we see multiple peaks at some large z:

3) When I propagate through lens with this methods, I'm afraid i'm not getting the right results as of the reason stated in 1.

Thanks in advance, Ivan

The code:

from LightPipes import*
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

wavelength = 632.8*nm
size = 2*mm
N = 200
N2 = int(N/2)
R = 0.3*mm
dz = 10*mm
f = 15*cm
n = (1.0 + 0.1j)*np.ones((N,N))
Z_len = 200
Icross = np.zeros((Z_len,N))
X = range(N)
Z = range(Z_len)
X, Z = np.meshgrid(X, Z)

F = Begin(size, wavelength, N)
F = GaussBeam(R, F)

############### Forvad / Fresnel ################
I=Intensity(F, 0)
plt.figure()
plt.imshow(I)
plt.show()

#F=Lens(F,f)

zs = np.linspace(0, 30, 10)
for z in zs:
    F_propogated = Fresnel(F, z*cm)
    # F_propogated = Forvard(F, z*cm)
    I=Intensity(F_propogated, 0)
    plt.figure()
    plt.plot(Intensity(F_propogated, 0)[N2, :])
    # plt.imshow(I)
    plt.show()

############### Steps ################
for i in range(0, Z_len):
  F=Steps(F, dz, 1, n)
  I=Intensity(F, 0)
  Icross[i][:N] = I[N2][:N]
  # plt.figure();plt.plot(Icross[i][:N]);plt.show()

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Z,
        Icross,
        rstride=1,
        cstride=1,
        cmap='rainbow',
        linewidth=0.0,
        )
plt.axis('off')
plt.show()

[ldoyle]

Dear Ivan,

I think you are on the right track. I'll try to comment each point:

1. The propagation with Fresnel/Forvard seems fine. For wavelength=632.8nm and w0=0.3mm the Rayleigh range is z_R=446.81mm according to an online Gauss beam calculator, so for your 30cm-200cm propagation indeed becomes significant.
2. The Steps result indeed looks weird after a longer propagation of ~1m. However I believe this is due to boundary effects, so this should be solveable by using a larger grid size (requiring more grid points if you want to keep the resolution, see comment below).
3. As long as the boundary problems etc. are solved, using Lens should work as expected. In general I recommend reading the manual and playing around with the examples to learn about the pitfalls of the calculations (e.g. edge effects or too hard focusing/ too short focal length compared to beam size)
4. You are right, it seems the full Gauss beam formula is not implemented. However the GaussBeam, GaussHermite and GaussLaguerre commands should be sufficient for most cases, as they provide the correct formula at beam waist (z=0). So for e.g. z=1mm, you can just use a propagation with Fresnel instead of using the analytical formula. Nevertheless, if you think it could be useful, please feel free to play around and submit a Pull Request if you're happy with your implementation.

I have added some printouts to your code and a final plot showing the weird intensity distribution after the final Steps propagation.

• At z~47cm, indeed the peak intensity drops to ~50% both for Fresnel and Steps, confirming the expected Rayleigh range.
• The intensity after the last Step shows some "square" behaviour, always an indicator you might be running into edge effects
• If you increase the size to size=8*mm and N to N=800 while keeping the beam size fixed, you can see the much better result after the final Step.

Let us know if this helps. All the best, Lenny

from LightPipes import *
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt

wavelength = 632.8*nm
size = 2*mm
N = 200
N2 = int(N/2)
R = 0.3*mm # if R means w0 -> at 633nm, z_Rayleigh is 446.67mm ~ 44.7cm
dz = 10*mm
f = 15*cm
n = (1.0 + 0.1j)*np.ones((N,N))
Z_len = 200
Icross = np.zeros((Z_len,N))
X = range(N)
Z = range(Z_len)
X, Z = np.meshgrid(X, Z)

F = Begin(size, wavelength, N)
F = GaussBeam(R, F)

############### Forvad / Fresnel ################
I=Intensity(F, 0)
print('Imax at z=0.0cm', I.max())
plt.figure()
plt.imshow(I)
plt.show()

# raise Exception('eof')

#F=Lens(F,f)

zs = np.linspace(0, 44.7, 10)
for z in zs:
    F_propogated = Fresnel(F, z*cm)
    # F_propogated = Forvard(F, z*cm)
    I=Intensity(F_propogated, 0)
    print(f'Imax at z={z:.1f}cm (Fresnel)', I.max())
    plt.figure()
    plt.plot(Intensity(F_propogated, 0)[N2, :])
    # plt.imshow(I)
    plt.show()

############### Steps ################
z_steps = 0.0
for i in range(0, Z_len):
  F=Steps(F, dz, 1, n)
  z_steps += dz
  I=Intensity(F, 0)
  print(f'Imax at z={z_steps/cm:.1f}cm (Steps)', I.max())
  Icross[i][:N] = I[N2][:N]
  # plt.figure();plt.plot(Icross[i][:N]);plt.show()

fig = plt.figure()
plt.imshow(I)
plt.title('Intensity after final Steps() call')

fig = plt.figure()
ax = fig.gca(projection='3d')
ax.plot_surface(X, Z,
        Icross,
        rstride=1,
        cstride=1,
        cmap='rainbow',
        linewidth=0.0,
        )
plt.axis('off')
plt.title('Lineout of intensity on axis along Steps()')
plt.show()

[IvanOstr]

Dear Lenny, Thank you very much for the comprehensive response!

1. , 2. , code fixes: You are right, that is indeed the case. I will definitely play with the code in the manual and add the full gaussian version if I will need to code it.

I'll update about successful propagation through lens with the same setup in the few days, if you don't close the issue by then.

[IvanOstr]

I have found out some of the effects you mentioned regarding too short or too hard focusing. 1) A strange behavior starts when the beam parameters are w0 = 10mm, wavelength = 800nm, focal length = 150cm, grid_dim = 1000 and continuous along all the propagation length. For example: This is taken from the code below. The situation gets worse with shorter wavelengths. However, I found out that increasing the grid dimension can help that a lot. What is the reason for that? I tried to investigate it and found out that by making the lens smaller - the strange effects seem to vanish more. For example, by manually setting the phase to be the same phase as before but only on a smaller square on the middle. Suppose I want to simulate a real life experiment with focal length = 15 cm, what would you recommend doing?

2) It seems that Steps does not agree with Fresnel when I change focal length to 250: The left picture is Fresnel and the right one is the same length in z, but with the Steps function. I guess it is a drawback of the method somehow? Increasing the grid size didn't work so well here. The steps code I used (written straight after the forvard code) is attached below.

3) Suppose I try to simulate a really small lens - say, 10 to 100 times smaller than the gaussian beam. Do you know of any limitations of the simulations when trying to simulate diffraction and refraction from really small objects? My method is currently to cut a small square of phase and apply it on the field.

Code for (2):

Z_len = 30# mm
Icross = np.zeros((Z_len, grid_dimension))
X = range(grid_dimension)
Z = range(Z_len)
X, Z = np.meshgrid(X, Z)
dz = 100*mm
n = (1.0 + 0j)*np.ones((grid_dimension, grid_dimension))
for i in range(0, Z_len):
  print(i)
  F=Steps(F, dz, 1, n)
  I=Intensity(F, 0)
  Icross[i][:grid_dimension] = I[N2][:grid_dimension]
  # plt.figure();plt.plot(Icross[i][:N]);plt.show()

plt.figure();plt.plot(x_arr,Icross[25][:grid_dimension]);plt.show()

Code for (1):

#Focus of a lens.
from LightPipes import*
import numpy as np
from mpl_toolkits.mplot3d import Axes3D
import matplotlib.pyplot as plt
from scipy.optimize import curve_fit

def gaussian(x, a, b, c):
    return a*np.exp(- ((  (x-b)/c )**2) / 2)

wavelength = 800 * nm
k = 2 * np.pi / wavelength
grid_dimension = 1000
N2 = int(grid_dimension / 2)
waist_in_mm = 10 # waist at z = 0 in mm
W = waist_in_mm * mm
size = W * 5 # mm
x_arr = np.linspace(0, waist_in_mm*5, grid_dimension) # mm

F = Begin(size, wavelength, grid_dimension)
F = GaussBeam(W, F)

focal_length = 150
f = focal_length * cm
F = Lens(F, f)

############### Forvad / Fresnel ################
I = Intensity(F, 0)
plt.figure()
plt.imshow(I)
plt.title('Intial Distribution')
plt.show()

zs = np.linspace(0, focal_length, 5)
for z in zs:
    F_propogated = Fresnel(F, z*cm)
    #F_propogated = Forvard(F, z*cm)
    I=Intensity(F_propogated, 0)
    plt.figure()
    plt.title('Propagation After ' + str(z) +'cm')
    plt.plot(x_arr, I[N2, :])
    # p0 = (1, 25, 1)
    # popt, pcov = curve_fit(gaussian, x_arr, I[N2, :], p0)
    # plt.plot(x_arr, gaussian(x_arr, *popt), 'r-')
    # plt.title('Propagation after ' + str(z) +'cm with waist = ' + str(2*popt[2]))
    # plt.imshow(I)
    plt.show()

[ldoyle]

Dear Ivan, I finally found the time to test your updated examples (by the way thanks for supplying a fully working code/ MWE, it helps a lot to investigate).

1. I do not have a complete or even scientific explanation, but I am sure it is related to the insufficient sampling of the field. You can plot the phase immediately after the lens (Warning: PhaseUnwrap or Phase(F, unwrap=True) take a while unfortunately). If you observe the wrapped phase (unwrap=False) you will notice the step from one grid point to the next will be delta_phi<< 2pi close to the center, but rapidly varying at the edges. Using unwrap=True the routine will try to "guess" the correct offset of 2pi to make the phase look smoother. This starts to fail at the edges (the square shape is an artefact of the unwrap routine).

-> In my experience, once the phase delta to neighbouring grid points becomes larger than 2pi (in regions where the field has non-negligible intensity), things start to get weird. In your case, Forvard seems to cope a little better than Fresnel (f=250cm vs. f=550cm before weird effects show up, resp.). However try both, in other cases it might be the other way around.

2. As you can see, you have very few grid points to sample the focus, so changes in the shape are hard to interpret. Very roughly, the focus spot size will be r0 = lambda*f/D. For your lens with f=150cm, D=2cm (est. from Gauss waist) and lambda=0.8um the focus spot radius will be ~60um, but your grid spacing is 50mm/1000=50um. To solve this and the above question about the real life experiment I recommend the LightPipes commands LensFresnel/LensForvard. Check the documentation for a simple example, I can elaborate on the f1/f2 trick used if necessary (I have been meaning to improve the documentation here anyway at some point). https://opticspy.github.io/lightpipes/manual.html#spherical-coordinates This effectively allows you to have a different grid spacing in the focus than at the lens. It also allows to use much shorter focal lengths compared to the beam diameter, however I have never investigated what the new limit is yet (I have joined LightPipes later in the game, the original implementation seems exactly as old as I am...) Note that Steps only works in normal coordinates, not spherical ones.

3. I'm not sure I follow your last question. Concretely, how large is the beam, the lens, and the focal length? If the beam is much larger than the lens, what do you do with the parts that lie outside, do you apply an aperture? Very generally, the underlying physics/maths use various methods to solve the diffractions integrals on a grid, so your grid spacing must be chosen to sufficiently account for the objects size/shape and potential phase shifts (see above). If you want to be sure the result makes sense you should start with a finer spacing and then decrease the sampling. If the result does not change much, you are fine. I would start by resolving your tiny objects by at least 10 grid points, e.g. a circular aperture should not be a cross of 3x3pixels but rather a smooth circle with 10x10pixels/gridpoints outer diameter.

Code for (1):

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

wavelength = 800 * nm
grid_dimension = 1000
N2 = int(grid_dimension / 2)
waist_in_mm = 10 # waist at z = 0 in mm
W = waist_in_mm * mm
size = W * 5 # mm
x_arr = np.linspace(0, waist_in_mm*5, grid_dimension) # mm

F = Begin(size, wavelength, grid_dimension)
F = GaussBeam(W, F)

# focal_length = 200 # effects start to show in Forvard
# focal_length = 250 # no effects in Forvard
# focal_length = 450 # effects start to show in Fresnel
focal_length = 550 # (almost) no effects in Fresnel
f = focal_length * cm
F = Lens(F, f)

############### Forvard / Fresnel ################
I = Intensity(F, 0)
do_unwrap = False
Ph = Phase(F, unwrap=do_unwrap)
PhLineout = Ph[N2, :]

fig, [axl, axm, axr] = plt.subplots(1, 3)
axl.imshow(I)
axm.imshow(Ph)
if do_unwrap:
    axr.plot(PhLineout)
else:
    axr.plot(PhLineout, 'x-') #crosses make jumps more visible
axl.set_title('Intensity at z=0')
axm.set_title('Phase at z=0')
axr.set_title('Phase lineout at z=0')
plt.show()

zs = np.linspace(0, focal_length, 5)
for z in zs:
    F_propogated = Fresnel(F, z*cm)
    # F_propogated = Forvard(F, z*cm)
    I=Intensity(F_propogated, 0)
    plt.figure()
    plt.title('Propagation After ' + str(z) +'cm')
    plt.plot(x_arr, I[N2, :])
    plt.show()

[jmmelko]

I agree that support for a general Gaussian beam would be useful. Should not be difficult to implement: just add parameters and compute the right intensity/phase functions. In the very least different waist sizes should be possible.

[FredvanGoor]

I propose new commands for propagating Gauss beams using analytical ABCD theory: user defined ABCD matrix, free space propagation and a lens command. I have add them to a new branch "develop_Fred" in this repo. I called them GaussABCD, GaussForvard and GaussLens, but I think that another name is better: maybe GBForvard and GBlens or GForvard and GLens. Comments are welcome! Please look at the test scripts in this branch.

[IvanOstr]

Dear Leonard: Thank you very much for the answer again. I hope the questions I ask are not too much :)

1. Interesting, I agree it has to do with the spacing. It is said in the manual that Fresnel should be more reliable. Maybe that's not the case in specific examples.
2. It is also a good explanation. I tried reading the part of "Spherical Coordinates" in the manual thoroughly and got quite confused. It seems to me like magic currently. Here is what I don't understand: a) I'm not sure what is meant by "The spherical coordinates follow the geometrical section of the divergent or convergent light beam". What is the geometrical section exactly? Is it the beam shape when looking "from above" (from the z/optical path direction). Also, the beam may not be divergent or convergent. b) How is following the beam in a floating coordinate system is implemented mathematically? I tried looking into the code of LensFresnel. I'm not sure I catch what's going on. I see the focal length and z which are passed as parameters are transformed by the lens makers law, but I'm not sure how it is equivalent to transforming a coordinate system and the other parts of the code. c) If I understand correctly, the intuition behind decreasing the grid spacing when transforming cartesian to spherical coordinate system is the same as if we take a flat paper and curve it with the hands. The paper now has less "area" if we look from above and thus the grid spacing changing to smaller one, although not evenly sampled anymore. Is that the direction of thought? d) Regarding the f1/f2 trick - please elaborate. From what I understand - the focal length f is created by combining two lenses - one real lens with focal length f1 and one lens (which is not really a lens but stems from the fact that you transform the coordinate system to spherical?) with focal length f2 and using the lens maker formula. This way your grid is not zero at the focal length of LensFresnel. Even if I did understand that, I'm not sure how Lens Fresnel is actually a lens (same point as b). e) The example looks really cool, do you have a hint why the phase becomes so flat in the focal length?
3. I'll upload an example later. Generally, I'm trying to understand if the diffraction by a really small lens is phenomenon that stems more from the fact that it is a small object that bends the wave front or from the accumulated phase when going through this object.

Dear Fredvan: Good idea! The problem i'll have with that is that the aperture (lens) I use may be much smaller then the beam, but the functions will be really helpful otherwise.

[photon-rider]

Dear Leonard,

Greetings..

Currently I am trying to realize a 4f-system, with LightPipes, and seem to have run into similar issues as discussed in this thread.

1. As mentioned in your reply above could you please elaborate on the trick of using f1 & f2 w.r.t f, while applying the LensFresnel command. I have followed this link https://opticspy.github.io/lightpipes/manual.html#spherical-coordinates and modified the example presented there to compare the transverse magnification (i.e. between my object(an aperture) and the corresponding image).

The modified code: (While running this piece I encounter the following error, which seems to be related to the values one assigns to f1 "LensFresnel: Behind focus")

from LightPipes import *
import numpy as np
import matplotlib.pyplot as plt
#(1) BASICS
N = 3001

labda=800*nm
size=50*cm

u=100*cm                    #Object Distance from the Lens
v=100*cm                    #Image Distance from the Lens
f=50*cm                       #Focal Length
f1=10*m                       # The tricks..
f2=f1*f/(f1-f)

frac=f/f1
newsize=frac*size

w=6*cm

#(2) OPTICAL SYSTEM
F=Begin(size,labda,N);
F=CircAperture(w,0,0,F);
ph_in=Phase(F);ph_in=PhaseUnwrap(ph_in)
I_in=Intensity(F,0)

F=Fresnel(u,F)                    #Propogate Object to Lens

F=Lens(f1,0,0,F)                 #ApplyLens
F=LensFresnel(f2,v,F)
F=Convert(F)
ph_out=Phase(F);ph_out=PhaseUnwrap(ph_out)
I_out=Intensity(F,0)

x=[]
for i in range(N):
    x.append((-newsize/2+i*newsize/N)/cm)

#(3) THE PLOTS
fig = plt.figure(figsize=(20,18))
ax1 = plt.subplot(311)
ax2 = plt.subplot(312)

ax1.plot(x,ph_in[int(N/2)],'r--',label='At Object')
ax1.plot(x,ph_out[int(N/2)],'b',label='At Image');
#ax1.set_xlim(-newsize/2/cm,newsize/2/cm)
#ax1.set_ylim(-2,4)
ax1.set_xlabel('x [cm]');
ax1.set_ylabel('phase [rad]');
ax1.set_title('Phase')
legend = ax1.legend(loc='upper center', shadow=True)

ax2.plot(x,I_in[int(N/2)],'r--',label='At Object')
ax2.plot(x,I_out[int(N/2)], 'b',label='At Image');
#ax2.set_xlim(-newsize/2/cm,newsize/2/cm)
#ax2.set_ylim(0,1000)
ax2.set_xlabel('x [cm]');
ax2.set_ylabel('Intensity [a.u.]');
ax2.set_title('Intensity')
legend = ax2.legend(loc='upper center', shadow=True)

plt.show()

2. Could you also comment upon the limitations with these methods (gridsize, propagation distances, etc..) in relation to the parameters in my task; my beam width is large (~12cm), and I am placing an obstacle ~_1cm_wide along its path. The task is to magnify the object by 10x or 100x (which is also large).

3. Am I correct in thinking it would be a mistake to not leave sufficient padding (of zeros) between my object and the simulation boundaries (i.e. size) , hence I probably cannot magnify arbitrarily.

4. Lastly, is there a simpler way to work this if I can ignore the distribution focal plane (after the lens).

I am new to LightPipes, and hope that it will be possible to find some guidance here :)

Regards -Hari