Unexpected intensity distribution for Fresnel propagation

[yung-p]

Hello!

This is more of a question about a result i got, rather than an issue i'm experiencing with LightPipes. If this is not the right place to ask this question please let me know!

I am simulating a 4f system and comparing the performance of the Angular Spectrum Method (ASM) and Fresnel propagation.
I have two figures, and in each of the figures in the first row the intensity distributions in three planes are shown; z=0, z=2f and z=4f. The second row show the intensity along the x-axis

My input plane is a unit amplitude plane plane wave masked by a circular aperture with radius r_beam = 3mm.
These are my input parameters:
N = 1000
size = 15 * mm
wave_length = 690 * nm
r_beam = 3 * mm
f1, f2 = 200 * mm, 200 * mm

This is the first figure, using Forvard

And the second figure, using Fresnel

I know that the intensity distributions, in the planes z=0 and z=4f should look pretty much identical for both ASM and Fresnel propagation, but the intensity for Fresnel seems to diminish laterally at z = 800mm. Is this maybe because the conditions i'm in are not paraxial enough?

Thanks in advance!

-Philip

[ldoyle]

Dear Philip,
this is the right place to ask. That is indeed interesting, maybe you can provide a minumum working example of code to reproduce the lower plot?

• From gut feeling, I would say this lens has a long enough focal length compared to beam diameter/ not too fast focusing, I would expect the problem to lie somewhere else (I might be wrong, though).
• You should play around with the grid size (while keeping object size, to investigate edge effects in Forvard/Fresnel) and the grid dimension N (both with same or larger grid size)
• The focal plane lineouts look similar, but one has a factor of 2 higher peak. How well is the integrated intensity conserved (I.sum())?

Hope this helps, Lenny

[yung-p]

Hi Lenny,

Thanks for the reply! I think i've made some progress on this.

• The Fresnel number Nf = r_beam^2/(lambda f) = 65. And I read it must be ~1 in order for the Fresnel approximation to be valid. My intuition also is that a focal length of 200mm was large enough for Fresnel to be valid (a cone with a radius of 3mm and length of 200mm looks pretty paraxial to me). How strictly should one adhere to the Fresnel number?
• Unfortunately, enlarging the grid size did not improve the intensity distribution (and same N) in the 4f plane for Fresnel. However, bigger N (and same grid size), significantly improves the intensity distribution in the 4f plane (shown in the figure below). Also, I'm starting to recognise an Airy disc in the 2f plane.
• The total energy for Fresnel indeed seems to be half of the of Forvard, Here is what I got: E_forv = 125629, E_fr = 64944. Energy should be conserved for both propagation methods right? Is violation of energy conservation also an indication that paraxiality is violated?

N = 2000
size = 15* mm
wave_length = 690 * nm
r_laser_beam = 3 * mm
f1, f2 = 200 * mm, 200 * mm

...
If i choose larger focal lengths, It gets better. The 4f plane is starts to look like more like a unit amplitude wavefront and the Airy disk also starts to become more visible.

N = 2000
size = 15* mm
wave_length = 690 * nm
r_laser_beam = 3 * mm
f1, f2 = 500 * mm, 500 * mm

clearly visible (I overlaid the cutout analytically computed Airy disk in orange).

So, my conclusion would still be that the parameters I initially choose were not paraxial enough for Fresnel to be valid. For Forvard however, the initial parameters do seem appropriate, because the 0, and 4f plane look practically the same. The 2f plane does unfortunately not show an Airy disk. I know this can be solved by a coordinate transform (as was posted in an earlier issue, about one month ago). I have unfortunately not been able to make this work.

Thank you for your suggestions. I'm curious what you think.
-Philip

[yung-p]

And here is the code to produce the plot Fresnel. I somehow wasn't able to link a .py file. I hope this works.

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

N = 1000
size = 15* mm
wave_length = 690 * nm
r_laser_beam = 3 * mm
f1, f2 = 200 * mm, 200 * mm

