Calculate the power( ratio) in the bucket and times the diffraction limit

[tianjiaoshanzai]

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.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

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()

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

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.