Comments (5)
ldoyle
commented on June 14, 2024
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
from lightpipes.
yung-p
commented on June 14, 2024
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
from lightpipes.
yung-p
commented on June 14, 2024
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()
from lightpipes.
FredvanGoor
commented on June 14, 2024
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()
from lightpipes.
FredvanGoor
commented on June 14, 2024
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.
from lightpipes.
Related Issues (20)
- Inconsistent definitions of waist w and xshift HOT 2
- Support for astigmatic beams HOT 2
- The grid_dimension and grid_size properties should be documented
- Cannot install LightPipes on Windows 7 - please help HOT 2
- Finite difference documentation typo? HOT 1
- installation problems
- Add Interpol option in Fresnel and Forvard HOT 1
- Warning when paraxiality is violated HOT 2
- Shift with direct integration Forward
- Diffraction From Thin Wire HOT 4
- Step by Step Propagation With a Short Focal Length Lens HOT 3
- LightPipes on an iPad
- LPtest() fails on some systems while installation is correct HOT 1
- Intensity plot gets cut up when attempting to propagate finite airy beam HOT 9
- Cite lightpipes in a paper HOT 2
- lens Following Field Produced by Forward Function
- Airy Beam HOT 6
- LG doughnut mode HOT 4
- Microwave frequencies HOT 3
Recommend Projects
-
React
A declarative, efficient, and flexible JavaScript library for building user interfaces.
-
Vue.js
🖖 Vue.js is a progressive, incrementally-adoptable JavaScript framework for building UI on the web.
-
Typescript
TypeScript is a superset of JavaScript that compiles to clean JavaScript output.
-
TensorFlow
An Open Source Machine Learning Framework for Everyone
-
Django
The Web framework for perfectionists with deadlines.
-
Laravel
A PHP framework for web artisans
-
D3
Bring data to life with SVG, Canvas and HTML. 📊📈🎉
-
Recommend Topics
-
javascript
JavaScript (JS) is a lightweight interpreted programming language with first-class functions.
-
web
Some thing interesting about web. New door for the world.
-
server
A server is a program made to process requests and deliver data to clients.
-
Machine learning
Machine learning is a way of modeling and interpreting data that allows a piece of software to respond intelligently.
-
Visualization
Some thing interesting about visualization, use data art
-
Game
Some thing interesting about game, make everyone happy.
Recommend Org
-
Facebook
We are working to build community through open source technology. NB: members must have two-factor auth.
-
Microsoft
Open source projects and samples from Microsoft.
-
Google
Google ❤️ Open Source for everyone.
-
Alibaba
Alibaba Open Source for everyone
-
D3
Data-Driven Documents codes.
-
Tencent
China tencent open source team.
from lightpipes.