distances = [0, 2f11000, 4f21000 ]
crop_factor = 5 #zooming in for the middle plot
x = np.linspace(-size / 2, size / 2, N)
y = np.linspace(-size / 2, size / 2, N)

def Fresnel_sim_4f_system(N, size, wave_length, r_laser_beam, f1, f2):
focal_lengths = np.array((f1, f2))
d_s = [focal_lengths[0], focal_lengths[0] + focal_lengths[1], focal_lengths[1]]

n_d1 = 2
n_d2 = 3
n_d3 = 2

d_1 = np.linspace(0, d_s[0], n_d1)
d_2 = np.linspace(0, d_s[1], n_d2)
d_3 = np.linspace(0, d_s[2], n_d3)


# preallocation of fields
fields_1 = np.zeros((N, N, n_d1)).astype(complex)
fields_2 = np.zeros((N, N, n_d2)).astype(complex)
fields_3 = np.zeros((N, N, n_d3)).astype(complex)

F_0 = Begin(size, wave_length, N)  # initial field
F_0 = CircAperture(F_0, r_laser_beam, x_shift=0.0, y_shift=0.0)

for i in range(0, n_d1):
    fields_1[:, :, i] = Fresnel(F_0, d_1[i]).field

P_1 = Fresnel(F_0, d_1[-1])
L_1 = Lens(P_1, f1, x_shift=0.0, y_shift=0.0)  # second normal lens

for i in range(0, n_d2):
    fields_2[:, :, i] = Fresnel(L_1, d_2[i]).field

P_2 = Fresnel(L_1,d_2[-1])
L_2 = Lens(P_2, f2, x_shift=0.0, y_shift=0.0)

for i in range(0, n_d3):
    fields_3[:, :, i] = Fresnel(L_2, d_3[i]).field

intensity_1 = abs(fields_1) ** 2
intensity_2 = abs(fields_2) ** 2
intensity_3 = abs(fields_3) ** 2

phases_1 = np.angle(fields_1)
phases_2 = np.angle(fields_2)
phases_3 = np.angle(fields_3)

two_f_energy = print(Intensity(Fresnel(L_1, d_2[1])).sum())

return intensity_1, intensity_2, intensity_3, phases_1, phases_2, phases_3,two_f_energy

intensity_1,
intensity_2,
intensity_3, phases_1, phases_2, phases_3, two_f_energy = Fresnel_sim_4f_system(N, size, wave_length, r_laser_beam, f1, f2)

def intenstiy_phase_cropper(intensity, phase, N,crop_factor): #function to crop middle plot
cropped_intensity_2 = intensity[:,:,1][int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor)), int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor))]
cropped_phase_2 = phase[int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor)),int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor))]

new_x = np.linspace(-((size/2) / crop_factor), ((size/2) / crop_factor), int(N / crop_factor) )*1000
new_y = np.linspace(-((size/2) / crop_factor), ((size/2) / crop_factor), int(N / crop_factor) )*1000
return cropped_intensity_2, cropped_phase_2, new_x, new_y

cropped_intensity_2,
cropped_phase_2,
new_x,
new_y = intenstiy_phase_cropper(intensity_2, phases_2[:,:,1], N, crop_factor)

fig1, axs = plt.subplots(2, 3, sharex='col', gridspec_kw={'hspace': 0.05}, figsize=(10, 7))
axs[0, 0].pcolormesh(x, y, intensity_1[:, :, 0], shading='auto',
vmin=-1, vmax=1, cmap='twilight')
axs[0, 1].pcolormesh(new_x, new_y, cropped_intensity_2, shading='auto',
vmin=-1, vmax=1, cmap='twilight')
axs[0, 2].pcolormesh(x, y, intensity_3[:, :, -1], shading='auto',
vmin=-1, vmax=1, cmap='twilight')
axs[1, 0].plot(x, intensity_1[:, :, 0][int(N / 2), :])
axs[1, 1].plot(new_x, (np.amax(cropped_intensity_2[int(N / crop_factor / 2), :])**-1)*cropped_intensity_2[int(N / crop_factor / 2), :])

