METHOD FOR PROCESSING ANISOTROPIC MAGNETIZED PLASMA MEDIUM AND SYSTEM THEREOF

Abstract

A method for processing anisotropic magnetized plasma medium and a system thereof are provided. The method includes obtaining the Maxwell equation and the polarization current density equation based on the electromagnetic characteristics of anisotropic magnetized plasma; processing the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed; based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, numerical iterative equations for electric field intensity, magnetic field intensity, and polarization current density in anisotropic magnetized plasma medium are obtained; a numerical modeling simulation electromagnetic model is used to determine the electromagnetic characteristics of the electromagnetic model. The use of the present disclosure to simulate the propagation of electromagnetic waves in anisotropic magnetized plasma medium has higher numerical calculation accuracy.

Description

TECHNICAL FIELD

The present disclosure relates to the field of plasma medium processing, in particular to a method for processing anisotropic magnetized plasma medium and a system thereof.

BACKGROUND

Near space hypersonic vehicles have the advantages of both fast flight speed and moderate flight altitude, which brings great challenges to radar warning. At the same time, near space hypersonic vehicles can maneuver to avoid missile defense systems, which is very beneficial for military penetration and has significant military applications and research value. However, with the continuous faster of aircraft flight speed, the plasma sheath generated by hypersonic vehicle during near space flight can have an impact on the reception of measurement and control signals, and even cause communication interruption, resulting in a blackout and incalculable consequences. Therefore, the study of electromagnetic wave propagation mechanism in plasma is of great significance to overcome or even solve the problem of communication when ionization blackout occurs.

In the study of electromagnetic wave propagation in plasma, analytical calculations are difficult to obtain due to the complexity of the target model, and time-domain numerical methods play an important role. The advantage of the time-domain method lies in selecting a suitable broadband excitation source, and the electromagnetic response within the required frequency band can be obtained through a single simulation. Especially, the finite-difference time-domain method is highly concerned, with numerical characteristics of simple, intuitive, and easy to implement, it is widely used in electromagnetic wave propagation problems in plasma. However, the traditional FDTD (2,2) method only has second-order accuracy in time and space, and often requires the use of finer Yee grids to reduce numerical dispersion and anisotropy errors. In particular, plasma medium are often anisotropic, and their electromagnetic characteristics differ significantly in space. Therefore, it is necessary to use a finer grid than general medium simulations to discretize the electromagnetic model and reduce numerical errors. This will undoubtedly increase the memory usage of traditional FDTD (2,2) methods while reducing computational efficiency. In order to solve this problem, some high-order FDTD (2,4) methods with second-order numerical accuracy in time and fourth-order numerical accuracy in space have been proposed and used for dispersion medium simulation. The FDTD (2,4) method not only retains the explicit solving characteristics of the traditional FDTD (2,2) method, but also improves the numerical calculation accuracy of the time-domain method. However, the FDTD (2,4) method still cannot avoid long-term error accumulation, and its numerical stability needs to meet stricter CFL conditions. Subsequently, a symplectic finite-difference time-domain method (symplectic FDTD (4,4) and SFDTD (4,4)) with fourth-order numerical accuracy in both time and space is proposed. Unlike the ordinary FDTD (2,4) method, the symplectic algorithm is based on the fundamental principles of Hamiltonian mechanics and preserves the symplectic structure of the original system. It preserves the original symplectic structure of the discretized difference equation, and the symplectic algorithm has long-term stability and energy conservation characteristics. At present, the symplectic algorithm has been used to simulate electromagnetic propagation and electromagnetic scattering problems in anisotropic magnetized plasma. Although the previous work provided a Hamiltonian equation containing the current density term, there was no Taylor series expansion of the exponential matrix algorithm. The numerical discretization form of the SFDTD (4,4) method for current density field component is similar to that of the traditional FDTD (2,2) method. In numerical calculations, each iterative solution is simply divided into multi-level iterative solution. This symplectic discrete form only holds when the split matrix U^α=0 (U is the split matrix, α≥2). Generally speaking, when dealing with more complex dispersion model functions, the higher-order terms of the derived splitting matrix are not equal to zero, that is, U^α0 (α≥2).

SUMMARY

The objective of the present disclosure is to provide a method for processing anisotropic magnetized plasma medium and a system thereof, which can improve higher numerical calculation accuracy.

To achieve the above objective, the present disclosure provides the following solution:

• • A method for processing anisotropic magnetized plasma medium, including:
  • Obtaining a Maxwell equation and a polarization current density equation based on electromagnetic characteristics of anisotropic magnetized plasma;
  • Processing the Maxwell equation and the polarization current density equation to obtain an electric field intensity, a magnetic field intensity, and a polarization current density after processed;
  • Using a matrix exponential time-domain finite difference method to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed;
  • Applying a numerical modeling simulation electromagnetic model to determine electromagnetic characteristics of the electromagnetic model according to the numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in the anisotropic magnetized plasma medium.

Optionally, the Maxwell equation is:

$$\begin{gathered} \nabla \times E=-{\mu }_0\frac{\partial H}{\partial t},\text{ \\ ∇×H=ε0⁢∂E∂t+;} \end{gathered}$$

The polarization current density equation is:

$\frac{dJ}{dt}+v={\epsilon }_0{\omega }_p^{2}E{+}_b\times ;$

Wherein, E is the electric field intensity, H is the magnetic field intensity, J is the polarization current density, ε₀, μ₀ are respectively a dielectric constant and a permeability in vacuum, ω_p is a plasma frequency, v is a plasma collision frequency, ω_b=B₀/m_e is a electron cyclotron frequency, B₀ is an external static magnetic field, e is an electron charge, and m_e is an electron mass.

Optionally, processing the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed, specifically including:

• • Performing multi-level symplectic discretization on the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density in a symplectic discretization format:

$E^{n+l/m}=E^{n+\left(l-1\right)/m}+\frac{d_l\Delta t}{{\epsilon }_0}\nabla \times H^{n+l/m}+\frac{-d_l\Delta t}{{\epsilon }_0}{J}^{n+l/m}{H}^{n+l/m}=H^{n+\left(l-1\right)/m}+\frac{-c_l\Delta t}{{\epsilon }_0}\nabla \times E^{n+l/m}{J}^{n+l/m}={\mathrm{DJ}}^{n+\left(l-1\right)/m}+{\epsilon }_0{\omega }_p^2{\mathrm{FE}}^{n+\left(l-1\right)/m},$

Wherein, d_l and c_l are coefficients of symplectic propagator,

$\Omega =\left[\begin{array}{ccc}-v& -{\omega }_{\mathrm{bz}}& {\omega }_{\mathrm{by}}\\ {\omega }_{\mathrm{bz}}& -v& -{\omega }_{bx}\\ -{\omega }_{\mathrm{by}}& {\omega }_{bx}& -v\end{array}\right],$

l is the order of a field component =J, , )^T, and m is the number of stages of non-dissipative p-order symplectic integral, D=e^CM, F=Ω⁻¹(D−I).

Optionally, using a matrix exponential time-domain finite difference method to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed specifically includes:

• • Assuming that the plasma medium is biased by a static magnetic field in the z-direction, based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, the matrix exponential time-domain finite difference method is used to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium:

$E_x^{n+l/m}\left(i+\frac{1}{2},j,k\right)=E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)+\text{⁠}\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{y1}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{1}{2},k\right)\right]\text{⁠}+{\gamma }_{z1}\left[-H_y^{n❘l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+H_y^{n❘l/m}\left(i+\frac{1}{2},j,k-\frac{1}{2}\right)\right]+{\gamma }_{y2}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{3}{2}{\rm A},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{3}{2},k\right)\right]+{\gamma }_{z2}\left[-H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{3}{2}\right)+H_y^{n+l/m}\left(i+\frac{1}{2},j,k-\frac{3}{2}\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)$ $E_y^{n❘l/m}\left(i,j+\frac{1}{2},k\right)=E_x^{n❘\left(l1\right)/m}\left(i,j+\frac{1}{2},k\right)+\text{⁠}\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{z1}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)-H_x^{n-l/m}\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)\right]+{\gamma }_{x1}\left[-H_z^{n❘l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)+H_z^{n❘l/m}\left(i-\frac{1}{2},j+\frac{1}{2},k\right)\right]+{\gamma }_{z2}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{3}{2}\right)-H_x^{n+l/m}\left(i,j+\frac{1}{2},k-\frac{3}{2}\right)\right]+{\gamma }_{x2}\left[-H_z^{n+l/m}\left(i+\frac{3}{2},j+\frac{1}{2},k\right)+H_z^{n+l/m}\left(i-\frac{3}{2},j+\frac{1}{2},k\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_y^{n-\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)$ $E_z^{n+l/m}\left(i,j,k+\frac{1}{2},k\right)=E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2},k\right)+\text{⁠}\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{x1}\left[H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)-H_y^{n+l/m}\left(i-+\frac{1}{2},j,k+\frac{1}{2}\right)\right]\text{⁠}+{\gamma }_{y1}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)+H_x^{n+l/m}\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)\right]+{\gamma }_{x2}\left[H_y^{n+l/m}\left(i+\frac{3}{2},j,k+\frac{1}{2}\right)-H_y^{n+l/m}\left(i-\frac{3}{2},j,k+\frac{1}{2}\right)\right]\right\}+{\gamma }_{y2}\left[-H_x^{n+l/m}\left(i,j+\frac{3}{2},k+\frac{1}{2}\right)+H_x^{n+l/m}\left(i,j-\frac{3}{2},k+\frac{1}{2}\right)\right]-\frac{d_l\Delta t}{{\epsilon }_0}{J}_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)$ $H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)=H_x^{n❘\left(l1\right)/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)]+\text{⁠}\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{z1}[E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+1\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)\text{⁠}+{\gamma }_{y1}\left[-E_z^{n-\left(l-1\right)/m}\left(i,j+1,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z2}\left[E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+2\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k-1\right)\right]+{\gamma }_{y2}[-E_z^{n-\left(l-1\right)/m}\left(i,j+2,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j-1,k+\frac{1}{2}\right)\right\}$ $H_y^{n-l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)=H_y^{n-\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{x1}\left[E_z^{n+\left(l-1\right)/m}\left(i+1,j+\frac{1}{2},k\right)-E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z1}\left[-E_x^{n❘\left(l1\right)/m}\left(i+\frac{1}{2},j,k+1\right)+E_x^{n❘\left(l1\right)/m}\left(i+\frac{1}{2},j,k\right)\right]+{\gamma }_{x2}\left[E_z^{n+\left(l-1\right)/m}\left(i+2,j,k+\frac{1}{2}\right)-E_z^{n+\left(l-1\right)/m}\left(i-1,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z2}\left[-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+2\right)+E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k-1\right)\right]\right\}$ $H_z^{n-l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)=H_z^{n-\left(l-1\right)/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)\text{⁠}+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{y1}\left[E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j+1,k\right)-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)\right]\text{⁠}+{\gamma }_{x1}\left[-E_y^{n❘\left(l1\right)/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n\left(l1\right)/m}\left(i,j+\frac{1}{2},k\right)\right]+{\gamma }_{y2}\left[E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j+2,k\right)-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j-1,k\right)\right]+{\gamma }_{x2}\left[-E_y^{n+\left(l-1\right)/m}\left(i+2,j+\frac{1}{2},k\right)+E_y^{n+\left(l-1\right)/m}\left(i-1,j+\frac{1}{2},k\right)\right]\right\}$ $J_x^{n+\frac{l}{m}}\left(i+\frac{1}{2},j,k\right)=\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)\text{⁠}{e}^{-{\mathrm{vc}}_l\mathrm{dt}}\right]\text{⁠}J_x^{n+\frac{l-1}{m}}\left(\text{⁠}i+\frac{1}{2},\text{⁠}j,\text{⁠}k\right)+$ $\left[\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0\text{⁠}{\omega }_p^2\text{⁠}\left\{\frac{v\left[1-\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_p^2}+\frac{{\omega }_b\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_p^2}\right\}{E}_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0\text{⁠}{\omega }_p^2\text{⁠}\left\{\frac{v\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}]}{v^2+{\omega }_p^2}+\frac{{\omega }_b[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}-1}{v^2+{\omega }_p^2}\right\}{E}_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)$

