A numerical simulation method for optical field transmission based on matrix operations

Abstract

This invention discloses a numerical simulation method for optical field transmission based on matrix operations, belonging to the field of optical field transmission research. First, the matrix of the optical transmission system is determined according to the task requirements. Then, the generalized Huygens-Fresnel diffraction integral is transformed into a summation form in discrete rectangular coordinates, and further transformed into a matrix product form. Using the determined optical transmission system matrix parameters and the determined numerical integration algorithm matrix, the target optical field is finally calculated using the input optical field and the algorithm matrix. This invention avoids the Gibbs phenomenon, confusion, and tailing effects caused by Fourier transform. This invention allows for arbitrary selection of the number of grid points in the target optical field, offering flexibility and economy.

Description

Technical Field

[0001] This invention belongs to the field of optical field transmission technology, specifically relating to a numerical simulation method for optical field transmission based on matrix operations. Background Technology

[0002] With the continuous development of laser technology, lasers are widely used in many fields such as industry, communication, medicine and national defense. Laser transmission is one of the important links in the application of laser systems, and mainly studies the various physical processes that occur when lasers are emitted from the system outlet and transmitted in the atmosphere or other media. Laser transmission problems can be studied through physical experiments and numerical simulations to obtain regular understanding. (Ren Bin, Chen Chunyi, Yang Huamin. A review of numerical simulation technology for light wave transmission in atmospheric turbulence [J]. Journal of System Simulation, 2017, 29(08): 1631-1640.) Numerical simulation is favored by researchers for its advantages of ease of implementation, controllable parameters, and cost-effectiveness. (Sun Quan, Lü Pin, Ning Yu, Xi Fengjie, Liu Wenguang, Xu Xiaojun. Application of Seelight optical system simulation software in adaptive optics [J]. Optoelectronic Engineering, 2018, 45(03): 130-136.) Among the numerical simulation methods for laser transmission, the angular spectrum propagation method is widely used because of its clear physical picture. However, the angular spectrum propagation method inevitably has the following limitations due to the involvement of Fourier transform: First, the limited discrete sampling frequency will cause signal oscillations at abrupt changes or boundaries in the optical field; second, the limited discrete Fourier transform will cause confusion and tailing effects; third, the angular spectrum propagation method has complex requirements for the computational grid parameters, which makes the calculation inconvenient and inflexible, and may even lead to incorrect calculation results due to improper grid selection; fourth, the angular spectrum propagation method requires that the number of grids in the input optical field and the target optical field must be the same, which is computationally expensive for simulation scenarios that only care about the Strelby ratio or the power in the barrel (Schmidt, Jason D. Numerical simulation of optical wave propagation with examples in MATLA[B], (SPIE Press Monograph Vol.PM199,2010)). Summary of the Invention

[0003] The purpose of this invention is to provide a numerical simulation method for optical field transmission based on matrix operations, which solves the technical problems of traditional angular spectrum propagation method in optical field transmission numerical simulation. The finite discrete sampling frequency leads to signal oscillation at abrupt changes or boundaries in the optical field, and the method has relatively complex requirements on the network parameters involved in the calculation, resulting in inflexible calculation or even erroneous results and high computational cost.

[0004] To achieve the above objectives and solve the above technical problems, the technical solution of the present invention is as follows:

[0005] A numerical simulation method for optical transmission based on matrix operations includes the following steps:

[0006] S1: Determine the ABCD matrix of the optical transmission system based on mission requirements.

[0007] S2: After converting the generalized Huygens-Fresnel diffraction integral of the target light field into a summation form in a discrete rectangular coordinate system, it prepares for conversion into a matrix product form;

[0008] The generalized Huygens-Fresnel integral of the target light field is shown in equation (1).

[0009]

[0010] Among them, (x m ,y n Let m, n ∈ (1, N1), and (x) be the Cartesian coordinates of the input light field. j ,y k ) represents the Cartesian coordinates of the target light field; j,k∈(1,N2), U is the target light field, u is the input light field; A, B, D are the parameters in the ABCD matrix of the optical transmission system determined in step 1;

[0011] The generalized Huygens-Fresnel diffraction integral is transformed into a summation form in discrete rectangular coordinates by numerical integration, as shown in equation (2).

[0012]

