98 lines
3.1 KiB
Python
98 lines
3.1 KiB
Python
clc; clear; close all;
|
||
|
||
% ================= 1. 物理参数定义 =================
|
||
r = 0.1; % 圆柱半径
|
||
eps_r = 5.0; % 相对介电常数
|
||
m = sqrt(eps_r); % 相对折射率 m = ~1.414
|
||
k0 = 50; % 背景真空中波数 (k=6)
|
||
k1 = k0 * m; % 圆柱内部波数
|
||
x_size = k0 * r; % 尺寸参数 x = k0*a
|
||
|
||
% ================= 2. 计算域网格设置 =================
|
||
x_range = 1;
|
||
y_range = 1;
|
||
Nx = 500;
|
||
Ny = 500;
|
||
x_vec = linspace(0, x_range, Nx);
|
||
y_vec = linspace(0, y_range, Ny);
|
||
[X, Y] = meshgrid(x_vec, y_vec);
|
||
|
||
xc = 0.5; yc = 0.5;
|
||
[Phi, R] = cart2pol(X - xc, Y - yc); % 转换为极坐标
|
||
|
||
% ================= 3. 场初始化 =================
|
||
E_scat = zeros(size(X)); % 散射场
|
||
E_int = zeros(size(X)); % 内部场
|
||
|
||
% Wiscombe 截断准则（决定级数展开需要算到第几阶）
|
||
N_trunc = round(x_size + 4.05 * x_size^(1/3) + 2);
|
||
|
||
% ================= 4. 2D Mie 级数展开计算 =================
|
||
% 2D 圆柱级数从 -N 到 +N
|
||
for n = -N_trunc : N_trunc
|
||
|
||
% 边界处的贝塞尔函数值
|
||
J_nx = besselj(n, x_size);
|
||
J_nmx = besselj(n, k1 * r);
|
||
H_nx = besselh(n, 1, x_size);
|
||
|
||
% 边界处的导数值 (利用递推公式 Z_n' = 0.5 * (Z_{n-1} - Z_{n+1}))
|
||
J_nx_p = 0.5 * (besselj(n-1, x_size) - besselj(n+1, x_size));
|
||
J_nmx_p = 0.5 * (besselj(n-1, k1*r) - besselj(n+1, k1*r));
|
||
H_nx_p = 0.5 * (besselh(n-1, 1, x_size) - besselh(n+1, 1, x_size));
|
||
|
||
% 计算 TM 偏振下的散射系数 a_n (对应 E_z)
|
||
num_a = m .* J_nx .* J_nmx_p - J_nx_p .* J_nmx;
|
||
den_a = J_nmx .* H_nx_p - m .* J_nmx_p .* H_nx;
|
||
a_n = num_a ./ den_a;
|
||
|
||
% 计算内部透射系数 c_n
|
||
num_c = J_nx .* H_nx_p - J_nx_p .* H_nx; % 这其实是 Wronskian
|
||
c_n = num_c ./ den_a;
|
||
|
||
% 空间相位因子: i^n * exp(i*n*phi)
|
||
phase = (1i)^n * exp(1i * n * Phi);
|
||
|
||
% 累加外部散射场 (仅在 R >= r 区域有效)
|
||
out_idx = R >= r;
|
||
E_scat(out_idx) = E_scat(out_idx) + a_n .* besselh(n, 1, k0 * R(out_idx)) .* phase(out_idx);
|
||
|
||
% 累加内部总场 (仅在 R < r 区域有效)
|
||
in_idx = R < r;
|
||
E_int(in_idx) = E_int(in_idx) + c_n .* besselj(n, k1 * R(in_idx)) .* phase(in_idx);
|
||
end
|
||
|
||
% ================= 5. 组装全场并绘图 =================
|
||
% 入射平面波: u_inc = exp(i*k0*x)
|
||
phase_shift = exp(1i * k0 * xc);
|
||
E_scat = E_scat .* phase_shift;
|
||
E_int = E_int .* phase_shift;
|
||
|
||
E_inc = exp(1i * k0 * X);
|
||
|
||
% 总场 = 外部(入射 + 散射) + 内部场
|
||
% 组装总场
|
||
E_total = zeros(size(X));
|
||
E_total(R >= r) = E_inc(R >= r) + E_scat(R >= r);
|
||
E_total(R < r) = E_int(R < r);
|
||
%
|
||
% % 提取最大场强做对比
|
||
% max_E_val = max(abs(E_total(:)));
|
||
% fprintf('2D 理论解析解中心区域最大场强 (max |E_total|): %.4f\n', max_E_val);
|
||
|
||
% 绘图
|
||
figure('Color','w');
|
||
pcolor(X, Y, abs(E_total-E_inc));
|
||
max_E_real = max(max(abs(E_total-E_inc)));
|
||
shading interp;
|
||
axis equal tight;
|
||
colorbar;
|
||
colormap jet;
|
||
title(sprintf('2D Cylinder Mie Scattering |E_{scatter}| (Max = %.4f)', max_E_real));
|
||
|
||
% 绘制圆柱边界
|
||
hold on;
|
||
theta_circle = linspace(0, 2*pi, 100);
|
||
plot(xc + r * cos(theta_circle), yc + r * sin(theta_circle), 'k--', 'LineWidth', 1.5);
|
||
hold off;
|