Because a current density position is the same as an electric field position, therefore, when calculating

$J_{\lambda }^{n+l/m}\left(i+\frac{1}{2},j,k\right),$

it is necessary to perform interpolation on

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)\mathrm{and}{E}_y^{n-l/m}\left(i+\frac{1}{2},j,k\right),$

because J_y and E_y are no value at a point

$\left(i+\frac{1}{2},j,k\right),$

so it needs to be averaged by four adjacent diagonal values to obtain:

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(J_y^{n-l/m}\left(i,j-\frac{1}{2},k\right)+J_y^{n-l/m}\left(i,j+\frac{1}{2},k\right)\text{⁠}\text{⁠}+J_y^{n-l/m}\left(i+1,j+\frac{1}{2},k\right)+J_y^{n-l/m}\left(i+1,j-\frac{1}{2},k\right)\right)$ $E_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(E_Y^{n-l/m}\left(i,j-\frac{1}{2},k\right)+E_y^{n-l/m}\left(i,j+\frac{1}{2},k\right)+E_y^{n-l/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n-l/m}\left(i+1,j-\frac{1}{2},k\right)\right)$ $J_y^{n+\frac{l}{m}}\left(i,j+\frac{1}{2},k\right)=\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]{J_y^{n+}}^{\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)\text{⁠}+$ $\left[\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)\text{⁠}+$ ${\epsilon }_0{\omega }_p^2\text{⁠}\left\{\frac{v\left[1-\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}+\frac{{\omega }_b\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}\right\}$ $E_y^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)-$ ${\epsilon }_0{\omega }_p^2\left\{\frac{v\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}+\frac{{\omega }_b\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}-1\right]}{v^2+{\omega }_b^2}\right\}{E}_x^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)$

Because the current density position is the same as the electric field position, therefore, when calculating

$J_y^{n+l/m}\left(i,j+\frac{1}{2},k\right),$

it is necessary to perform interpolation on

$J_x^{n-l/m}\left(i,j+\frac{1}{2},k\right)\mathrm{and}{E}_x^{n+l/m}\left(i,j+\frac{1}{2},k\right),$

because J_x and E_x are no value at a point

$\left(i,j+\frac{1}{2},k\right),$

so it needs to be averaged by four adjacent diagonal values to obtain:

$$\begin{gathered} J_x^{n-l/m}\left(i,j+\frac{1}{2},k\right)=\frac{1}{4}\left(J_x^{n-l/m}\left(i-\frac{1}{2},j,k\right)+J_x^{n-l/m}\left(i+\frac{1}{2},j,k\right)+\text{ \\ Jxn-l/m(i+12,j+1,k)+Jxn-l/m(i-12,j+1,k)}\right)E_x^{n+l/m}\left(i,j+\frac{1}{2},k\right)=\frac{1}{4}\left(E_x^{n-l/m}\left(i-\frac{1}{2},j,k\right)+E_x^{n-l/m}\left(i+\frac{1}{2},j,k\right)+\text{ \\ Exn-l/m(i+12,j+1,k)+Exn-l/m(i-12,j+1,k)}\right)J_z^{n+\frac{l}{m}}\left(i,j,k+\frac{1}{2}\right)=e^{{\mathrm{vc}}_l\mathrm{dt}}{J}_z^{n+\frac{l-1}{m}}\left(i,j,k+\frac{1}{2}\right)-\text{ \\ e-vcl⁢dt-1v⁢Ezn+l-1m(i,j,k+12)} \end{gathered}$$

Wherein, i, j, k represent spatial nodes of electrical field, magnetic field, and current density:

${\gamma }_{x1}=\frac{9}{8}\frac{1}{\Delta x},{\gamma }_{y1}=\frac{9}{8}\frac{1}{\Delta y},{\gamma }_{z1}=\frac{9}{8}\frac{1}{\Delta z},{\gamma }_{x2}=\frac{-1}{24}\frac{1}{\Delta x},{\gamma }_{y2}=\frac{-1}{24}\frac{1}{\Delta y},{\gamma }_{z2}=\frac{-1}{24}\frac{1}{\Delta z}.$

Optionally, the electromagnetic model comprises an anisotropic magnetized plasma slab model, a blunt cone model, and a sphere model.

A system for processing anisotropic magnetized plasma medium, including:

• • A Maxwell equation and polarization current density equation determination module, configured to obtain a Maxwell equation and a polarization current density equation based on electromagnetic characteristics of anisotropic magnetized plasma;
  • A Maxwell equation and polarization current density equation processing module, configured to process the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed;
  • A matrix exponential time-domain finite difference processing module, configured to obtain numerical iterative equations for electric field intensity, magnetic field intensity, and polarization current density in anisotropic magnetized plasma medium using the matrix exponential time-domain finite difference method based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed;
  • An electromagnetic characteristic determination module, configured to determine electromagnetic characteristics of the electromagnetic model using a numerical modeling simulation electromagnetic model based on the numerical iterative equations of the electric field intensity, the magnetic field intensity, and the polarization current density in the anisotropic magnetized plasma medium.

Optionally, the Maxwell equation is:

$\nabla \times E=-{\mu }_0\frac{\partial H}{\partial t};$ $\nabla \times H={\epsilon }_0\frac{\partial E}{\partial t}+;$

The polarization current density equation is:

$\begin{array}{cc}\frac{dJ}{dt}+v={\epsilon }_0{\omega }_p^{2}E+& { }_b\times \end{array};$

Wherein, E is a electric field intensity, H is a magnetic field intensity, J is a polarization current density, ε₀ and μ₀ are respectively a dielectric constant and a permeability in vacuum, ω_p is a plasma frequency, v is a plasma collision frequency, ω_b=B₀/m_e is an electron cyclotron frequency, B₀ is an external static magnetic field, e is an electron charge, and m_e is an electron mass.

Optionally, the Maxwell equation and polarization current density equation processing module includes:

• • A Maxwell equation and polarization current density equation processing unit, configured to perform multi-level symplectic discretization on the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density in the symplectic discretization format:

$\begin{array}{c}{E}^{n+l/m}=E^{n+\left(l-1\right)/m}+\frac{d_l\Delta t}{{\epsilon }_0}\nabla \times H^{n+l/m}+\frac{-d_l\Delta t}{{\epsilon }_0}{J}^{n+l/m}\\ H^{n+l/m}=H^{n+\left(l-1\right)/m}+\frac{-c_l\Delta t}{{\mu }_0}\nabla \times E^{n+\left(l-1\right)/m}\\ J^{n+l/m}={\mathrm{DJ}}^{n+\left(l-1\right)/m}+{\epsilon }_0{\omega }_p^2{\mathrm{FE}}^{n+\left(l-1\right)/m}\end{array},$

Wherein, d_l and c_l are coefficients of the symplectic propagator,

$\Omega =\left[\begin{array}{ccc}-v& -{\omega }_{bz}& {\omega }_{by}\\ {\omega }_{bz}& -v& -{\omega }_{bx}\\ -{\omega }_{\mathrm{by}}& {\omega }_{bx}& -v\end{array}\right],$

l is the order of a field component =J, , )^T, and m is the number of stages of non-dissipative p-order symplectic integral, D=e^CM, F=Ω⁻¹(D−I).

Optionally, the matrix exponential time-domain finite difference processing module includes:

• • A Matrix exponential time-domain finite difference processing unit, assuming that the plasma medium is biased by a static magnetic field in the z-direction, based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, the matrix exponential time-domain finite difference method is used to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium:

