Open yung-p opened 3 years ago
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?
Forvard
/Fresnel
) and the grid dimension N (both with same or larger grid size)I.sum()
)?Hope this helps, Lenny
Hi Lenny,
Thanks for the reply! I think i've made some progress on this.
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
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()
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()
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.
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![4f_fresnel_intensity_Nis1000](https://user-images.githubusercontent.com/71134818/106911324-4c24ae80-6702-11eb-9969-ce0ffbe6c40f.png)
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