Python和Matplotlib电磁一维二维三维模拟

自由空间的时间相关麦克斯韦旋度方程为

∂ E ∂ t = 1 ε 0 ∇ × H ∂ H ∂ t = − 1 μ 0 ∇ × E . begin{aligned}frac{partial boldsymbol{E}}{partial t} &=frac{1}{varepsilon_{0}} nabla times boldsymbol{H} \frac{partial boldsymbol{H}}{partial t} &=-frac{1}{mu_{0}} nabla times boldsymbol{E} .end{aligned} ∂t∂E​∂t∂H​​=ε0​1​∇×H=−μ0​1​∇×E.​

E 和 H 是三个维度的向量，因此，一般来说，等式 (1.1a) 和 (1.1b) 分别代表三个方程。 我们将从仅使用 Ex 和 Hy 的简单一维案例开始，因此 等式(1.1a) 和 (1.1b) 变为

∂ E x ∂ t = − 1 ε 0 ∂ H y ∂ z , ∂ H y ∂ t = − 1 μ 0 ∂ E x ∂ z . begin{aligned}&frac{partial E_{x}}{partial t}=-frac{1}{varepsilon_{0}} frac{partial H_{y}}{partial z}, \&frac{partial H_{y}}{partial t}=-frac{1}{mu_{0}} frac{partial E_{x}}{partial z} .end{aligned} ​∂t∂Ex​​=−ε0​1​∂z∂Hy​​,∂t∂Hy​​=−μ0​1​∂z∂Ex​​.​

这些是沿z方向行驶的平面波的方程，其电场在X方向上定向，并以y方向定向磁场。

对时间和空间导数进行中心差分近似，给出

E x n + 1 / 2 ( k ) − E x n − 1 / 2 ( k ) Δ t = − 1 ε 0 H y n ( k + 1 2 ) − H y n ( k − 1 2 ) Δ x H y n + 1 ( k + 1 2 ) − H y n ( k + 1 2 ) Δ t = − 1 μ 0 E x n + 1 / 2 ( k + 1 ) − E x n + 1 / 2 ( k ) Δ x begin{aligned}frac{E_{x}^{n+1 / 2}(k)-E_{x}^{n-1 / 2}(k)}{Delta t} &=-frac{1}{varepsilon_{0}} frac{H_{y}^{n}left(k+frac{1}{2}right)-H_{y}^{n}left(k-frac{1}{2}right)}{Delta x} \frac{H_{y}^{n+1}left(k+frac{1}{2}right)-H_{y}^{n}left(k+frac{1}{2}right)}{Delta t} &=-frac{1}{mu_{0}} frac{E_{x}^{n+1 / 2}(k+1)-E_{x}^{n+1 / 2}(k)}{Delta x}end{aligned} ΔtExn+1/2​(k)−Exn−1/2​(k)​ΔtHyn+1​(k+21​)−Hyn​(k+21​)​​=−ε0​1​ΔxHyn​(k+21​)−Hyn​(k−21​)​=−μ0​1​ΔxExn+1/2​(k+1)−Exn+1/2​(k)​​

import numpy as np
from math import exp
from matplotlib import pyplot as plt

ke = 200
ex = np.zeros(ke)
hy = np.zeros(ke)

# Pulse parameters
kc = int(ke / 2)
t0 = 40
spread = 12

nsteps = 100

# Main FDTD Loop
for time_step in range(1, nsteps + 1):

    # Calculate the Ex field
    for k in range(1, ke):
        ex[k] = ex[k] + 0.5 * (hy[k - 1] - hy[k])

    # Put a Gaussian pulse in the middle
    pulse = exp(-0.5 * ((t0 - time_step) / spread) ** 2)
    ex[kc] = pulse

    # Calculate the Hy field
    for k in range(ke - 1):
        hy[k] = hy[k] + 0.5 * (ex[k] - ex[k + 1])

# Plot the outputs as shown in Fig. 1.2
plt.rcParams['font.size'] = 12
plt.figure(figsize=(8, 3.5))

plt.subplot(212)
plt.plot(hy, color='k', linewidth=1)
plt.ylabel('H$_y$', fontsize='14')
plt.xlabel('FDTD cells')
plt.xticks(np.arange(0, 201, step=20))
plt.xlim(0, 200)
plt.yticks(np.arange(-1, 1.2, step=1))
plt.ylim(-1.2, 1.2)

plt.subplots_adjust(bottom=0.2, hspace=0.45)
plt.show()