$E_x^{n+l/m}\left(i+\frac{1}{2},j,k\right)=E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{y1}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{1}{2},k\right)\right]+{\gamma }_{z1}\left[-H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+H_y^{n+l/m}\left(i+\frac{1}{2},j,k-\frac{1}{2}\right)\right]+{\gamma }_{y2}\left[H_2^{n+l/m}\left(i+\frac{1}{2},j+\frac{3}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{3}{2},k\right)\right]+{\gamma }_{z2}\left[-H_y^{n❘l/m}\left(i+\frac{1}{2},j,k+\frac{3}{2}\right)+H_y^{n❘l/m}\left(i+\frac{1}{2},j,k-\frac{3}{2}\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)$ $E_y^{n❘l/m}\left(i,j+\frac{1}{2},k\right)=E_x^{n❘\left(l1\right)/m}\left(i,j+\frac{1}{2},k\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{z1}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)-H_x^{n-1/m}\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)\right]+{\gamma }_{x1}\left[-H_z^{n❘l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)+H_z^{n❘l/m}\left(i-\frac{1}{2},j+\frac{1}{2},k\right)\right]+{\gamma }_{z2}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{3}{2}\right)-H_x^{n+l/m}\left(i,j+\frac{1}{2},k-\frac{3}{2}\right)\right]+{\gamma }_{x2}\left[-H_z^{n+l/m}\left(i+\frac{3}{2},j+\frac{1}{2},k\right)+H_z^{n+l/m}\left(i-\frac{3}{2},j+\frac{1}{2},k\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_y^{n-\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)$ $E_z^{n+l/m}\left(i,j,k+\frac{1}{2}\right)=E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{x1}\left[H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)-H_y^{n+l/m}\left(i-\frac{1}{2},j,k+\frac{1}{2}\right)\right]+{\gamma }_{y1}\left[-H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)+H_x^{n+l/m}\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)\right]+{\gamma }_{x2}\left[H_y^{n+l/m}\left(i+\frac{3}{2},j,k+\frac{1}{2}\right)-H_y^{n+l/m}\left(i-\frac{3}{2},j,k+\frac{1}{2}\right)\right]+{\gamma }_{y2}\left[-H_x^{n+l/m}\left(i,j+\frac{3}{2},k+\frac{1}{2}\right)+H_x^{n+l/m}\left(i,j-\frac{3}{2},k+\frac{1}{2}\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)$ $H_x^{n❘l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)=H_x^{n❘\left(l1\right)/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{z1}\left[E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+1\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)\right]+{\gamma }_{y1}\left[-E_z^{n-\left(l-1\right)/m}\left(i,j+1,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z2}\left[E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+2\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k-1\right)\right]+{\gamma }_{y2}\left[-E_z^{n-\left(l-1\right)/m}\left(i,j+2,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j-1,k+\frac{1}{2}\right)\right]\right\}$ $H_y^{n-l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)=H_y^{n-\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{x1}\left[E_z^{n+\left(l-1\right)/m}\left(i+1,j+\frac{1}{2},k\right)-E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z1}\left[-E_x^{n❘\left(l1\right)/m}\left(i+\frac{1}{2},j,k+1\right)+E_x^{n❘\left(l1\right)/m}\left(i+\frac{1}{2},j,k\right)\right]+{\gamma }_{x2}\left[E_z^{n+\left(l-1\right)/m}\left(i+2,j,k+\frac{1}{2}\right)-E_z^{n+\left(l-1\right)/m}\left(i-1,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z2}\left[-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+2\right)+E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k-1\right)\right]\right\}$ $H_z^{n-l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)=H_z^{n-\left(l-1\right)/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{y1}\left[E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j+1,k\right)-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)\right]+{\gamma }_{x1}\left[-E_y^{n❘\left(l1\right)/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n\left(l1\right)/m}\left(i,j+\frac{1}{2},k\right)\right]+{\gamma }_{y2}\left[E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j+2,k\right)-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j-1,k\right)\right]+{\gamma }_{x2}\left[-E_y^{n+\left(l-1\right)/m}\left(i+2,j+\frac{1}{2},k\right)+E_y^{n+\left(l-1\right)/m}\left(i-1,j+\frac{1}{2},k\right)\right]\right\}$ $J_x^{n+\frac{l}{m}}\left(i+\frac{1}{2},j,k\right)=\left[\cos \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+\left[\sin \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0{\omega }_p^2\left\{\frac{v\left[1-\cos \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}+\frac{{\omega }_b\sin \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}\right\}{E}_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0{\omega }_p^2\left\{\frac{v\sin \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}+\frac{{\omega }_b\left[\cos \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}-1\right]}{v^2+{\omega }_b^2}\right\}{E}_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)$

Because a current density position is the same as an electric field position, therefore, when calculating

$J_x^{n+l/m}\left(i+\frac{1}{2},j,k\right),$

it is necessary to perform interpolation because J, and E, are no value at a point on

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)\mathrm{and}{E}_y^{n-l/m}\left(i+\frac{1}{2},j,k\right),$

because J_y and E_y are no value at point

$\left(i+\frac{1}{2},j,k\right),$

so it needs to be averaged by four adjacent diagonal values to obtain:

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(J_y^{n-l/m}\left(i,j-\frac{1}{2},k\right)+J_y^{n-l/m}\left(i,j+\frac{1}{2},k\right)+J_y^{n-l/m}\left(i+1,j+\frac{1}{2},k\right)+J_y^{n-l/m}\left(i+1,j-\frac{1}{2},k\right)\right)$ $E_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(E_y^{n-l/m}\left(i,j-\frac{1}{2},k\right)+E_y^{n-l/m}\left(i,j+\frac{1}{2},k\right)+E_y^{n-l/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n-l/m}\left(i+1,j-\frac{1}{2},k\right)\right)$ $J_y^{n+\frac{l}{m}}\left(i,j+\frac{1}{2},k\right)=\left[\cos \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_y^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)+\left[\sin \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0{\omega }_p^2\left\{\frac{v\left[1-\cos \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}+\frac{{\omega }_b\sin \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}\right\}{E}_y^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)-{\epsilon }_0{\omega }_p^2\left\{\frac{v\sin \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}+\frac{{\omega }_b\left[\cos \left({\omega }_b{c}_{l}dt\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}-1\right]}{v^2+{\omega }_b^2}\right\}{E}_x^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)$

Because the current density position is the same as the electric field position, therefore, when calculating

$J_y^{n+l/m}\left(i,j+\frac{1}{2},k\right),$

it is necessary to perform interpolation on

$J_x^{n-l/m}\left(i,j+\frac{1}{2},k\right)\mathrm{and}{E}_x^{n+l/m}\left(i,j+\frac{1}{2},k\right),$

because J_x and E_x are no value at a point

$\left(i,j+\frac{1}{2},k\right),$

so it needs to be averaged by four adjacent diagonal values to obtain:

$J_x^{n-l/m}\left(i,j+\frac{1}{2},k\right)=\frac{1}{4}\left(J_x^{n-1/m}\left(i-\frac{1}{2},j,k\right)+J_x^{n-l/m}\left(i+\frac{1}{2},j,k\right)+J_x^{n-l/m}\left(i+\frac{1}{2},j+1,k\right)+J_x^{n-l/m}\left(i-\frac{1}{2},j+1,k\right)\right)$ $E_x^{n+l/m}\left(i,j+\frac{1}{2},k\right)=\frac{1}{4}\left(E_K^{n-l/m}\left(i-\frac{1}{2},j,k\right)+E_x^{n-l/m}\left(i+\frac{1}{2},j,k\right)+E_x^{n-l/m}\left(i+\frac{1}{2},j+1,k\right)+E_x^{n-l/m}\left(i-\frac{1}{2},j+1,k\right)\right)$ $J_z^{n+\frac{l}{m}}\left(i,j,k+\frac{1}{2}\right)=e^{-{\mathrm{vc}}_l\mathrm{dt}}{J}_z^{n+\frac{l-1}{m}}\left(i,j,k+\frac{1}{2}\right)-\frac{e^{-{\mathrm{vc}}_l\mathrm{dt}}-1}{v}{E}_z^{n+\frac{l-1}{m}}\left(i,j,k+\frac{1}{2}\right)$

Wherein, i, j, k represents spatial nodes of electrical field, magnetic field, and current density:

${\gamma }_{x1}=\frac{9}{8}\frac{1}{\Delta x},{\gamma }_{y1}=\frac{9}{8}\frac{1}{\Delta y},{\gamma }_{z1}=\frac{9}{8}\frac{1}{\Delta z},{\gamma }_{x2}=\frac{-1}{24}\frac{1}{\Delta x},{\gamma }_{x2}=\frac{-1}{24}\frac{1}{\Delta y},{\gamma }_{z2}=\frac{-1}{24}\frac{1}{\Delta z}.$

Optionally, the electromagnetic model includes an anisotropic magnetized plasma slab model, a blunt cone model, and a sphere model.

According to the specific embodiments provided by the present disclosure, the following technical effects are disclosed:

• • The present disclosure provides a method for processing anisotropic magnetized plasma medium, which obtains Maxwell equation and polarization current density equation based on the electromagnetic characteristics of each anisotropic magnetized plasma; process the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed; based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, numerical iterative equations for electric field intensity, magnetic field intensity, and polarization current density in anisotropic magnetized plasma medium are obtained using the matrix exponential time-domain finite difference method; according to the numerical iterative equations for the electric field intensity, the magnetic field intensity, and polarization current density in anisotropic magnetized plasma medium, a numerical modeling simulation electromagnetic model is used to determine the electromagnetic characteristics of the electromagnetic model. The use of the present disclosure ME-FDTD (4,4) to simulate the propagation of electromagnetic waves in anisotropic magnetized plasma medium has higher numerical calculation accuracy.

BRIEF DESCRIPTION OF THE DRAWINGS

In order to provide a clearer illustration of the embodiments of the present disclosure or the technical solutions in the prior art, a brief introduction will be given to the accompanying drawings required in the embodiments. It is evident that the accompanying drawings in the following description are only some embodiments of the present disclosure. For ordinary skilled person in the art, other accompanying drawings can also be obtained based on these drawings without any creative effort.

FIG. 1 is a flowchart of a method for processing anisotropic magnetized plasma medium of the present disclosure;

FIG. 2 shows the variation trend of the modulus value of growth factor ζ with knumz;

FIG. 3 is a change curve of dispersion error of RCP waves with electromagnetic wave frequency;

FIG. 4 is a change curve of dispersion error of LCP waves with electromagnetic wave frequency;

FIG. 5 is a geometric schematic diagram of anisotropic magnetized plasma slab;

FIG. 6 is a change curve of the reflection coefficient of RCP waves with frequency;

FIG. 7 is a change curve of the reflection coefficient of LCP waves with frequency;

FIG. 8 is a simulation model of a three-dimensional blunt cone;

FIG. 9 is the RCS result of the three-dimensional blunt cone;

FIG. 10 shows the relative calculation errors of the two methods;

FIG. 11 is the RCS result of a three-dimensional sphere;

FIG. 12 shows the relative calculation errors of the two methods;

FIG. 13 is a structural diagram of the anisotropic magnetized plasma medium processing system of the present disclosure.

DETAILED DESCRIPTION OF THE EMBODIMENTS

