二维Z-FDTD代码C++

以下是一个基于C++的二维Z-FDTD代码示例：

#include <iostream>
#include <cmath>

using namespace std;

const double c = 299792458.0; // 光速
const double mu0 = 4 * M_PI * 1e-7; // 真空磁导率
const double eps0 = 1 / (c * c * mu0); // 真空电容率

const int nx = 100; // x方向网格数
const int ny = 100; // y方向网格数
const double dx = 1e-3; // x方向网格间距
const double dy = 1e-3; // y方向网格间距
const double dt = dx / (2 * c); // 时间步长
const double T = 1e-8; // 模拟时长
const int nt = T / dt; // 时间步数

double Ez[nx][ny] = {0}; // 电场
double Hy[nx][ny] = {0}; // 磁场

double source(double t) { // 激励源
    return exp(-(t - 5e-9) * (t - 5e-9) / (2 * 1e-9 * 1e-9));
}

int main() {
    for (int n = 0; n < nt; ++n) { // 时间步循环
        // 更新Ez
        for (int i = 1; i < nx - 1; ++i) {
            for (int j = 1; j < ny - 1; ++j) {
                Ez[i][j] = Ez[i][j] + (dt / eps0) * ((Hy[i][j] - Hy[i - 1][j]) / dx - (Hy[i][j] - Hy[i][j - 1]) / dy);
            }
        }
        // 加入激励源
        Ez[nx / 2][ny / 2] = Ez[nx / 2][ny / 2] + source(n * dt);
        // 更新Hy
        for (int i = 1; i < nx - 1; ++i) {
            for (int j = 1; j < ny - 1; ++j) {
                Hy[i][j] = Hy[i][j] + (dt / mu0) * ((Ez[i][j] - Ez[i][j + 1]) / dy - (Ez[i + 1][j] - Ez[i][j]) / dx);
            }
        }
    }
    return 0;
}

该代码使用了Z-FDTD方法模拟了一个正方形区域内的电磁波传播情况。其中，Ez和Hy分别表示电场和磁场，通过时间步循环和Maxwell方程的离散化，不断更新它们的值，并最终输出结果。代码中的source函数表示激励源，它在模拟过程中为电场添加了一个高斯型脉冲信号。