axs[1, 2].plot(x, intensity_3[:, :, -1][int(N / 2), :])
axs[0, 0].set_ylabel('Intensity Patterns \n y[mm]', fontweight='bold')

for j in range(0, 3):
axs[0, j].set_title(str(distances[j]) + 'mm', fontsize=11)
axs[1, j].set_xlabel('x [mm]', fontweight='bold');
axs[1, 0].yaxis.set_visible(True)
axs[1, 0].set_ylabel('Amplitude [arb.]', fontweight='bold')

plt.show()

[FredvanGoor]

Hi yung-p,

I tried your script.
Put your code between two ``` lines.

You should use Forvard if the Fresnel number >> 1 (angular spectrum method). This is the case in your problem:
N_Fresnel = (r_laser_beam)^2/(wave_length * f1) = 65.22

Fred van Goor.

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

N = 1000
size = 15* mm
wave_length = 690 * nm
r_laser_beam = 3 * mm
f1, f2 = 200 * mm, 200 * mm

distances = [0, 2*f1*1000, 4*f2*1000 ]
crop_factor = 5 #zooming in for the middle plot
x = np.linspace(-size / 2, size / 2, N)
y = np.linspace(-size / 2, size / 2, N)

def Fresnel_sim_4f_system(N, size, wave_length, r_laser_beam, f1, f2):
    focal_lengths = np.array((f1, f2))
    d_s = [focal_lengths[0], focal_lengths[0] + focal_lengths[1], focal_lengths[1]]
    
    n_d1 = 2
    n_d2 = 3
    n_d3 = 2
    
    d_1 = np.linspace(0, d_s[0], n_d1)
    d_2 = np.linspace(0, d_s[1], n_d2)
    d_3 = np.linspace(0, d_s[2], n_d3)
    
    
    # preallocation of fields
    fields_1 = np.zeros((N, N, n_d1)).astype(complex)
    fields_2 = np.zeros((N, N, n_d2)).astype(complex)
    fields_3 = np.zeros((N, N, n_d3)).astype(complex)
    
    F_0 = Begin(size, wave_length, N)  # initial field
    F_0 = CircAperture(F_0, r_laser_beam, x_shift=0.0, y_shift=0.0)
    
    for i in range(0, n_d1):
        fields_1[:, :, i] = Fresnel(F_0, d_1[i]).field
    
    P_1 = Fresnel(F_0, d_1[-1])
    L_1 = Lens(P_1, f1, x_shift=0.0, y_shift=0.0)  # second normal lens
    
    for i in range(0, n_d2):
        fields_2[:, :, i] = Fresnel(L_1, d_2[i]).field
    
    P_2 = Fresnel(L_1,d_2[-1])
    L_2 = Lens(P_2, f2, x_shift=0.0, y_shift=0.0)
    
    for i in range(0, n_d3):
        fields_3[:, :, i] = Fresnel(L_2, d_3[i]).field
    
    intensity_1 = abs(fields_1) ** 2
    intensity_2 = abs(fields_2) ** 2
    intensity_3 = abs(fields_3) ** 2
    
    phases_1 = np.angle(fields_1)
    phases_2 = np.angle(fields_2)
    phases_3 = np.angle(fields_3)
    
    two_f_energy = print(Intensity(Fresnel(L_1, d_2[1])).sum())
    
    return intensity_1, intensity_2, intensity_3, phases_1, phases_2, phases_3,two_f_energy



intensity_1,intensity_2,intensity_3, phases_1, phases_2, phases_3, two_f_energy = Fresnel_sim_4f_system(N, size, wave_length, r_laser_beam, f1, f2)

def intenstiy_phase_cropper(intensity, phase, N,crop_factor): #function to crop middle plot
    cropped_intensity_2 = intensity[:,:,1][int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor)), int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor))]
    cropped_phase_2 = phase[int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor)),int(0.5*(N-N/crop_factor)):-int(0.5*(N-N/crop_factor))]
    
    new_x = np.linspace(-((size/2) / crop_factor), ((size/2) / crop_factor), int(N / crop_factor) )*1000
    new_y = np.linspace(-((size/2) / crop_factor), ((size/2) / crop_factor), int(N / crop_factor) )*1000
    return cropped_intensity_2, cropped_phase_2, new_x, new_y

cropped_intensity_2,cropped_phase_2,new_x,new_y = intenstiy_phase_cropper(intensity_2, phases_2[:,:,1], N, crop_factor)

fig1, axs = plt.subplots(2, 3, sharex='col', gridspec_kw={'hspace': 0.05}, figsize=(10, 7))
axs[0, 0].pcolormesh(x, y, intensity_1[:, :, 0], shading='auto',
vmin=-1, vmax=1, cmap='twilight')
axs[0, 1].pcolormesh(new_x, new_y, cropped_intensity_2, shading='auto',
vmin=-1, vmax=1, cmap='twilight')
axs[0, 2].pcolormesh(x, y, intensity_3[:, :, -1], shading='auto',
vmin=-1, vmax=1, cmap='twilight')
axs[1, 0].plot(x, intensity_1[:, :, 0][int(N / 2), :])
axs[1, 1].plot(new_x, (np.amax(cropped_intensity_2[int(N / crop_factor / 2), :])**-1)*cropped_intensity_2[int(N / crop_factor / 2), :])

axs[1, 2].plot(x, intensity_3[:, :, -1][int(N / 2), :])
axs[0, 0].set_ylabel('Intensity Patterns \n y[mm]', fontweight='bold')

for j in range(0, 3):
    axs[0, j].set_title(str(distances[j]) + 'mm', fontsize=11)
    axs[1, j].set_xlabel('x [mm]', fontweight='bold');
    axs[1, 0].yaxis.set_visible(True)
    axs[1, 0].set_ylabel('Amplitude [arb.]', fontweight='bold')

plt.show()

Alternative code that does the same, but using Forvard.
Increase the grid for better results.
Use Forvard with high Fresnel number.

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

N = 3000
N2=int(N/2)
size = 15* mm
wavelength = 690 * nm
r_aperture = 3 * mm
f1, f2 = 2000 * mm, 2000 * mm

F=Begin(size,wavelength,N)
F=CircAperture(F,r_aperture)

print(r_aperture*r_aperture/wavelength/f1) #Fresnel number

I0=Intensity(F)
F=Forvard(F,f1)
F=Lens(F,f1)
F=Forvard(F,f1)
I1=Intensity(F)
F=Forvard(F,f2)
F=Lens(F,f2)
F=Forvard(F,f2)
I2=Intensity(F)

fig=plt.figure(figsize=(11,6))
fig.suptitle("4f relay imaging")
ax1 = fig.add_subplot(231);ax1.axis('off')
ax2 = fig.add_subplot(232);ax2.axis('off')
ax3 = fig.add_subplot(233);ax3.axis('off')
ax4 = fig.add_subplot(234)
ax5 = fig.add_subplot(235)
ax6 = fig.add_subplot(236)

ax1.imshow(I0,cmap='jet')

ax2.imshow(I1,cmap='jet')

ax2.set_xlim(N2-100,N2+100)
ax2.set_ylim(N2-100,N2+100)

ax3.imshow(I0,cmap='jet')

X=np.linspace(-size/2,size/2,N)
ax4.plot(X/mm,I0[N2]); ax4.set_xlabel('x[mm]')

ax5.plot(X/mm,I1[N2]); ax5.set_xlabel('x[mm]')
ax5.set_xlim(-0.8, 0.8)

ax6.plot(X/mm,I2[N2]); ax6.set_xlabel('x[mm]')

plt.show()

[FredvanGoor]

Maybe we should add a new propagation command (let's call it Propagate(F,z) or something like that) It calculates the Fresnel number and decides to use Forvard or Fresnel. Another method to decide is to use the so-called Gauss pilot beam method. First a Gauss beam propagates through the system and based on the radius of curvature of the beam the propagation method is chosen.
We could also introduce a propagation command using ABCD matrices and pure Gauss beams. I have tried that in the past for the Mathcad (and or?) Matlab versions of LightPipes. It works fine.