[0013] Where W is a vector related to numerical integration, and the vector W contains N1 elements, and Δ is the grid spacing in the rectangular coordinate system of the input light field;

[0014] The numerical integration method chosen is the gradient algorithm, where W = (0.5, 1, 1, ..., 0.5) × Δ

[0015] Alternatively, Simpson's algorithm or other suitable numerical integration methods can be chosen. For Simpson's algorithm:

[0016] have

[0017] By rearranging terms in formula (2), we obtain formula (3) as shown below.

[0018]

[0019] because and Because they are highly similar, it is specifically stipulated that:

[0020] For (x) m ,y n Given x1 = y1, x2 = y2, x3 = y3, x4 = y4, ..., x N1 =y N1 That is, the Cartesian coordinate system in which the input light field u is located is a matrix with N1×N1 grids;

[0021] For (x) j ,y k Given x1 = y1, x2 = y2, x3 = y3, x4 = y4, ..., x N2 =y N2 That is, the Cartesian coordinate system in which the output light field U is located is a matrix with a grid size of N2×N2;

[0022] By using this rule, the abscissa interval and ordinate interval of the rectangular coordinate system where the input light field is located are equal, so the numerical integration vector W is completely equal to the abscissa and ordinate of the rectangular coordinate system where the input light field is located.

[0023] Therefore, formula (3) can be transformed into the form of matrix multiplication as shown in formula (4):

[0024]

[0025] Where H is the introduced algorithm matrix, with a size of N2×N1, and H is defined as follows:

[0026]

[0027] Where x m ,y n These are the horizontal and vertical coordinate vectors of the input light field, x and y, respectively. j ,y k These are the horizontal and vertical coordinate vectors of the target light field, respectively.

[0028] This transforms the generalized Huygens-Fresnel diffraction integral into a summation form in discrete rectangular coordinates.

[0029] S3: Substitute the ABCD matrix parameters determined in S1 and W, which is related to the numerical integration method, into equation (4) to solve for H;

[0030] S4: Substitute the input light field u and the matrix H obtained in S3 into Equation 3 to calculate the target light field U in the form of the following matrix product, as shown in Equation (6):

[0031]

[0032] Wherein, the input light field u is an N1×N1 matrix, the target light field U is an N2×N2 matrix, and H is an N1×N2 matrix.

[0033] The beneficial effects of this invention are as follows:

[0034] This invention transforms the generalized Huygens-Fresnel integral into a simple matrix multiplication operation by selecting a suitable numerical integration algorithm and utilizing the commutativity of the ordinate and abscissa in the generalized Huygens-Fresnel integral. This method avoids the Gibbs phenomenon, confusion, and tailing effects associated with Fourier transforms. Furthermore, this method allows for arbitrary selection of the number of grid cells in the target light field, offering flexibility and cost-effectiveness.

[0035] Compared to traditional optical field transmission simulation methods—the Fourier transform method (or angular spectrum propagation method)—this invention offers advantages in terms of speed, accuracy, simplicity, and flexibility. Specifically, it avoids the Gibbs phenomenon—oscillations at abrupt changes or boundaries in the optical field caused by the limited sampling frequency in the Fourier transform; secondly, it avoids the confusion and tailing effects caused by improper frequency sampling in the discrete Fourier transform; thirdly, the flexible operation of the numerical integration method simplifies the selection of mesh parameters in the optical transmission numerical simulation process; and fourthly, the calculation of the target optical field can be limited to a specific point or region without calculating the entire target optical field. For numerical simulation scenarios where the focus is on specific physical parameters of the optical field—such as the Strelby ratio or the power in the barrel—the simulation speed can be further improved.

[0036] This invention can replace the traditional angular spectrum propagation algorithm for fast and accurate simulation of light field propagation in vacuum, atmosphere or other media. Attached Figure Description

[0037] Figure 1 This is a schematic diagram of the technical solution of the present invention;

[0038] Figure 2 This is a grayscale image showing the near-field intensity distribution of a rectangular light beam.

[0039] Figure 3 A comparison of numerical simulation results of the present invention and the Fourier transform method in a rectangular beam vacuum transmission scenario;

[0040] Figure 4 Grayscale image of the fundamental Gaussian beam;

[0041] Figure 5 This paper compares the numerical simulation results of the present invention and the Fourier transform method in the vacuum transmission scenario of the fundamental mode Gaussian beam. Detailed Implementation