The following will provide a clear and complete description of the technical solution in the embodiments of the present disclosure, in conjunction with the accompanying drawings. Obviously, the described embodiments are only a part of the embodiments of the present disclosure, not all of them. Based on the embodiments in the present disclosure, all other embodiments obtained by ordinary skilled person in the art without creative labor fall within the scope of the present disclosure.

The objective of the present disclosure is to provide a method for processing anisotropic magnetized plasma medium and a system thereof, which can improve higher numerical calculation accuracy.

In order to make the above objectives, features, and advantages of the present disclosure more apparent and understandable, the following will provide further detailed explanations of the present disclosure in conjunction with the accompanying drawings and specific implementation methods.

FIG. 1 is a flowchart of anisotropic magnetized plasma medium processing methods of the present disclosure. As shown in FIG. 1, a method for processing anisotropic magnetized plasma medium includes:

• • Step 101: obtaining a Maxwell equation and a polarization current density equation based on electromagnetic characteristics of anisotropic magnetized plasma;
  • The Maxwell equation is:

$\begin{array}{cc}\nabla \times E=-{\mu }_0\frac{\partial H}{\partial t}& \left(1\right)\end{array}$ $\begin{array}{cc}\nabla \times H={\epsilon }_0\frac{\partial E}{\partial t}+& \left(2\right)\end{array}$

The polarization current density equation is:

$\begin{array}{cc}\frac{\mathrm{dJ}}{\mathrm{dt}}+v={\epsilon }_0{\omega }_p^{2}E{+}_b\times & \left(3\right)\end{array}$

• • wherein, E is the electric field intensity, H is the magnetic field intensity, J is the polarization current density, ε₀, μ₀ are respectively a dielectric constant and a permeability in vacuum, ω_p is a plasma frequency, v is a plasma collision frequency, ω_b=B₀/m_e is a electron cyclotron frequency, B₀ is an external static magnetic field, e is an electron charge, and m_e is an electron mass.

The equation can be expressed in the following form

$\begin{array}{cc}\frac{\mathrm{dJ}}{\mathrm{dt}}={\epsilon }_0{\omega }_p^{2}E+\Omega J& \left(4\right)\end{array}$

Wherein

$\begin{array}{cc}\Omega =\left[\begin{array}{ccc}-v& -{\omega }_{\mathrm{bz}}& {\omega }_{\mathrm{by}}\\ {\omega }_{\mathrm{bz}}& -v& -{\omega }_{\mathrm{bx}}\\ -{\omega }_{\mathrm{by}}& {\omega }_{\mathrm{bx}}& -v\end{array}\right]& \left(5\right)\end{array}$

Then, write equations (1), (2), and (3) in matrix form:

$\begin{array}{cc}\frac{\partial }{\partial t}\left[\begin{array}{c}H\\ E\\ J\end{array}\right]=\left[\begin{array}{ccc}0& -\frac{\nabla \times }{{\mu }_0}& 0\\ \frac{\nabla \times }{{\epsilon }_0}& 0& -\frac{1}{{\epsilon }_0}\\ 0& {\epsilon }_0{\omega }_p^2& \Omega \end{array}\right]\left[\begin{array}{c}H\\ E\\ J\end{array}\right]& \left(6\right)\end{array}$

Let =J, , )^T, the equation can be expressed as:

$\begin{array}{cc}\frac{\partial \Psi }{\partial t}=& \left(7\right)\end{array}$

Wherein

$\begin{array}{cc}A=\left[\begin{array}{ccc}0& -\frac{\nabla \times }{{\mu }_0}& 0\\ \frac{\nabla \times }{{\epsilon }_0}& 0& -\frac{1}{{\epsilon }_0}\\ 0& {\epsilon }_0{\omega }_p^2& \Omega \end{array}\right],\nabla \times =\left[\begin{array}{ccc}0& -\frac{\partial }{\partial z}& \frac{\partial }{\partial y}\\ \frac{\partial }{\partial z}& 0& -\frac{\partial }{\partial x}\\ -\frac{\partial }{\partial y}& \frac{\partial }{\partial x}& 0\end{array}\right]& \left(8\right)\end{array}$

• • Step 102: processing the Maxwell equation and the polarization current density equation to obtain an electric field intensity, a magnetic field intensity, and a polarization current density after processed.
  • performing multi-level symplectic discretization on the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density in a symplectic discretization format:

$E^{n+l/m}=E^{n+\left(l-1\right)/m}+\frac{d_l\Delta t}{{\epsilon }_0}\nabla \times H^{n+l/m}+\frac{-d_l\Delta t}{{\epsilon }_0}{J}^{n+l/m}$ $H^{n+l/m}=H^{n+\left(l-1\right)/m}+\frac{-c_l\Delta t}{{\mu }_0}\nabla \times E^{n+\left(l-1\right)/m}$ $J^{n+l/m}={\mathrm{DJ}}^{n+\left(l-1\right)/m}+{\epsilon }_0{\omega }_p^2{\mathrm{FE}}^{n+\left(l-1\right)/m},$

• • wherein, d_l and c_l are coefficients of symplectic propagator,

$\Omega =\left[\begin{array}{ccc}-v& -{\omega }_{\mathrm{bz}}& {\omega }_{\mathrm{by}}\\ {\omega }_{\mathrm{bz}}& -v& -{\omega }_{\mathrm{bx}}\\ -{\omega }_{\mathrm{by}}& {\omega }_{\mathrm{bx}}& -v\end{array}\right],$

• •  l is the order of a field component =J, , )^T, and m is the number of stages of non-dissipative p-order symplectic integral, D=e^CM, F=Ω⁻¹(D−I).

According to the theory of symplectic propagation operators, the field at the time t=Δt can be represented by an exponential operator exp(AΔt) as

$\begin{array}{cc}\Psi \left(\Delta t\right)=e^{A\Delta t}\psi \left(0\right)& \left(9\right)\end{array}$

Let

$\begin{array}{cc}A=U+V& \left(10\right)\end{array}$

Wherein

$\begin{array}{cc}U=\left[\begin{array}{ccc}0& -\frac{\nabla \times }{{\mu }_0}& 0\\ 0& 0& 0\\ 0& {\epsilon }_0{\omega }_p^2& \Omega \end{array}\right],V=\left[\begin{array}{ccc}0& 0& 0\\ \frac{\nabla \times }{{\epsilon }_0}& 0& -\frac{1}{{\epsilon }_0}\\ 0& 0& 0\end{array}\right]& \left(11\right)\end{array}$

Using a non-dissipative m-stage and p-order exponential integral method to approximate the time evolution matrix exp (AΔt),

$\begin{array}{cc}{e}^{A\Delta t}=e^{\left(U+V\right)\Delta t}=\prod _{l=1}^m{e}^{d_l\Delta \mathrm{tV}}{e}^{c_l\Delta \mathrm{tU}}+O\left({\Delta }_t^{p+1}\right)& \left(12\right)\end{array}$

Wherein, d_l and c_l are the symplectic propagator coefficients, and it can be observed that U^α0 V^α=0 (α≥2). Perform Taylor series expansion to e^dlΔtV in equation (12) to obtain:

$\begin{array}{cc}{e}^{d_l\Delta \mathrm{tV}}=I+d_l\Delta \mathrm{tV}+\frac{{\left(d_l\Delta \mathrm{tV}\right)}^2}{2!}++\frac{{\left(d_l\Delta \mathrm{tV}\right)}^n}{n!}+=I+d_l\Delta \mathrm{tV}& \left(13\right)\end{array}$

Perform Taylor series expansion to e^clΔtU in equation (12) to obtain:

$\begin{array}{cc}{e}^{c_l\Delta \mathrm{tU}}=I+c_l\Delta \mathrm{tU}+\frac{{\left(c_l\Delta \mathrm{tU}\right)}^2}{2!}++\frac{{\left(c_l\Delta \mathrm{tU}\right)}^n}{n!}+=G& \left(14\right)\end{array}$

According to the form of matrix U in equation (11), it can be written as follows:

$\begin{array}{cc}U=\left[\begin{array}{ccc}0& u_{12}& 0\\ 0& 0& 0\\ 0& u_{32}& u_{33}\end{array}\right]& \left(15\right)\end{array}$

Then

$\begin{array}{cc}{U}^2=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& u_{32}{u}_{33}& u_{33}^2\end{array}\right]& \left(16\right)\end{array}$ $$ $U^n=\left[\begin{array}{ccc}0& 0& 0\\ 0& 0& 0\\ 0& u_{32}{u}_{33}^{n‐1}& u_{33}^n\end{array}\right]$

Except for zero elements in matrix G, the other elements are:

$\begin{array}{cc}{b}_{11}=b_{22}=1& \left(17\right)\end{array}$ $\begin{array}{cc}{b}_{12}=\alpha u_{12}& \left(18\right)\end{array}$ $\begin{array}{cc}\begin{array}{c}{b}_{32}=\alpha u_{32}+\frac{{\alpha }^2{u}_{32}{u}_{33}}{2!}++\frac{{\alpha }^n{u}_{32}{u}_{33}^{n-1}}{n!}\\ =\frac{u_{32}}{u_{33}}\left(1+\alpha u_{33}+\frac{{\alpha }^2{u}_{33}}{2!}++\frac{{\alpha }^n{u}_{33}^{n-1}}{n!}-1\right)\\ =\frac{u_{32}}{u_{33}}\left(e^{\alpha u^{33}}-1\right)\end{array}& \left(19\right)\end{array}$ $\begin{array}{cc}{b}_{33}=1+\alpha u_{33}+\frac{{\alpha }^2{u}_{33}}{2!}++\frac{{\left(\alpha u_{33}\right)}^n}{n!}=e^{{\mathrm{\alpha \iota \iota }}_{33}}& \left(20\right)\end{array}$

Therefore, according to the Taylor series expansion form of the exponential matrix e^dlΔtV and e^clΔtU (equations (13) and (14)), the symplectic discretization form of the field component =J, , )^T in order l is:

