%% 电偶极辐射仿真（有限差分法） clear; clc; close all; %% 参数设置 lambda = 0.5; % 波长(m) c = 3e8; % 光速(m/s) f = c / lambda; % 频率(Hz) omega = 2*pi*f; % 角频率 p0 = 1e-9; % 电偶极矩振幅(C·m) sigma = 0.1; % 脉冲宽度(s) t0 = 0.5/sigma; % 脉冲中心时间(s) dx = lambda/50; % 空间步长(m) dy = dx; % 二维仿真时dy=dx dz = dx; % 三维仿真时dz=dx dt = dx/c/1.5; % 时间步长(s)（满足CFL条件） x_range = [-1 1]; % 计算区域x范围(m) y_range = [-1 1]; % 计算区域y范围(m) z_range =[0 0] ; % 二维仿真时z=0 epsilon0 = 8.854187817e-12; % 真空介电常数(F/m) mu0 = 4*pi*1e-7; % 真空磁导率(H/m) %% 初始化网格 [x,y,z] = meshgrid(... x_range(1):dx:x_range(2), ... y_range(1):dy:y_range(2), ... z_range(1):dz:z_range(2)); N_x = size(x,1); % 网格点数x方向 N_y = size(y,2); % 网格点数y方向 N_z = size(z,3); % 网格点数z方向 %% 初始化场量 Ex = zeros(N_x,N_y,N_z); % x方向电场 Ey = zeros(N_x,N_y,N_z); % y方向电场 Hz = zeros(N_x,N_y,N_z); % z方向磁场 %% 边界条件设置（PML吸收边界） pml_order = 4; % PML阶数 pml_thickness = 10; % PML厚度(网格单元数) % 计算PML参数（这里仅简化实现） sigma_x = zeros(N_x,N_y,N_z); sigma_y = zeros(N_x,N_y,N_z); for i = 1:pml sigma_x(i,:,:) = i/pml * 1e-3; sigma_x(N_x-i+1,:,:) = i/pml* 1e-3; end %% 时间迭代参数 t_max = 2*t0; % 总仿真时间(s) n_t = round(t_max/dt); % 总时间步数 E_record= cell(n_t,1); % 存储电场数据 %% 仿真主循环 for n = 1:n_t t = n*dt; % 当前时间 % 电偶极激励（高斯脉冲） p = p0 * exp(-(t - t0)^2/(2*sigma^2)); % 电场更新（FDTD公式） Hz(2:end-1,2:end-1,2:end-1) = ... Hz(2:end-1,2:end-1,2:end-1) ... + (dt/(mu0*dy)) * (Ey(2:end-1,3:end,2:end-1) - Ey(2:end-1,2:end-1,2:end-1)) ... - (dt/(mu0*dx)) * (Ex(3:end,2:end-1,2:end-1) - Ex(2:end-1,2:end-1,2:end-1)); % 电场更新（考虑电偶极源） Ex(floor(N_x/2),floor(N_y/2),floor(N_z/2)) = ... Ex(floor(N_x/2),floor(N_y/2),floor(N_z/2)) ... + p/(epsilon0*dx*dy*dz) * dt; % 边界处理（PML） Hz = Hz .* exp(-sigma_x*dt/(2*mu0)); % 记录电场数据 E_record{n} = abs(Ex(:,:,1)); % 记录 % 可视化更新 if mod(n,50) == 0 figure(1); surf(x(:,:,1), y(:,:,1), E_record{n}, 'EdgeColor','none'); shading interp; title(['t = ',num2str(t,'%.2f'),' s']); axis([-1 1 -1 1 0 max(E_record{n}(:))]); drawnow; end end %% 三维场分布可视化 figure(2); p = isosurface(x,y,z,abs(Ex),max(abs(Ex(:)))*0.5); isonormals(x,y,z,abs(Ex),p); patch(p,'FaceColor','interpolated','EdgeColor','none'); isonormals(x,y,z,abs(Ex),p); % 优化表面法线 daspect([1 1 1 ]); view(45,30); % 设置视角 camlight('headlight'); % 添加光源 lighting phong; xlabel('x (m)'); ylabel('y (m)'); zlabel('z (m)'); title('三维电场分布'); colormap jet; %% 近场特性分析 figure(3); near_field = Ex(floor(N_x/2)-10:floor(N_x/2)+10, ... floor(N_y/2)-10:floor(N_y/2)+10, ... floor(N_z/2)); surf(abs(near_field)); title('近场分布'); xlabel('x方向'); ylabel('y方向'); zlabel('电场强度');中的for i = 1:pml显示未定义函数或变量 'pml'，请告诉我怎么样修改