[0042] To better understand this invention, the design principles of this invention are as follows: Figure 1 As shown, Figure 1The right side shows the near-field spot pattern of the Gaussian beam, and the left side shows the far-field spot pattern of the Gaussian beam. The propagation distance from the near field to the far field is 300 meters. The formula in the middle is the detailed process of converting the generalized Huygens-Fresnel diffraction integral into matrix multiplication.

[0043] The present invention will now be further described with reference to the accompanying drawings and specific examples.

[0044] A numerical simulation method for optical transmission based on matrix operations includes the following steps:

[0045] S1: Determine the ABCD matrix of the optical transmission system based on mission requirements.

[0046] In geometric optics, light rays, optical systems, and transmission media have specific matrix representation methods. The transformation matrix of a transmission system can be represented as the product of the transformation matrices of each subsystem. This method belongs to the category of matrix optics.

[0047] S2: After converting the generalized Huygens-Fresnel diffraction integral of the target light field into a summation form in a discrete rectangular coordinate system, it is further converted into a matrix product form;

[0048] The core concept of this invention is to avoid using the Discrete Fourier Transform (DFT), thereby avoiding the Gibbs phenomenon of oscillations at abrupt changes or boundaries in the optical field caused by the limited sampling frequency in the DFT. It also avoids the confusion and tailing effects caused by improper frequency sampling in the DFT. To solve the aforementioned technical problems, the generalized Huygens-Fresnel diffraction integral was chosen. However, directly applying the generalized Huygens-Fresnel diffraction integral is neither appropriate nor practical due to the large computational cost. Therefore, this invention sets the x and y coordinates of the input and target optical fields to be consistent, and then utilizes the interchangeable nature of their horizontal and vertical coordinates to transform the generalized Huygens-Fresnel integral into matrix multiplication using a numerical integration algorithm.

[0049] Therefore, this invention transforms the generalized Huygens-Fresnel integral of the target light field into matrix multiplication by selecting a suitable numerical integration algorithm, including but not limited to gradient algorithm and Simpson algorithm;

[0050] The generalized Huygens-Fresnel integral of the target light field is shown in equation (1).

[0051]

[0052] Among them, (x m ,y n Let m, n ∈ (1, N1), and (x) be the Cartesian coordinates of the input light field. j ,y k) represents the Cartesian coordinates of the target light field; j,k∈(1,N2), U is the target light field, u is the input light field; a, b, c are the parameters in the ABCD matrix of the optical transmission system determined in step 1;

[0053] The generalized Huygens-Fresnel diffraction integral is transformed into a summation form in discrete rectangular coordinates by numerical integration, as shown in equation (2).

[0054]

[0055] Where W is a vector related to numerical integration, and the vector W contains N1 elements, and Δ is the grid spacing in the rectangular coordinate system of the input light field;

[0056] Applying a phase-shifting operation to formula (2) yields formula (3) as shown below.

[0057]

[0058] It is noted here that and Because they are highly similar, it is specifically stipulated that:

[0059] For (x) m ,y n Given x1 = y1, x2 = y2, x3 = y3, x4 = y4, ..., x N1 =y N1 That is, the Cartesian coordinate system in which the input light field u is located is a matrix with N1×N1 grids;

[0060] For (x) j ,y k Given x1 = y1, x2 = y2, x3 = y3, x4 = y4, ..., x N2 =y N2 That is, the Cartesian coordinate system in which the output light field U is located is a matrix with a grid size of N2×N2;

[0061] By using this rule, we can obtain that the interval between the horizontal and vertical coordinates of the rectangular coordinate systems where the input light field and the target light field are located is equal. Therefore, the horizontal and vertical coordinates of the numerical integration vector W are completely equal with respect to the rectangular coordinate system where the target light field is located.

[0062] Therefore, formula (3) can be transformed into the form of matrix multiplication as shown in formula (4):

[0063]

[0064] Where H is the introduced algorithm matrix, with a size of N2×N1, and H is defined as follows:

[0065]

[0066] Where x m ,y n These are the horizontal and vertical coordinate vectors of the input light field, x and y, respectively. j ,y k These are the horizontal and vertical coordinate vectors of the target light field, respectively.