$\begin{array}{cc}{E}^{n+l/m}=E^{n+\left(l-1\right)/m}+\frac{d_l\Delta t}{{\epsilon }_0}\nabla \times H^{n+l/m}+\frac{-d_l\Delta t}{{\epsilon }_0}{J}^{n+l/m}& \left(21\right)\end{array}$ $\begin{array}{cc}{H}^{n+l/m}=H^{n+\left(l-1\right)/m}+\frac{-c_l\Delta t}{{\mu }_0}\nabla \times E^{n+\left(l-1\right)/m}& \left(22\right)\end{array}$ $\begin{array}{cc}{J}^{n+l/m}=e^{c_l\mathrm{\Omega \Delta }t}{J}^{n+\left(l-1\right)/m}+{\epsilon }_0{\omega }_p^2{\Omega }^{-1}\left(e^{c_l\mathrm{\Omega \Delta }t}-I\right)E^{n+\left(l-1\right)/m}& \left(23\right)\end{array}$

Equations (21), (22), and (23) represent the symplectic discretization form of anisotropic magnetized plasma in a symplectic framework. equations (21) and (22) can be directly used for numerical iteration to solve E and H. But for equation (23), it can be seen from equation (5) that Ω is a 3×3 matrix, and further processing of the exponential matrix operator e^clΩΔt is required, to obtain its equivalent 3×3 matrix form before obtaining the specific numerical iteration equation for equation (23).

For the exponential matrix coefficients e^clΩΔt, let C=c_lΩ, then according to equation (5), it can be obtained that

$\begin{array}{cc}C=\left[\begin{array}{ccc}-{\mathrm{vc}}_l& -{\omega }_{\mathrm{bz}}{c}_l& {\omega }_{\mathrm{by}}{c}_l\\ {\omega }_{\mathrm{bz}}{c}_l& -{\mathrm{vc}}_l& -{\omega }_{\mathrm{bx}}{c}_l\\ -{\omega }_{\mathrm{by}}{c}_l& {\omega }_{bx}{c}_l& -{\mathrm{vc}}_l\end{array}\right]& \left(24\right)\end{array}$

Due to the fact that the elements in matrix C are all non-zero, it is difficult to obtain the equivalent polynomial function of the exponential matrix operator e^Ω by Taylor series expansion.

The matrix index e^Ct can be expanded into the following form:

$\begin{array}{cc}{e}^{Ct}=x_1\left(t\right)I+x_2\left(t\right)C+x_3\left(t\right)C^2& \left(25\right)\end{array}$

Wherein, λ₁, λ₂, and λ₃ are the characteristic value of the matrix C,

$\begin{array}{cc}{\lambda }_1=-{\mathrm{vc}}_l& \left(26\right)\end{array}$ $\begin{array}{cc}{\lambda }_2=-{\mathrm{vc}}_l+{c_l\left(-{\omega }_{\mathrm{bx}}^2-{\omega }_{\mathrm{by}}^2-{\omega }_{\mathrm{bz}}^2\right)}^{1/2}& \left(27\right)\end{array}$ $\begin{array}{cc}{\lambda }_3=-{\mathrm{vc}}_l-{c_l\left(-{\omega }_{\mathrm{bx}}^2-{\omega }_{\mathrm{by}}^2-{\omega }_{\mathrm{bz}}^2\right)}^{1/2}& \left(28\right)\end{array}$

Since λ₁, λ₂, and λ₃ are not equal to each other, the general solution can be written as:

$\begin{array}{cc}{x}_i\left(t\right)=a_{i1}{e}^{{\lambda }_{1}t}+a_{i2}{e}^{{\lambda }_{2}t}+a_{i3}{e}^{{\lambda }_{3}t}\infty & \left(29\right)\end{array}$

Wherein, a_i1, a_i2 and a_i3 (i=1, 2, 3) represent coefficients. In addition, the general solution satisfies the following initial conditions:

• • (1) x(0)=1, x′(0)=0, x″(0)=0, substitute it into equation (29) to obtain

$\begin{array}{cc}{x}_1\left(t\right)=a_{11}{e}^{{\lambda }_{1}t}+a_{12}{e}^{{\lambda }_{2}t}+a_{13}{e}^{{\lambda }_{3}t}& \left(30\right)\end{array}$

Wherein

$\begin{array}{cc}{a}_{11}=\frac{{\lambda }_2{\lambda }_3}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_1-{\lambda }_3\right)}& \left(31\right)\end{array}$ $\begin{array}{cc}{a}_{12}=\frac{-{\lambda }_1{\lambda }_3}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_2-{\lambda }_3\right)}& \left(32\right)\end{array}$ $\begin{array}{cc}{a}_{13}=\frac{{\lambda }_1{\lambda }_2}{{\lambda }_1{\lambda }_2-{\lambda }_1{\lambda }_3-{\lambda }_2{\lambda }_3+{\lambda }_3^2}& \left(33\right)\end{array}$

• • (2) x(0)=0, x′(0)=1, x″(0)=0, then obtain:

$\begin{array}{cc}{x}_2\left(t\right)=a_{21}{e}^{{\lambda }_{1}t}+a_{22}{e}^{{\lambda }_{2}t}+a_{23}{e}^{{\lambda }_{3}t}& \left(34\right)\end{array}$

Wherein

$\begin{array}{cc}{a}_{21}=\frac{-{\lambda }_2-{\lambda }_3}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_1-{\lambda }_3\right)}& \left(35\right)\end{array}$ $\begin{array}{cc}{a}_{22}=\frac{{\lambda }_1+{\lambda }_3}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_2-{\lambda }_3\right)}& \left(36\right)\end{array}$ $\begin{array}{cc}{a}_{23}=\frac{-{\lambda }_1-{\lambda }_2}{{\lambda }_1{\lambda }_2-{\lambda }_1{\lambda }_3-{\lambda }_2{\lambda }_3+{\lambda }_3^2}& \left(37\right)\end{array}$

• • (3) x(0)=0, x′(0)=0, x″(0)=1, then obtain:

$\begin{array}{cc}{x}_3\left(t\right)=a_{31}{e}^{{\lambda }_{1}t}+a_{32}{e}^{{\lambda }_{2}t}+a_{33}{e}^{{\lambda }_{3}t}& \left(38\right)\end{array}$

Wherein

$\begin{array}{cc}{a}_{31}=\frac{1}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_1-{\lambda }_3\right)}& \left(39\right)\end{array}$ $\begin{array}{cc}{a}_{32}=\frac{-1}{\left({\lambda }_1-{\lambda }_2\right)\left({\lambda }_2-{\lambda }_3\right)}& \left(40\right)\end{array}$ $\begin{array}{cc}{a}_{33}=\frac{1}{{\lambda }_1{\lambda }_2-{\lambda }_1{\lambda }_3-{\lambda }_2{\lambda }_3+{\lambda }_3^2}& \left(41\right)\end{array}$

Substitute equations (26)-(41) into equation (25), and the exponential matrix e^Ct can be represented as a 3×3 matrix. Due to the complexity representation of the matrix e^Ct and its ease of calculation, it is not provided in this disclosure for brevity. Next, change t to Δt, let D=e^CΔt, and then equation (25) becomes:

$\begin{array}{cc}D=e^{C\Delta t}=x_1\left(\Delta t\right)I+x_2\left(\Delta t\right)C+x_3\left(\Delta t\right)C^2& \left(42\right)\end{array}$

Substitute equation (42) into equation (23) and let F=Ω⁻¹(D−I), then the symplectic discretization form of the current source J is:

$\begin{array}{cc}{J}^{n+l/m}={\mathrm{DJ}}^{n+\left(l-1\right)/m}+{\epsilon }_0{\omega }_p^2{\mathrm{FE}}^{n+\left(l-1\right)/m}& \left(43\right)\end{array}$

Finally, equations (21), (22), and (43) represent the symplectic iterative form of electromagnetic field components in anisotropic magnetized plasma medium under multi-level symplectic integration approximation.

• • Step 103: using a matrix exponential time-domain finite difference method to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed.

For the convenience of description, taking the x-direction as an example, the iterative equation of the ME-SFDTD (4,4) method for the electromagnetic field components in the x-direction involved is (spatial partial derivatives are approximated using fourth-order difference):

