I am simulating the calculation of a circular aperture beam passing through a lens and converging at the focal plane. I want to calculate what is the radius of the spot when the energy in the barrel at the focal plane is 83.9%. And calculate its times of the diffraction limit.
If the diameter of the center of the first dark ring (Airy disk) of a Fraunhofer diffraction pattern is taken as the diffraction limit diameter.Theoretically, When focused through an ideal lens to the focal plane, the measured ideal diffraction limit diameter should be $d{difflimit} =\theta {difflimit} f$ = $2.44 λ * f / D$, where λ is the wavelength, f is the focal length, and D is the diameter of the circular aperture.[1]
According to Formula 1, it can be calculated that $u{ref}$ =0.839,when $d{u,\mathrm{ref}}$ is taken as diffraction limit diameter metioned above
There are some details for $u_{ref}$
According to Fraunhofer diffraction theory, the far-field intensity distribution of the reference beam with radius $r$ is calculated in paraximal approximation see A.1
$$I\left(r\right)=I_0\left\lbrace\frac{2J_{1}\left[\pi d_{u,near}r/\left(f\lambda\right)\right]}{\pi d_{u,near}r/\left(f\lambda\right)}-\varepsilon^2\frac{2J_{1}\left[\pi\varepsilon d_{u,near}r/\left(f\lambda\right)\right]}{\pi\varepsilon d_{u,near}r/\left(f\lambda\right)}\right\rbrace^2\cdots\cdots\cdots\cdots(\text{ A.}1)$$
$I_{0}$ -- the reference beam far-field peak intensity;
$J_{1}$ -- first-order Bessel function;
$d_\mathrm{u.near}$ -- reference beam near-field beam diameter;
$f$ -- device equivalent focal length;
$\lambda$ -- actual beam laser wavelength;
$\varepsilon$ - center blocking ratio. (center circle area as a proportion of total area)
Let $I(r)=0$, the positional equation to be satisfied by the zero radius of the circular beam intensity is given in Eq. (A.2).
$J_1[\pi d_\text{u.near}r/(f\lambda)]-\varepsilon J_1[\pi\varepsilon d_\text{u,near}r/(f\lambda)]=0\quad........................( A.2 )$
According to Formula 1, when $\varepsilon$ = 0, it can be calculated that c=2.44, $u_{ref}$ =0.839
So the calculation result we would expect should be $\beta = d{simulate}/d{therory} =1 $,where $d{simulate}$ is the circular diameter of the region that accounts for 83.9% of the total energy in the simulation results，and $d{therory}$ = $2.44 F.lam f / w$
But in practice, the code didn't work as I thought it would.
Code Parts
At first, I tried combining the Lens() and Propagate() functions.
Then ,I tried LensFarfield() Fuction
And now,I am trying use Spherical coordinates
calRediuswithPow.py is the base function I used to calculate it, there are some things wrong with its algorithm, such as the center of mass calculating itself, and the energy calculating itself, but I don't think it's going to be particularly wrong for the same field and it's propagation, and during debugging, the energy percentage was not more than 0.3% wrong
They all produce beta results that are way off compared to one, generally speaking aspherical coordinate calculations will be about 30% off, and will change greatly with focal length, and diameter of the circular hole,even if I use "ideal sampling": $dx=\frac{\lambda *z}{L}$ [2]73-74. Spherical coordinates in my example are calculated to be 0.3-0.4, which is off by 60-70%.
I'm not sure what went wrong, the radius calculations should not logically produce such a large deviation. Is there something I am not taking into account in my simulation method.
calRediuswithPow.py
### use calRediuswithPow() function to calculate the radius corresponding to energy in specify bucket
```python
from LightPipes import *
import matplotlib as mpl
import numpy as np
import copy as _copy
def calRediuswithPow(fieldIn,energyPercent=0.839):
'''
Calculate the radius corresponding to the ring energy
Parameters
----------
energyPercent : double
fieldIn : lp.Field
The input field.
Returns
-------
The Redius correspend to the energyPercent power.
"""
----------
'''
fieldIn=_copy.copy(fieldIn)
I=Intensity(fieldIn)
sumEnergy=0
# sumEnergy=sum(I)
sumEnergy1=0
sumEnergy1=0
F_row=np.size(fieldIn.field,0)
F_col=np.size(fieldIn.field,1)
for i in range(F_row):
for j in range(F_col):
sumEnergy+=I[i][j]
#sumEnergy2=Power(fieldIn)
targetEnergy=energyPercent*sumEnergy
# Y, X = fieldIn.mgrid_cartesian
#Determine spot center of mass Coordinate system with array indexed x,y coordinates
x_power_sum=0
y_power_sum=0
for i in range(F_row):
for j in range(F_col):
x_power_sum+=I[i][j]*i
for i in range(F_row):
for j in range(F_col):
y_power_sum+=I[i][j]*j
x_center=x_power_sum/sumEnergy
y_center=y_power_sum/sumEnergy
#Calculate the approximation of the true radius using the Bisection method.
# First set the maximum value of the radius to the diagonal of the matrix
r_max=np.sqrt(F_row**2+F_col**2)
r_min=0
girdR=binarySearchR(I,r_max,r_min,targetEnergy,x_center,y_center,F_row,F_col)
Rres=girdR*fieldIn.dx
return Rres
def binarySearchR(intensityIn,r_max,r_min,targetEnergy,x_center,y_center,F_row,F_col):
'''
Bisection method to approximate the real radius (the fitting effect is a bit poor, later can be changed to a point by point cascade outward expansion)
The error is close to 0.3%
'''
while r_min <= r_max:
r_mid=(r_min+(r_max-r_min)/2)
midEnergy=calPowerWithRadius(intensityIn,F_row,F_col,x_center,y_center,r_mid)
if midEnergy == targetEnergy:
return r_mid
elif midEnergy>targetEnergy:
r_max=r_mid-1
else:
r_min=r_mid+1
#close
return r_mid
#TODO:Sort, compute distance from center of mass for all cells, put into set sort, set map
def calPowerWithRadius(intensityIn,F_row,F_col,x_center,y_center,radius):
'''
Calculate the power within a specified radius
'''
powerRes=0
for i in range(F_row):
for j in range(F_col):
if ((i-x_center)**2+(j-y_center)**2)<=radius**2:
powerRes+=intensityIn[i][j]
return powerRes
```
LensAndLensFarfieldTry.py
### try combine Lens() and Propagate() And try LensFarfield()
```python
from LightPipes import *
import matplotlib as mpl
import numpy as np
import calRediuswithPow
mpl.use('TkAgg')
import matplotlib.pyplot as plt
#import Utils
# wavelength = 500*nm
# size = 7*mm
# N = 100
# w0 = 1*mm
# f = 1*m
# z = 1*m
wavelength = 1000 * nm
N = 3000 #2840
size = 125 * mm
f = 5.5*m
z = f
w0 = 1*mm
# dx = wavelength*z/size
# N = int(size/dx)
F = Begin(size, wavelength, N)
F= CircAperture(R=50*mm,Fin=F)
F = Lens(F, f)
# F = Fresnel(F, z)
# F = Forvard(F, z)
F = Propagate(F, f)
# F2=LensFresnel(f2,f,F2);
# F=LensFarField(F,f)
# F=Forward(F)
#############################Beita calculate start###############################################
d = 100*mm
FToCalB=F
Dpow=calRediuswithPow.calRediuswithPow(FToCalB,0.839)*2
DTherory=2.44*FToCalB.lam*f/d
beita=Dpow/DTherory
print('The Beita is :',beita)
#############################Beita calculate end###############################################
```
TrySphericalCoordinates.py
### Try use Spherical Coordinatesto calculate
```python
from LightPipes import *
import matplotlib.pyplot as plt
#import Utils
import calRediuswithPow
def TheExample(N):
fig = plt.figure(figsize=(15, 9.5))
ax1 = fig.add_subplot(221)
ax2 = fig.add_subplot(222)
ax3 = fig.add_subplot(223)
ax4 = fig.add_subplot(224)
labda = 1000 * nm
size = 10 * mm
f = 100 * cm
f1 = 10 * m
f2 = f1 * f / (f1 - f)
frac = f / f1
newsize = frac * size
w = 10 * mm
w=w/5
#something wrong#########################################
# wavelength = 1000 * nm
# # N = 500
# size = 125 * mm
# f = 500 * mm
# z = 500 * mm
# w0 = 1 * mm
############################################
F = Begin(size, labda, N)
F = CircAperture(F, R=w)
# 1) Using Lens and Fresnel:
F1 = Lens(f, 0, 0, F)
# F1 = Fresnel(f, F1)
F1 = Propagate(F1,f)
phi1 = Phase(F1)
phi1 = PhaseUnwrap(phi1)
I1 = Intensity(0, F1)
x1 = []
for i in range(N):
x1.append((-size / 2 + i * size / N) / mm)
# 2) Using Lens + LensFresnel and Convert:
F2 = Lens(f1, 0, 0, F)
F2 = LensFresnel(f2, f, F2)
F2 = Convert(F2)
phi2 = Phase(F2)
phi2 = PhaseUnwrap(phi2)
I2 = Intensity(0, F2)
#Utils.show(F2,"F2 Inten")
x2 = []
for i in range(N):
x2.append((-newsize / 2 + i * newsize / N) / mm)
ax1.plot(x1, phi1[int(N / 2)], 'k--', label='Without spherical coordinates')
ax1.plot(x2, phi2[int(N / 2)], 'k', label='With spherical coordinates')
ax1.set_xlim(-newsize / 2 / mm, newsize / 2 / mm)
ax1.set_ylim(-2, 4)
ax1.set_xlabel('x [mm]')
ax1.set_ylabel('phase [rad]')
ax1.set_title('phase, N = %d' % N)
legend = ax1.legend(loc='upper center', shadow=True)
ax2.plot(x1, I1[int(N / 2)], 'k--', label='Without spherical coordinates')
ax2.plot(x2, I2[int(N / 2)], 'k', label='With spherical coordinates')
ax2.set_xlim(-newsize / 2 / mm, newsize / 2 / mm)
ax2.set_ylim(0, 1000)
ax2.set_xlabel('x [mm]')
ax2.set_ylabel('Intensity [a.u.]')
ax2.set_title('intensity, N = %d' % N)
legend = ax2.legend(loc='upper center', shadow=True)
ax3.imshow(I1)
ax3.axis('off')
ax3.set_title('Without spherical coordinates')
ax3.set_xlim(int(N / 2) - N * frac / 2, int(N / 2) + N * frac / 2)
ax3.set_ylim(int(N / 2) - N * frac / 2, int(N / 2) + N * frac / 2)
ax4.imshow(I2)
ax4.axis('off')
ax4.set_title('With spherical coordinates')
plt.figtext(0.3, 0.95, 'Spherical Coordinates, f = 100cm lens\nGrid dimension is: %d x %d pixels' % (N, N),
fontsize=18, color='red')
#############################Beita calculate start###############################################
d = w
FToCalB=F2
Dpow=calRediuswithPow.calRediuswithPow(FToCalB,0.839)*2
DTherory=2.44*FToCalB.lam*f/d
beita=Dpow/DTherory
print('The Beita is :',beita)
#############################Beita calculate end###############################################
# TheExample(100) # 100 x 100 grid
TheExample(1000) # 1000 x 1000 grid
plt.show()
```
Hello, @FredvanGoor
I am simulating the calculation of a circular aperture beam passing through a lens and converging at the focal plane. I want to calculate what is the radius of the spot when the energy in the barrel at the focal plane is 83.9%. And calculate its times of the diffraction limit.
If the diameter of the center of the first dark ring (Airy disk) of a Fraunhofer diffraction pattern is taken as the diffraction limit diameter.Theoretically, When focused through an ideal lens to the focal plane, the measured ideal diffraction limit diameter should be $d{difflimit} =\theta {difflimit} f$ = $2.44 λ * f / D$, where λ is the wavelength, f is the focal length, and D is the diameter of the circular aperture.[1]
$$u{\mathrm{~ref}}=\frac{\int{0}^{2\pi}\int{0}^{d{u,\mathrm{ref}}/2 }I\left(r\right)\mathrm{d}r\mathrm{d}\theta}{\iint I\left(r\right)\mathrm{d}r\mathrm{d}\theta}$$
Formula 1
According to Formula 1, it can be calculated that $u{ref}$ =0.839,when $d{u,\mathrm{ref}}$ is taken as diffraction limit diameter metioned above
There are some details for $u_{ref}$
According to Fraunhofer diffraction theory, the far-field intensity distribution of the reference beam with radius $r$ is calculated in paraximal approximation see A.1 $$I\left(r\right)=I_0\left\lbrace\frac{2J_{1}\left[\pi d_{u,near}r/\left(f\lambda\right)\right]}{\pi d_{u,near}r/\left(f\lambda\right)}-\varepsilon^2\frac{2J_{1}\left[\pi\varepsilon d_{u,near}r/\left(f\lambda\right)\right]}{\pi\varepsilon d_{u,near}r/\left(f\lambda\right)}\right\rbrace^2\cdots\cdots\cdots\cdots(\text{ A.}1)$$ $I_{0}$ -- the reference beam far-field peak intensity; $J_{1}$ -- first-order Bessel function; $d_\mathrm{u.near}$ -- reference beam near-field beam diameter; $f$ -- device equivalent focal length; $\lambda$ -- actual beam laser wavelength; $\varepsilon$ - center blocking ratio. (center circle area as a proportion of total area) Let $I(r)=0$, the positional equation to be satisfied by the zero radius of the circular beam intensity is given in Eq. (A.2). $J_1[\pi d_\text{u.near}r/(f\lambda)]-\varepsilon J_1[\pi\varepsilon d_\text{u,near}r/(f\lambda)]=0\quad........................( A.2 )$ According to Formula 1, when $\varepsilon$ = 0, it can be calculated that c=2.44, $u_{ref}$ =0.839So the calculation result we would expect should be $\beta = d{simulate}/d{therory} =1 $,where $d{simulate}$ is the circular diameter of the region that accounts for 83.9% of the total energy in the simulation results，and $d{therory}$ = $2.44 F.lam f / w$
But in practice, the code didn't work as I thought it would.
Code Parts
At first, I tried combining the Lens() and Propagate() functions. Then ,I tried LensFarfield() Fuction And now,I am trying use Spherical coordinates
calRediuswithPow.py is the base function I used to calculate it, there are some things wrong with its algorithm, such as the center of mass calculating itself, and the energy calculating itself, but I don't think it's going to be particularly wrong for the same field and it's propagation, and during debugging, the energy percentage was not more than 0.3% wrong
They all produce beta results that are way off compared to one, generally speaking aspherical coordinate calculations will be about 30% off, and will change greatly with focal length, and diameter of the circular hole,even if I use "ideal sampling": $dx=\frac{\lambda *z}{L}$ [2]73-74. Spherical coordinates in my example are calculated to be 0.3-0.4, which is off by 60-70%.
I'm not sure what went wrong, the radius calculations should not logically produce such a large deviation. Is there something I am not taking into account in my simulation method.
calRediuswithPow.py
### use calRediuswithPow() function to calculate the radius corresponding to energy in specify bucket ```python from LightPipes import * import matplotlib as mpl import numpy as np import copy as _copy def calRediuswithPow(fieldIn,energyPercent=0.839): ''' Calculate the radius corresponding to the ring energy Parameters ---------- energyPercent : double fieldIn : lp.Field The input field. Returns ------- The Redius correspend to the energyPercent power. """ ---------- ''' fieldIn=_copy.copy(fieldIn) I=Intensity(fieldIn) sumEnergy=0 # sumEnergy=sum(I) sumEnergy1=0 sumEnergy1=0 F_row=np.size(fieldIn.field,0) F_col=np.size(fieldIn.field,1) for i in range(F_row): for j in range(F_col): sumEnergy+=I[i][j] #sumEnergy2=Power(fieldIn) targetEnergy=energyPercent*sumEnergy # Y, X = fieldIn.mgrid_cartesian #Determine spot center of mass Coordinate system with array indexed x,y coordinates x_power_sum=0 y_power_sum=0 for i in range(F_row): for j in range(F_col): x_power_sum+=I[i][j]*i for i in range(F_row): for j in range(F_col): y_power_sum+=I[i][j]*j x_center=x_power_sum/sumEnergy y_center=y_power_sum/sumEnergy #Calculate the approximation of the true radius using the Bisection method. # First set the maximum value of the radius to the diagonal of the matrix r_max=np.sqrt(F_row**2+F_col**2) r_min=0 girdR=binarySearchR(I,r_max,r_min,targetEnergy,x_center,y_center,F_row,F_col) Rres=girdR*fieldIn.dx return Rres def binarySearchR(intensityIn,r_max,r_min,targetEnergy,x_center,y_center,F_row,F_col): ''' Bisection method to approximate the real radius (the fitting effect is a bit poor, later can be changed to a point by point cascade outward expansion) The error is close to 0.3% ''' while r_min <= r_max: r_mid=(r_min+(r_max-r_min)/2) midEnergy=calPowerWithRadius(intensityIn,F_row,F_col,x_center,y_center,r_mid) if midEnergy == targetEnergy: return r_mid elif midEnergy>targetEnergy: r_max=r_mid-1 else: r_min=r_mid+1 #close return r_mid #TODO:Sort, compute distance from center of mass for all cells, put into set sort, set map def calPowerWithRadius(intensityIn,F_row,F_col,x_center,y_center,radius): ''' Calculate the power within a specified radius ''' powerRes=0 for i in range(F_row): for j in range(F_col): if ((i-x_center)**2+(j-y_center)**2)<=radius**2: powerRes+=intensityIn[i][j] return powerRes ```LensAndLensFarfieldTry.py
### try combine Lens() and Propagate() And try LensFarfield() ```python from LightPipes import * import matplotlib as mpl import numpy as np import calRediuswithPow mpl.use('TkAgg') import matplotlib.pyplot as plt #import Utils # wavelength = 500*nm # size = 7*mm # N = 100 # w0 = 1*mm # f = 1*m # z = 1*m wavelength = 1000 * nm N = 3000 #2840 size = 125 * mm f = 5.5*m z = f w0 = 1*mm # dx = wavelength*z/size # N = int(size/dx) F = Begin(size, wavelength, N) F= CircAperture(R=50*mm,Fin=F) F = Lens(F, f) # F = Fresnel(F, z) # F = Forvard(F, z) F = Propagate(F, f) # F2=LensFresnel(f2,f,F2); # F=LensFarField(F,f) # F=Forward(F) #############################Beita calculate start############################################### d = 100*mm FToCalB=F Dpow=calRediuswithPow.calRediuswithPow(FToCalB,0.839)*2 DTherory=2.44*FToCalB.lam*f/d beita=Dpow/DTherory print('The Beita is :',beita) #############################Beita calculate end############################################### ```TrySphericalCoordinates.py
### Try use Spherical Coordinatesto calculate ```python from LightPipes import * import matplotlib.pyplot as plt #import Utils import calRediuswithPow def TheExample(N): fig = plt.figure(figsize=(15, 9.5)) ax1 = fig.add_subplot(221) ax2 = fig.add_subplot(222) ax3 = fig.add_subplot(223) ax4 = fig.add_subplot(224) labda = 1000 * nm size = 10 * mm f = 100 * cm f1 = 10 * m f2 = f1 * f / (f1 - f) frac = f / f1 newsize = frac * size w = 10 * mm w=w/5 #something wrong######################################### # wavelength = 1000 * nm # # N = 500 # size = 125 * mm # f = 500 * mm # z = 500 * mm # w0 = 1 * mm ############################################ F = Begin(size, labda, N) F = CircAperture(F, R=w) # 1) Using Lens and Fresnel: F1 = Lens(f, 0, 0, F) # F1 = Fresnel(f, F1) F1 = Propagate(F1,f) phi1 = Phase(F1) phi1 = PhaseUnwrap(phi1) I1 = Intensity(0, F1) x1 = [] for i in range(N): x1.append((-size / 2 + i * size / N) / mm) # 2) Using Lens + LensFresnel and Convert: F2 = Lens(f1, 0, 0, F) F2 = LensFresnel(f2, f, F2) F2 = Convert(F2) phi2 = Phase(F2) phi2 = PhaseUnwrap(phi2) I2 = Intensity(0, F2) #Utils.show(F2,"F2 Inten") x2 = [] for i in range(N): x2.append((-newsize / 2 + i * newsize / N) / mm) ax1.plot(x1, phi1[int(N / 2)], 'k--', label='Without spherical coordinates') ax1.plot(x2, phi2[int(N / 2)], 'k', label='With spherical coordinates') ax1.set_xlim(-newsize / 2 / mm, newsize / 2 / mm) ax1.set_ylim(-2, 4) ax1.set_xlabel('x [mm]') ax1.set_ylabel('phase [rad]') ax1.set_title('phase, N = %d' % N) legend = ax1.legend(loc='upper center', shadow=True) ax2.plot(x1, I1[int(N / 2)], 'k--', label='Without spherical coordinates') ax2.plot(x2, I2[int(N / 2)], 'k', label='With spherical coordinates') ax2.set_xlim(-newsize / 2 / mm, newsize / 2 / mm) ax2.set_ylim(0, 1000) ax2.set_xlabel('x [mm]') ax2.set_ylabel('Intensity [a.u.]') ax2.set_title('intensity, N = %d' % N) legend = ax2.legend(loc='upper center', shadow=True) ax3.imshow(I1) ax3.axis('off') ax3.set_title('Without spherical coordinates') ax3.set_xlim(int(N / 2) - N * frac / 2, int(N / 2) + N * frac / 2) ax3.set_ylim(int(N / 2) - N * frac / 2, int(N / 2) + N * frac / 2) ax4.imshow(I2) ax4.axis('off') ax4.set_title('With spherical coordinates') plt.figtext(0.3, 0.95, 'Spherical Coordinates, f = 100cm lens\nGrid dimension is: %d x %d pixels' % (N, N), fontsize=18, color='red') #############################Beita calculate start############################################### d = w FToCalB=F2 Dpow=calRediuswithPow.calRediuswithPow(FToCalB,0.839)*2 DTherory=2.44*FToCalB.lam*f/d beita=Dpow/DTherory print('The Beita is :',beita) #############################Beita calculate end############################################### # TheExample(100) # 100 x 100 grid TheExample(1000) # 1000 x 1000 grid plt.show() ```Others
Spherical coordinates is a little hard to understand, for example I don't know how to set f1 and z. @ldoyle 's xplanation confused me even more
End
Thanks for your reading! I've been puzzling over these calculations for a long time, so if you have time please help me!
References
[1] Hecht E. Optics[M]. 5.Boston: Pearson Education, Inc.,2017: vi, 480-484. [2] Voelz D G.Computational fourier optics: a MATLAB tutorial[J],2011.