在FDTD中绘制某个点的场强随波长变化的曲线的脚本

以下是在MATLAB中绘制某个点的场强随波长变化的FDTD脚本示例：

%% 参数设置
c = 299792458; % 光速
f0 = c/1e-6; % 中心频率
lambda0 = 1e-6; % 中心波长
n = 1.5; % 介质折射率
dx = lambda0/20; % 空间步长
dt = dx/c/2; % 时间步长
tmax = 600*dt; % 最大模拟时间
x = -10*dx:dx:10*dx; % 空间网格
y = -10*dx:dx:10*dx;
[X, Y] = meshgrid(x, y); % 空间网格矩阵
z = dx/2; % 某个点的位置
epsilon = n^2*ones(size(X)); % 介质介电常数

%% 激励光源
source = zeros(size(X));
source(size(X,1)/2,size(X,2)/2) = 1;

%% 初始化场分量
Ex = zeros(size(X));
Ey = zeros(size(X));
Hz = zeros(size(X));

%% FDTD时间步进
for t = 1:round(tmax/dt)
    % 更新Ex、Ey分量
    Ex(:,2:end-1) = Ex(:,2:end-1) + dt/(dx*epsilon(:,2:end-1)).*(Hz(:,2:end-1)-Hz(:,1:end-2));
    Ey(2:end-1,:) = Ey(2:end-1,:) - dt/(dx*epsilon(2:end-1,:)).*(Hz(2:end-1,:)-Hz(1:end-2,:));
    % 更新Hz分量
    Hz(2:end-1,2:end-1) = Hz(2:end-1,2:end-1) + dt/(dx*n).^2.*(Ex(2:end-1,3:end)-Ex(2:end-1,2:end-1)+Ey(3:end,2:end-1)-Ey(2:end-1,2:end-1));
    % 激励光源
    Hz(z,z) = Hz(z,z) + sin(2*pi*f0*t*dt);
end

%% 绘制场强随波长变化曲线
lambda = linspace(0.5, 1.5, 101)*lambda0; % 波长范围
freq = c./lambda; % 频率范围
E = zeros(size(lambda)); % 场强
for k = 1:length(lambda)
    n_lambda = n*ones(size(X)); % 介质折射率
    n_lambda(abs(X)<dx/2 & abs(Y)<dx/2) = n/sqrt(epsilon(z,z)); % 在某个点处改变折射率
    dx_lambda = lambda(k)/20; % 新的空间步长
    dt_lambda = dx_lambda/c/2; % 新的时间步长
    tmax_lambda = 600*dt_lambda; % 新的最大模拟时间
    x_lambda = -10*dx_lambda:dx_lambda:10*dx_lambda; % 新的空间网格
    y_lambda = -10*dx_lambda:dx_lambda:10*dx_lambda;
    [X_lambda, Y_lambda] = meshgrid(x_lambda, y_lambda); % 新的空间网格矩阵
    epsilon_lambda = n_lambda.^2*ones(size(X_lambda)); % 新的介质介电常数
    % 初始化场分量
    Ex_lambda = zeros(size(X_lambda));
    Ey_lambda = zeros(size(X_lambda));
    Hz_lambda = zeros(size(X_lambda));
    % FDTD时间步进
    for t_lambda = 1:round(tmax_lambda/dt_lambda)
        % 更新Ex、Ey分量
        Ex_lambda(:,2:end-1) = Ex_lambda(:,2:end-1) + dt_lambda/(dx_lambda*epsilon_lambda(:,2:end-1)).*(Hz_lambda(:,2:end-1)-Hz_lambda(:,1:end-2));
        Ey_lambda(2:end-1,:) = Ey_lambda(2:end-1,:) - dt_lambda/(dx_lambda*epsilon_lambda(2:end-1,:)).*(Hz_lambda(2:end-1,:)-Hz_lambda(1:end-2,:));
        % 更新Hz分量
        Hz_lambda(2:end-1,2:end-1) = Hz_lambda(2:end-1,2:end-1) + dt_lambda/(dx_lambda*n_lambda).^2.*(Ex_lambda(2:end-1,3:end)-Ex_lambda(2:end-1,2:end-1)+Ey_lambda(3:end,2:end-1)-Ey_lambda(2:end-1,2:end-1));
        % 激励光源
        Hz_lambda(abs(X_lambda)<dx_lambda/2 & abs(Y_lambda)<dx_lambda/2) = Hz(z,z);
    end
    % 计算场强
    E(k) = abs(Hz_lambda(z,z));
end
figure;
plot(freq, E);
xlabel('频率 (Hz)');
ylabel('场强');

在脚本中，首先设置了模拟的参数，包括光速、中心频率、中心波长、介质折射率、空间步长、时间步长、最大模拟时间、空间网格等。然后，定义了激励光源并初始化场分量。接着，进行FDTD时间步进，直到达到最大模拟时间。在某个点处改变介质折射率，然后重新进行FDTD模拟，计算该点处的场强，并绘制场强随波长变化曲线。