$\begin{array}{cc}\begin{array}{c}{E}_x^{n+l/m}\left(i+\frac{1}{2},j,k\right)=E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \\ \{{\gamma }_{y1}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{1}{2},k\right)\right]+\\ {\gamma }_{z1}\left[-H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+H_y^{n+l/m}\left(i+\frac{1}{2},j,k-\frac{1}{2}\right)\right]+\\ {\gamma }_{y2}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{3}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{3}{2},k\right)\right]+\\ {\gamma }_{z2}\left[-H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{3}{2}\right)+H_y^{n+l/m}\left(i+\frac{1}{2},j,k-\frac{3}{2}\right)\right]\}-\\ \frac{d_l\Delta t}{{\epsilon }_0}{J}_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)\end{array}& \left(11\right)\end{array}$ $\begin{array}{cc}\begin{array}{c}{J}_x^{n+\frac{l}{m}}\left(i+\frac{1}{2},j,k\right)=\\ \left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+\\ \left[\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+\\ \frac{v\left[1-\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}{E}_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)-\\ \frac{v\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}{E}_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)\end{array}& \left(45\right)\end{array}$ $\begin{array}{cc}\begin{array}{c}{H}_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)=H_x^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)+\frac{c_l\Delta t}{{\mu }_0}\times \\ \{{\gamma }_{z1}\left[H_z^{n+l/m}\left(i,j+\frac{1}{2},k+1\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)\right]+\\ {\gamma }_{y1}\left[-H_y^{n+l/m}\left(i,j+1,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+\\ {\gamma }_{z2}\left[H_z^{n+l/m}\left(i,j+\frac{1}{2},k+2\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k-1\right)\right]+\\ {\gamma }_{y2}\left[-H_y^{n+l/m}\left(i,j+2,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j-1,k+\frac{1}{2}\right)\right]\}\end{array}& \left(46\right)\end{array}$

Wherein,

${\gamma }_{y1}=\frac{9}{8}\frac{1}{\Delta y},{\gamma }_{y2}=\frac{-1}{24}\frac{1}{\Delta y},{\gamma }_{z1}=\frac{9}{8}\frac{1}{\Delta z}\mathrm{and}{\gamma }_{z2}=\frac{-1}{24}\frac{1}{\Delta z},$

the numerical discretization forms of field components in other direction are similar.

It should be noted that, the current density position is the same as the electric field position. Therefore, when calculating

$J_x^{n+l/m}\left(i+\frac{1}{2},j,k\right),$

it is necessary to perform interpolation on

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)\mathrm{and}{E}_y^{n-l/m}\left(i+\frac{1}{2},j,k\right),$

because J_y and E_y are no value at the spatial coordinate

$\left(i+\frac{1}{2},j,k\right),$

so it is necessary to obtain an average weighted approximation from the values of four adjacent field components:

$\begin{array}{cc}{J}_y^{n+l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(J_y^{n+l/m}\left(i,j-\frac{1}{2},k\right)+J_y^{n+l/m}\left(i,j+\frac{1}{2},k\right)+J_y^{n+l/m}\left(i+1,j+\frac{1}{2},k\right)+J_y^{n+l/m}\left(i+1,j-\frac{1}{2},k\right)\right)& \left(47\right)\end{array}$ $\begin{array}{cc}{E}_y^{n+l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(E_y^{n+l/m}\left(i,j-\frac{1}{2},k\right)+E_y^{n+l/m}\left(i,j+\frac{1}{2},k\right)+E_y^{n+l/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n+l/m}\left(i+1,j-\frac{1}{2},k\right)\right)& \left(48\right)\end{array}$

For the convenience of description, consider the numerical stability in one-dimensional cases. When simulating anisotropic magnetized plasma medium, define each field component in the spatial domain as

$\begin{array}{cc}{\varphi }_s^n={\varphi }_0{\zeta }^n{\widehat{\exp (\mathrm{jk}}}_{\mathrm{num}}s\Delta )& \left(49\right)\end{array}$

Wherein, j=√{square root over (−1)}, represents the modulus of the field component, ζ represents the growth factor, n is the time factor, s is the spatial factor, k_num represents the numerical wavenumber, and Δ represents the spatial step size.

By substituting (49) into equations (21), (22), and (23), the iterative equations for the field components of anisotropic magnetized plasma medium can be obtained (from time n to n+l/m):

$\begin{array}{cc}{L1}_{n+l/m}^{}=M_{R1}^n& \left(50\right)\end{array}$

Wherein Ψ=[H_x,H_y,E_x,E_y,J_x,J_y]^T.

According to equation (50), the numerical iteration equation for the m-stage and p-order ME-SFDTD method from n to n+1 time steps can be written in the following matrix form:

$\begin{array}{cc}\prod _{l=0}^{m-1}{L\left(m-l\right)}_{n+1}^{}={\prod _{l=0}^{m-l}}_{R\left(m-l\right)}^n& \left(51\right)\end{array}$ $\begin{array}{cc}{M}_{\mathrm{Li}}=\left[\begin{array}{cccccc}{\zeta }^{c_i}& 0& 0& 0& 0& 0\\ 0& {\zeta }^{\mathrm{ci}}& 0& 0& 0& 0\\ 0& d_i{F}_2{D}_z{\zeta }^{c_i}& {\zeta }^{d_i}& 0& d_i{F}_3{\zeta }^{c_i}& 0\\ -d_i{F}_2{D}_z{\zeta }^{c_i}& 0& 0& {\zeta }^{d_i}& 0& d_i{F}_3{\zeta }^{c_i}\\ 0& 0& 0& 0& {\zeta }^{c_i}& 0\\ 0& 0& 0& 0& 0& {\zeta }^{c_i}\end{array}\right]& \left(52\right)\end{array}$ $\begin{array}{cc}{M}_{\mathrm{Ri}}=\left[\begin{array}{cccccc}1& 0& 0& c_i{F}_1{D}_z& 0& 0\\ 0& 1& -c_i{F}_1{D}_z& 0& 0& 0\\ 0& 0& 1& 0& 0& 0\\ 0& 0& 0& 1& 0& 0\\ 0& 0& \epsilon {\omega }_p^2{P}_i& -\epsilon {\omega }_p^2{Q}_i& A_i& -B_i\\ 0& 0& \epsilon {\omega }_p^2{Q}_i& \epsilon {\omega }_p^2{P}_i& B_i& A_i\end{array}\right]& \left(53\right)\end{array}$

Wherein,

$F_1=\frac{\Delta t}{\mathrm{\mu \Delta }z},F_2=\frac{\Delta t}{\mathrm{\epsilon \Delta }z},F_3=\frac{\Delta t}{\epsilon },A_i=\cos \left({\omega }_b{c}_i\mathrm{dt}\right)e^{-{\mathrm{vc}}_i\mathrm{dt}},B_i=\sin \left({\omega }_b{c}_i\mathrm{dt}\right)e^{-{\mathrm{vc}}_i\mathrm{dt}},D_z=\widehat{\frac{j}{4}}\left[9\sin \left(\frac{k_{\mathrm{num}}\Delta }{2}\right)-\frac{1}{3}\sin \left(\frac{3k_{\mathrm{num}}\Delta }{2}\right)\right],j=\sqrt{-1},P_i=\frac{v\left[1-\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}+\frac{{\omega }_b\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2},Q_i=\frac{\nu \sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{{\nu }^2+{\omega }_b^2}+\frac{{\omega }_b\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}-1\right]}{v^2+{\omega }_b^2},$

i (i=1, 2, 3, . . . , m) stands for the sub-stage.

Substitute equations (52) and (53) into equation (51) to obtain:

$\begin{array}{cc}{n}^{}={\left(\prod _{l=0}^{m-1}{M}_{L\left(m-l\right)}-{\prod _{l=0}^{m-1}}_{R\left(m-l\right)}\right)}^n=0& \left(54\right)\end{array}$

To ensure that equation (54) has a non-zero solution, the determinant of its coefficient matrix M must be 0 (the expression of the coefficient matrix M is easily derived). Meanwhile, in order to ensure the numerical stability of the SFDTD (4,4) method, the modulus of the growth factor ζ must be less than or equal to 1. Although the expression for the growth factor ζ can be obtained by solving the determinant of matrix M, the expression for the growth factor ζ is complex and difficult to determine value range of ζ through the expression. Another feasible solution is to substitute the values of the parameters into the matrix M to calculate the numerical size of ζ. According to equations (7) and (8), it is known that the parameters in matrix M contain ω_b, ω_p, Δz, Δt, v, and σ₂. In numerical calculations, let Δz=75 μm, Δt=0.125 ps, v=20 GHZ, ω_b=100 Grad/s, and ω_p=15.53 Grad/s. Due to

$\widehat{{\sigma }_z=\frac{j}{4}}\left(9\sin \left(\frac{k_{\mathrm{num}}\Delta z}{2}\right)-\frac{1}{3}\sin \left(\frac{3k_{\mathrm{num}}\Delta z}{2}\right)\right),$

considering the symmetry of trigonometric functions, the range of values for k_numΔz is from 0 to π.

As shown in FIG. 2, at kΔz=[0, π], The modulus of c is within and on the unit circle, indicating that the ME-SFDTD (4,4) method proposed in the present disclosure can maintain stability in numerical iterations.

When calculating the dispersion error of the numerical method, it is necessary to obtain the value of the numerical wavenumber k_num, and the calculation process is similar to the solving process in part A. Let

$\begin{array}{cc}\zeta =e^{j\omega \Delta t}& \left(55\right)\end{array}$

Substitute equation (55) into (54) and let det(M)=0, then the expression for the numerical wavenumber k_num of the SFDTD (4,4) method can be solved (The parameter values of ω_b, ω_p, Δz, Δt and v are the same as the settings in the stability analysis section.

The analytical dispersion relationship between left-handed circularly polarized waves and right-handed circularly polarized waves in magnetized plasma is as follows:

$\begin{array}{cc}{\left(\frac{k_{R,L}c}{\omega }\right)}^2=1-\frac{{\omega }_p^2}{{\omega }^2\left[\left(1-\mathrm{jv}/\omega \right)\pm \left({\omega }_b/\omega \right)\right]}& \left(56\right)\end{array}$

Wherein, K_R,L is the complex wave number, and the subscripts R and L represent the right-handed circularly polarized waves (RCP) and the left-handed circularly polarized waves (LCP), respectively. The dispersion error is defined as

$\begin{array}{cc}\mathrm{Re}\left(k_{\mathrm{num}}-k_{R,L}\right)/\mathrm{Re}\left(k_{R,L}\right)& \left(57\right)\end{array}$

FIG. 3 and FIG. 4 respectively show the relationship between the dispersion error of FDTD (2,2) method and SFDTD (4,4) method and the frequency of electromagnetic waves, showing the numerical dispersion errors of the right-handed circularly polarized waves and the left-handed circularly polarized waves. It can be seen that the SFDTD (4,4) method has lower dispersion error compared to the FDTD (2,2) method.

• • Step 104: applying a numerical modeling simulation electromagnetic model to determine electromagnetic characteristics of the electromagnetic model according to the numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in the anisotropic magnetized plasma medium. The electromagnetic model includes an anisotropic magnetized plasma slab model, a blunt cone model, and a sphere model.

For the plasma slab model:

• • The reflection characteristics of anisotropic magnetized plasma slab are calculated through FDTD (2,2) method, and ME-SFDTD (4,4) method using analytical method. As shown in FIG. 4, transverse electromagnetic waves propagate in the z-direction. The simulation area contains 400×9×9 Yee grids, with the plasma slab occupying 200 grids in the x-direction. The size of the grid is Δx=Δy=Δz=75 μm. The calculation area is set as periodic boundary conditions in the y-direction and z-direction, and CPML technology is used to truncate at both ends of the x-direction. The parameters of magnetized plasma are: v=20 GHZ, ω_b=100 Grad/s, ω_p=15.53 Grad/s. The incident electromagnetic wave source adopts the parameters introduced in reference [1], wherein

$\begin{array}{cc}\mathrm{Gaussian}\left(t\right)=\exp \left(-\frac{4{\pi \left(t-t_0\right)}^2}{{\tau }^2}\right)& \left(58\right)\end{array}$

Wherein, τ=150Δt, Δt=0.125 ps, t₀=0.87.

FIG. 6 and FIG. 7 respectively show the reflection coefficients of the right-handed circularly polarized waves and the left-handed circularly polarized waves. It can be clearly observed that compared to the FDTD (2,2) method, the results obtained by the ME-SFDTD (4,4) method are closer to the analytical method, verifying the correctness of the proposed method.

For blunt cone model:

• • Generally speaking, the tip of a high-altitude flying target is similar to a blunt cone, so the anisotropic magnetized plasma blunt cone is modeled and simulated, and its single station RCS is calculated. As shown in FIG. 8, the blunt cone is placed in the total field area, with the adjacent area being the scattering field area and the outermost layer being the absorption layer. In order to better absorb scattered electromagnetic waves, a convolutional perfectly matched layer is used as the absorbing boundary condition. The simulation area contains 80×80×80 grids, with a grid size of Δx=Δy=Δz=0.25 cm. The bottom radius of the blunt cone model is 0.75 cm and the length of the blunt cone model is 7.5 cm. The parameters of magnetized plasma are: v=20 GHZ, ω_b=20 Grad/s, and ω_p=180.32 Grad/s. The Gaussian plane waves are introduced by the total field/scattered field technique, and the time-domain form of the Gaussian waveform is the same as formula (58), wherein τ=80Δt, Δt=4.167 ps, t₀=0.87. In numerical simulation, in order to compare and verify the correctness of the proposed method, the RCS results calculated using the fine grid FDTD (2,2) method are used as the benchmark. FIG. 9 shows the RCS results calculated using several different numerical methods. It can be observed that compared to the coarse grid FDTD (2,2) method, the ME-SFDTD (4,4) method results are closer to the coarse grid FDTD (2,2) method. In addition, in order to quantitatively analyze the calculation error, formula (59) is used to calculate the calculation error of numerical methods:

$\begin{array}{cc}{E}_r\left(f\right)=\frac{❘{E\left(f\right)}_{\mathrm{num}}-{E\left(f\right)}_{\mathrm{Dense}\mathrm{Grid}}❘}{\max ❘{E\left(f\right)}_{\mathrm{Dense}\mathrm{Grid}}❘}\times 100\%& \left(59\right)\end{array}$

Wherein, E(f)_num represents the results obtained by the coarse grid FDTD (2,2) method and the coarse grid ME-SFDTD (4,4) method, while E(f)_DenseGrid represents the results obtained by the fine grid FDTD (2,2) method. FIG. 10 shows the relative numerical calculation errors of the two methods, and it can be observed that the ME-SFDTD (4,4) method has lower numerical calculation errors, once again verifying the correctness of the proposed method.

For sphere model:

• • Finally, the RCS of the plasma sphere is also calculated. In numerical simulation, the settings for time step, space step, plasma parameters, and calculation space size are the same as in Example B. Just replace the blunt cone in FIG. 8 with a sphere, with a radius of R=5 cm. FIG. 11 and FIG. 12 show the RCS calculation results and relative numerical calculation errors, respectively. It can be observed that the proposed ME-SFDTD (4,4) method provides more accurate results compared to the traditional FDTD (2,2) method.

The present disclosure combines the matrix exponential method and the SFDTD (4,4) method with high order numerical accuracy in both time and space to simulate the propagation characteristics of electromagnetic waves in anisotropic magnetized plasma medium. Firstly, a unified matrix form and a time multi-level symplectic discretization process of the Maxwell equation and the polarization current density equation are established. Then, the exponential matrix method is used to solve the matrix exponential coefficient terms included in the derived discrete equation. Finally, combined with the spatial fourth-order difference approximation method, a numerical iteration equation for the ME-SFDTD (4,4) method for processing anisotropic magnetized plasma medium is derived. In numerical simulation, compared to the traditional FDTD (2,2) method, the ME-SFDTD (4,4) method can more accurately simulate the electromagnetic characteristics of the anisotropic magnetized plasma slab, the blunt cone, and the sphere, fully verifying the correctness of this method.

FIG. 13 shows the structural diagram of the anisotropic magnetized plasma medium processing system of the present disclosure. As shown in FIG. 13, a anisotropic magnetized plasma medium processing system includes:

• • A system for processing anisotropic magnetized plasma medium, comprising:
  • A Maxwell equation and polarization current density equation determination module 201, configured to obtain a Maxwell equation and a polarization current density equation based on electromagnetic characteristics of anisotropic magnetized plasma;
  • A Maxwell equation and polarization current density equation processing module 202, configured to process the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed;
  • A matrix exponential time-domain finite difference processing module 203, configured to obtain numerical iterative equations for electric field intensity, magnetic field intensity, and polarization current density in anisotropic magnetized plasma medium using the matrix exponential time-domain finite difference method based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed; and

An electromagnetic characteristic determination module 204, configured to determine electromagnetic characteristics of the electromagnetic model using a numerical modeling simulation electromagnetic model based on the numerical iterative equations of the electric field intensity, the magnetic field intensity, and the polarization current density in the anisotropic magnetized plasma medium.

The Maxwell equation in step 201 is:

$\nabla \times E=-{\mu }_0\frac{\partial H}{\partial t};\nabla \times H={\epsilon }_0\frac{\partial E}{\partial t}+;$

• • the polarization current density equation is:

$\frac{\mathrm{dJ}}{\mathrm{dt}}+v={\epsilon }_0{\omega }_p^{2}E{+}_b\times ;$

• • wherein, E is a electric field intensity, H is a magnetic field intensity, J is a polarization current density, ε₀ and μ₀ are respectively a dielectric constant and a permeability in vacuum, ω_p is a plasma frequency, v is a plasma collision frequency, ω_b=B₀/m_e is an electron cyclotron frequency, B₀ is an external static magnetic field, e is an electron charge, and m_e is an electron mass.

The Maxwell equation and polarization current density equation processing module 202 specifically includes:

• • a Maxwell equation and polarization current density equation processing unit, configured to perform multi-level symplectic discretization on the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density in the symplectic discretization format:

$E^{n+l/m}=E^{n+\left(l-1\right)/m}+\frac{d_l\Delta t}{{\epsilon }_0}\nabla \times H^{n+l/m}+\frac{-d_l\Delta t}{{\epsilon }_0}{J}^{n+l/m}{H}^{n+l/m}=H^{n+\left(l-1\right)/m}+\frac{-c_l\Delta t}{{\mu }_0}\nabla \times E^{n+\left(l-1\right)/m}{J}^{n+l/m}={\mathrm{DJ}}^{n+\left(l-1\right)/m}+{\epsilon }_0{\omega }_p^2{\mathrm{FE}}^{n+\left(l-1\right)/m},$

• • wherein, d_l and c_l are coefficients of the symplectic propagator,

$\Omega =\left[\begin{array}{ccc}-v& -{\omega }_{\mathrm{bz}}& {\omega }_{\mathrm{by}}\\ {\omega }_{\mathrm{bz}}& -v& -{\omega }_{\mathrm{bx}}\\ -{\omega }_{\mathrm{by}}& {\omega }_{\mathrm{bx}}& -v\end{array}\right],$

• •  l is the order of a field component =J, , )^T, and m is the number of stages of non-dissipative p-order symplectic integral, D=e^CM, F=Ω⁻¹(D−I).

The matrix exponential time-domain finite difference processing module 203 specifically includes:

• • a Matrix exponential time-domain finite difference processing unit, assuming that the plasma medium is biased by a static magnetic field in the z-direction, based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, the matrix exponential time-domain finite difference method is used to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium:

$E_x^{n+l/m}\left(i+\frac{1}{2},j,k\right)=E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{y1}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{1}{2},k\right)\right]+{\gamma }_{z1}\left[-H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+H_y^{n+l/m}\left(i+\frac{1}{2},j,k-\frac{1}{2}\right)\right]+{\gamma }_{y2}\left[H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{3}{2},k\right)-H_z^{n+l/m}\left(i+\frac{1}{2},j-\frac{3}{2},k\right)\right]+{\gamma }_{z2}\left[-H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{3}{2}\right)+H_y^{n+l/m}\left(i+\frac{1}{2},j,k-\frac{3}{2}\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)$ $E_y^{n+l/m}\left(i,j+\frac{1}{2},k\right)=E_x^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{z1}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)-H_x^{n-l/m}\left(i,j+\frac{1}{2},k-\frac{1}{2}\right)\right]+{\gamma }_{x1}\left[-H_z^{n+l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)+H_z^{n+l/m}\left(i-\frac{1}{2},j+\frac{1}{2},k\right)\right]+{\gamma }_{z2}\left[H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{3}{2}\right)-H_x^{n+l/m}\left(i,j+\frac{1}{2},k-\frac{3}{2}\right)\right]+{\gamma }_{x2}\left[-H_z^{n+l/m}\left(i+\frac{3}{2},j+\frac{1}{2},k\right)+H_z^{n+l/m}\left(i-\frac{3}{2},j+\frac{1}{2},k\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_y^{n-\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)$ $E_z^{n+l/m}\left(i,j,k+\frac{1}{2}\right)=E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)+\frac{d_l\Delta t}{{\epsilon }_0}\times \left\{{\gamma }_{x1}\left[H_y^{n+l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)-H_y^{n+l/m}\left(i-\frac{1}{2},j,k+\frac{1}{2}\right)\right]+{\gamma }_{y1}\left[-H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)+H_x^{n+l/m}\left(i,j-\frac{1}{2},k+\frac{1}{2}\right)\right]+{\gamma }_{x2}\left[H_y^{n+l/m}\left(i+\frac{3}{2},j,k+\frac{1}{2}\right)-H_y^{n+l/m}\left(i-\frac{3}{2},j,k+\frac{1}{2}\right)\right]+{\gamma }_{y2}\left[-H_x^{n+l/m}\left(i,j+\frac{3}{2},k+\frac{1}{2}\right)+H_x^{n+l/m}\left(i,j-\frac{3}{2},k+\frac{1}{2}\right)\right]\right\}-\frac{d_l\Delta t}{{\epsilon }_0}{J}_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)$ $H_x^{n+l/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)=H_x^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+\frac{1}{2}\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{z1}\left[E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+1\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)\right]+{\gamma }_{y1}\left[-E_z^{n-\left(l-1\right)/m}\left(i,j+1,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z2}\left[E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k+2\right)-E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k-1\right)\right]+{\gamma }_{y2}\left[-E_z^{n-\left(l-1\right)/m}\left(i,j+2,k+\frac{1}{2}\right)+E_z^{n+\left(l-1\right)/m}\left(i,j-1,k+\frac{1}{2}\right)\right]\right\}$ $H_y^{n-l/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)=H_y^{n-\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+\frac{1}{2}\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{x1}\left[E_z^{n+\left(l-1\right)/m}\left(i+1,j+\frac{1}{2},k\right)-E_z^{n+\left(l-1\right)/m}\left(i,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z1}\left[-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+1\right)+E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)\right]+{\gamma }_{x2}\left[E_z^{n+\left(l-1\right)/m}\left(i+2,j,k+\frac{1}{2}\right)-E_z^{n+\left(l-1\right)/m}\left(i-1,j,k+\frac{1}{2}\right)\right]+{\gamma }_{z2}\left[-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k+2\right)+E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k-1\right)\right]\right\}$ $H_z^{n-l/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)=H_z^{n-\left(l-1\right)/m}\left(i+\frac{1}{2},j+\frac{1}{2},k\right)+\frac{c_l\Delta t}{{\mu }_0}\times \left\{{\gamma }_{y1}\left[E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j+1,k\right)-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j,k\right)\right]+{\gamma }_{x1}[-E_y^{n+\left(l-1\right)/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n+\left(l-1\right)/m}\left(i,j+\frac{1}{2},k\right)+{\gamma }_{y2}\left[E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j+2,k\right)-E_x^{n+\left(l-1\right)/m}\left(i+\frac{1}{2},j-1,k\right)\right]+{\gamma }_{x2}\left[-E_y^{n+\left(l-1\right)/m}\left(i+2,j+\frac{1}{2},k\right)+E_y^{n+\left(l-1\right)/m}\left(i-1,j+\frac{1}{2},k\right)\right]\right\}$ $J_x^{n+\frac{l}{m}}\left(i+\frac{1}{2},j,k\right)=\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+$ $\left[\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0{\omega }_p^2\{\frac{v\left[1-\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}+$ $\frac{{\omega }_b\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}\}{E}_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0{\omega }_p^2\{\frac{v\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}+$ $\frac{{\omega }_b\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_{l}d}-1\right]}{v^2+{\omega }_b^2}\}{E}_y^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)$

