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;