[0067] This transforms the generalized Huygens-Fresnel diffraction integral of the target square into a summation form in a discrete rectangular coordinate system.

[0068] S3: Substitute the ABCD matrix parameters (a, b, d) determined in S1 and the W related to the numerical integration method determined in S2 into equation (4) to solve for H;

[0069] For example, the gradient algorithm can be used to demonstrate the process of calculating H, where W = (0.5, 1, 1, ..., 0.5) × Δ.

[0070]

[0071] Δ represents the grid spacing, (x1, x2, ..., x N1 Let (y1, y2, ..., y3) be the abscissa of the input light field. N1 (x'1, x'2, ..., x') represents the x and y coordinates of the input light field. N2 Let (y'1, y'2, ..., y') be the abscissa of the target light field. N2 () represents the ordinate of the target light field;

[0072] S4: Substitute the input light field u and the matrix H obtained in S3 into equation (4) to calculate the target light field U in the form of the following matrix product, as shown in equation (6):

[0073]

[0074] Wherein, the input light field u is an N1×N1 matrix, the target light field U is an N2×N2 matrix, and H is an N1×N2 matrix.

[0075] In step S2, the numerical integration method selected is the gradient algorithm, or Simpson's algorithm or other suitable numerical integration methods can be chosen. For Simpson's algorithm:

[0076] have

[0077] Example 1: Vacuum Transmission of Matrix Beams

[0078] A numerical simulation method for optical field transmission based on matrix operations is used to simulate the vacuum transmission of a rectangular beam, including the following steps:

[0079] S1: Determine the ABCD matrix of the transmission system.

[0080] The ABCD matrix of vacuum transmission Fresnel diffraction is

[0081]

[0082] S2: Convert the generalized Huygens-Fresnel integral into matrix multiplication using a gradient algorithm.

[0083]

[0084] S3: Substitute the parameters (a, b, d) of the ABCD matrix in S1 and the W related to the numerical integration method in S2 into Equation 2 to solve for H. N1×N2

[0085]

[0086] in Input the x-coordinates of the light field grid points. denoted as the x-coordinate of the target light field grid point.

[0087] In this example, N1 and N2 are both selected as 512, Δ is selected as 0.1mm, and the transmission distance z is selected as 300m.

[0088] S4: Substitute the input light field u (rectangular beam) and the matrix H obtained in step 3 into Equation 3 to calculate the target light field. In this example, the rectangular beam is as follows: Figure 2 As shown, the light wavelength is 1μm, the light power is 1W, and the light spot size is 8mm.

[0089] Figure 3 This paper compares the numerical simulation results of the matrix method and the Fourier transform method of this invention in a rectangular beam vacuum transmission scenario. It can be seen that the matrix method agrees better with the analytical solution than the Fourier transform method. In this example, matrix operations are implemented using the BLAS numerical computation toolkit. For a transmission distance of 300m, the Fourier transform method takes 2.239s; while the matrix algorithm of this invention takes only 0.343s, approximately six times faster. This example demonstrates the advantages of the matrix operation-based numerical simulation method for light field transmission of this invention, which is fast and accurate.

[0090] Example 2: Vacuum transmission of fundamental mode Gaussian beam

[0091] A numerical simulation method for optical field transmission based on matrix operations is used to simulate the vacuum transmission of a fundamental mode Gaussian beam, including the following steps:

[0092] S1: Determine the ABCD matrix of the transmission system.

[0093] The ABCD matrix of vacuum transmission Fresnel diffraction is

[0094]

[0095] S2: Convert the generalized Huygens-Fresnel integral into matrix multiplication using a gradient algorithm.

[0096]

[0097]

[0098] W = (0.5, 1, 1, ..., 0.5) × Δ

[0099] S3: Substitute the ABCD matrix parameters (a, b, d) from S1 and the W related to the numerical integration method from S2 into Equation 2 to solve for H. N1×N2

[0100]

[0101] in Input the x-coordinates of the light field grid points. denoted as the x-coordinate of the target light field grid point.

[0102] In this example, N1 and N2 are both selected as 512, Δ is selected as 0.1mm, and the transmission distance z is selected as 360m.

[0103] S4: Substitute the input light field u (fundamental mode Gaussian beam) and the matrix H from step 3 into Equation 3 to calculate the target light field. In this example, the fundamental mode Gaussian beam is as follows: Figure 3 As shown, the light wavelength is 1μm, the light power is 1W, and the light spot size is 8mm.

[0104] Figure 5 This paper compares the numerical simulation results of the matrix method and the Fourier transform method of this invention in a vacuum transmission scenario of a fundamental Gaussian beam. It can be seen that the matrix method is highly consistent with the analytical solution, while the Fourier transform method exhibits the Gibbs phenomenon, resulting in signal oscillation. In this example, matrix operations are implemented using the BLAS numerical computation toolkit. For a transmission distance of 360m, the Fourier transform method takes 2.199s; while the matrix method of this invention takes 0.632s, approximately three times faster. This example also demonstrates the advantages of the matrix operation-based numerical simulation method for optical field transmission of this invention, which is fast and accurate.

Claims

1. A numerical simulation method for optical transmission based on matrix operations, characterized in that, Includes the following steps: S1: Determine the ABCD matrix of the optical transmission system based on mission requirements. S2: After converting the generalized Huygens-Fresnel diffraction integral of the target light field into a summation form in a discrete rectangular coordinate system, it prepares for conversion into a matrix product form; The generalized Huygens-Fresnel integral of the target light field is shown in equation (1). Among them, (x m ,y n Let m, n be the Cartesian coordinates of the input light field; m, n ∈ (1, N1). (x j ,y k ) represents the Cartesian coordinates of the target light field; j,k∈(1,N2), U is the target light field, u is the input light field; A, B, D are the parameters in the ABCD matrix of the optical transmission system determined in step 1; The generalized Huygens-Fresnel diffraction integral is transformed into a summation form in discrete rectangular coordinates by numerical integration, as shown in equation (2). Where W is a vector related to numerical integration, and the vector W contains N1 elements, and Δ is the grid spacing in the rectangular coordinate system of the input light field; By rearranging terms in formula (2), we obtain formula (3) as shown below. because and Because they are highly similar, it is specifically stipulated that: For (x) m ,y n Given x1 = y1, x2 = y2, x3 = y3, x4 = y4, ..., x N1 =y N1 That is, the Cartesian coordinate system in which the input light field u is located is a matrix with N1×N1 grids; For (x) j ,y k Given x1 = y1, x2 = y2, x3 = y3, x4 = y4, ..., x N2 =u N2 That is, the Cartesian coordinate system in which the output light field U is located is a matrix with a grid size of N2×N2; By using this rule, the interval between the horizontal and vertical coordinates of the rectangular coordinate system in which the input light field is located is equal, so the numerical integration vector W is completely equal to the horizontal and vertical coordinates of the rectangular coordinate system in which the input light field is located. Therefore, formula (3) can be transformed into the form of matrix multiplication as shown in formula (4): Where H is the introduced algorithm matrix, with a size of N2×N1, and H is defined as follows: Where x m ,y n These are the horizontal and vertical coordinate vectors of the input light field, x and y, respectively. j ,y k These are the abscissa and ordinate vectors of the target light field, respectively; This transforms the generalized Huygens-Fresnel diffraction integral into a summation form in discrete rectangular coordinates. S3: Substitute the ABCD matrix parameters determined in S1 and W, which is related to the numerical integration method, into equation (4) to solve for H; S4: Substitute the input light field u and the matrix H obtained in S3 into Equation 3 to calculate the target light field U in the form of the following matrix product, as shown in Equation (6): Wherein, the input light field u is an N1×N1 matrix, the target light field U is an N2×N2 matrix, and H is an N1×N2 matrix.

2. The numerical simulation method for optical transmission based on matrix operations according to claim 1, characterized in that, The calculation of the target light field u can be performed by selecting only a point or a region without calculating the entire target light field, which saves computational costs and effectively improves computational efficiency.

3. A numerical simulation method for optical transmission based on matrix operations according to claim 1 or 2, characterized in that, In step S2, the numerical integration methods selected are gradient algorithm and Simpson algorithm.

4. A numerical simulation method for optical transmission based on matrix operations according to claim 1 or 2, characterized in that, For Simpson's algorithm, there is 5. A numerical simulation method for optical transmission based on matrix operations according to claim 1 or 2, characterized in that, For the gradient algorithm, we have W = (0.5, 1, 1, ... 0.5) × Δ.

Images