• • because a current density position is the same as an electric field position, therefore, when calculating

$J_x^{n+l/m}\left(i+\frac{1}{2},j,k\right),$

• •  it is necessary to perform interpolation on

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)\mathrm{and}{E}_y^{n-l/m}\left(i+\frac{1}{2},j,k\right),$

• •  because J_y and E_y are no value at a point

$\left(i+\frac{1}{2},j,k\right),$

• •  so it needs to be averaged by four adjacent diagonal values to obtain:

$J_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(J_y^{n-l/m}\left(i,j-\frac{1}{2},k\right)+J_y^{n-l/m}\left(i,j+\frac{1}{2},k\right)+J_y^{n-l/m}\left(i+1,j+\frac{1}{2},k\right)+J_y^{n-l/m}\left(i+1,j-\frac{1}{2},k\right)\right)$ $E_y^{n-l/m}\left(i+\frac{1}{2},j,k\right)=\frac{1}{4}\left(E_y^{n-l/m}\left(i,j-\frac{1}{2},k\right)+E_y^{n-l/m}\left(i,j+\frac{1}{2},k\right)+E_y^{n-l/m}\left(i+1,j+\frac{1}{2},k\right)+E_y^{n-l/m}\left(i+1,j-\frac{1}{2},k\right)\right)$ $J_y^{n+\frac{l}{m}}\left(i,j+\frac{1}{2},k\right)=\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_y^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)+$ $\left[\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]J_x^{n+\frac{l-1}{m}}\left(i+\frac{1}{2},j,k\right)+{\epsilon }_0{\omega }_p^2\{\frac{v\left[1-\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}\right]}{v^2+{\omega }_b^2}+$ $\frac{{\omega }_b\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}\}{E}_y^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)-{\epsilon }_0{\omega }_p^2\{\frac{v\sin \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}}{v^2+{\omega }_b^2}+$ $\frac{{\omega }_b\left[\cos \left({\omega }_b{c}_l\mathrm{dt}\right)e^{-{\mathrm{vc}}_l\mathrm{dt}}-1\right]}{v^2+{\omega }_b^2}\}{E}_x^{n+\frac{l-1}{m}}\left(i,j+\frac{1}{2},k\right)$

• • because the current density position is the same as the electric field position, therefore, when calculating

$J_y^{n+l/m}\left(i,j+\frac{1}{2},k\right),$

• •  it is necessary to perform interpolation on

$J_x^{n-l/m}\left(i,j+\frac{1}{2},k\right)\mathrm{and}{E}_x^{n+l/m}\left(i,j+\frac{1}{2},k\right),$

• •  because J_x and E_x are no value at a point

$\left(i,j+\frac{1}{2},k\right),$

• •  so it needs to be averaged by four adjacent diagonal values to obtain:

$J_x^{n-l/m}\left(i,j+\frac{1}{2},k\right)=\frac{1}{4}\left(J_x^{n-l/m}\left(i-\frac{1}{2},j,k\right)+J_x^{n-l/m}\left(i+\frac{1}{2},j,k\right)+J_x^{n-l/m}\left(i+\frac{1}{2},j+1,k\right)+J_x^{n-l/m}\left(i-\frac{1}{2},j+1,k\right)\right)$ $E_x^{n+l/m}\left(i,j+\frac{1}{2},k\right)=\frac{1}{4}\left(E_x^{n-l/m}\left(i-\frac{1}{2},j,k\right)+E_x^{n-l/m}\left(i+\frac{1}{2},j,k\right)+E_x^{n-l/m}\left(i+\frac{1}{2},j+1,k\right)+E_x^{n-l/m}\left(i-\frac{1}{2},j+1,k\right)\right)$ $J_z^{n+\frac{l}{m}}\left(i,j,k+\frac{1}{2}\right)=e^{-{\mathrm{vc}}_l\mathrm{dt}}{J}_z^{n+\frac{l-1}{m}}\left(i,j,k+\frac{1}{2}\right)-\frac{e^{-{\mathrm{vc}}_l\mathrm{dt}}-1}{v}{E}_z^{n+\frac{l-1}{m}}\left(i,j,k+\frac{1}{2}\right)$

• • wherein, i, j, k represent spatial nodes of electric field, magnetic field, and current density:

${\gamma }_{x1}=\frac{9}{8}\frac{1}{\Delta x},{\gamma }_{y1}=\frac{9}{8}\frac{1}{\Delta y},{\gamma }_{z1}=\frac{9}{8}\frac{1}{\Delta z},{\gamma }_{x2}=\frac{-1}{24}\frac{1}{\Delta x},{\gamma }_{y2}=\frac{-1}{24}\frac{1}{\Delta y},{\gamma }_{z2}=\frac{-1}{24}\frac{1}{\Delta z}.$

Optionally, the electromagnetic model includes an anisotropic magnetized plasma slab model, a blunt cone model, and a sphere model.

The various embodiments in the present disclosure are described in a progressive manner, and each embodiment focuses on the differences from other embodiments. The same and similar parts between each embodiment can be referred to each other. For the system disclosed in the embodiments, due to its correspondence with the methods disclosed in the embodiments, the description is relatively simple, please refer to the method section for relevant information.

The present disclosure has provided specific examples to explain the principles and implementation methods of the present disclosure, and the above embodiments are only used to help understand the methods and core ideas of the present disclosure. At the same time, for general skilled person in the art, there may be changes in specific implementation methods and application scope based on the ideas of the present disclosure. In summary, the content of this disclosure should not be understood as a limitation of the present disclosure.

Claims

1. A method for processing anisotropic magnetized plasma medium, comprising: obtaining a Maxwell equation and a polarization current density equation based on electromagnetic characteristics of anisotropic magnetized plasma;processing the Maxwell equation and the polarization current density equation to obtain an electric field intensity, a magnetic field intensity, and a polarization current density after processed;using a matrix exponential time-domain finite difference method to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed;applying a numerical modeling simulation electromagnetic model to determine electromagnetic characteristics of the electromagnetic model according to the numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in the anisotropic magnetized plasma medium.

2. The method for processing anisotropic magnetized plasma medium according to claim 1, wherein the Maxwell equation is:

3. The method for processing anisotropic magnetized plasma medium according to claim 2, wherein processing the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed, specifically comprising: performing multi-level symplectic discretization on the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density in a symplectic discretization format:

4. The method for processing anisotropic magnetized plasma medium according to claim 1, wherein using a matrix exponential time-domain finite difference method to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed specifically comprises: assuming that the plasma medium is biased by a static magnetic field in the z-direction, based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, the matrix exponential time-domain finite difference method is used to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium:

5. The method for processing anisotropic magnetized plasma medium according to claim 1, wherein the electromagnetic model comprises an anisotropic magnetized plasma slab model, a blunt cone model, and a sphere model.

6. A system for processing anisotropic magnetized plasma medium, comprising: a Maxwell equation and polarization current density equation determination module, configured to obtain a Maxwell equation and a polarization current density equation based on electromagnetic characteristics of anisotropic magnetized plasma;a Maxwell equation and polarization current density equation processing module, configured to process the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density after processed;a matrix exponential time-domain finite difference processing module, configured to obtain numerical iterative equations for electric field intensity, magnetic field intensity, and polarization current density in anisotropic magnetized plasma medium using the matrix exponential time-domain finite difference method based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed; andan electromagnetic characteristic determination module, configured to determine electromagnetic characteristics of the electromagnetic model using a numerical modeling simulation electromagnetic model based on the numerical iterative equations of the electric field intensity, the magnetic field intensity, and the polarization current density in the anisotropic magnetized plasma medium.

7. The system for processing anisotropic magnetized plasma medium according to claim 6, wherein the Maxwell equation is:

8. The system for processing anisotropic magnetized plasma medium according to claim 6, wherein the Maxwell equation and polarization current density equation processing module comprises: a Maxwell equation and polarization current density equation processing unit, configured to perform multi-level symplectic discretization on the Maxwell equation and the polarization current density equation to obtain the electric field intensity, the magnetic field intensity, and the polarization current density in the symplectic discretization format:

9. The system for processing anisotropic magnetized plasma medium according to claim 8, wherein the matrix exponential time-domain finite difference processing module comprises: a Matrix exponential time-domain finite difference processing unit, assuming that the plasma medium is biased by a static magnetic field in the z-direction, based on the electric field intensity, the magnetic field intensity, and the polarization current density after processed, the matrix exponential time-domain finite difference method is used to obtain numerical iterative equations for the electric field intensity, the magnetic field intensity, and the polarization current density in anisotropic magnetized plasma medium:

10. The system for processing anisotropic magnetized plasma medium according to claim 6, wherein the electromagnetic model comprises an anisotropic magnetized plasma slab model, a blunt cone model, and a sphere model.