光学紧聚焦与光力计算M代码

%% 光学紧聚焦与光力计算 - 优化版本
% 作者: Hamilton
% 日期: 2025-03-25
% 功能: 计算高NA物镜紧聚焦矢量光束的光场分布及对金纳米粒子的光学力

%% 1. 初始化与参数设置
clear all; close all; clc;
tic; % 开始计时

% 光学参数
n = 1.332;          % 介质折射率(水)
lambda = 1047e-9;   % 入射波长[m](1047nm近红外)
k = 2 * pi * n / lambda; % 介质中的波数
NA = 0.95 * n;      % 数值孔径(NA ≤ n)
alfy = asin(NA/n);  % 物镜最大半孔径角[rad]

% 光束参数
P = 0.5;            % 入射功率[W]
f = 3e-3;           % 物镜焦距[m]
w0 = 1e-4/sqrt(5.29); % 入射光束束腰半径[m](5.29为经验参数)
E0 = sqrt(P/(pi*w0^2)); % 入射电场振幅[V/m]

% 金纳米粒子参数
a = 19.1e-9;        % 金粒子半径[m]
eps0 = 8.854e-12;   % 真空介电常数[F/m]
eps1 = n^2;         % 介质相对介电常数
eps2 = -54 + 5.9i;  % 金的复介电常数(λ=1047nm)

%% 2. 计算网格设置
grid_size = 500;    % 网格点数(原1000，优化为500提高速度)
r = linspace(-3*lambda, 3*lambda, grid_size); % 径向坐标[m]
z = linspace(-3*lambda, 3*lambda, grid_size); % 轴向坐标[m]
[rr, zz] = meshgrid(r, z); % 创建二维网格

%% 3. 角度分段参数(用于Richards-Wolf积分)
Nx = 500;           % 角度采样点数
x1 = 5.43*pi/180;   % 第一区间上限[rad]
x2 = 35.36*pi/180;  % 第二区间上限[rad]

% 角度采样
theta = linspace(0, alfy, Nx); % 均匀角度采样
dtheta = theta(2) - theta(1);  % 角度步长

% 分段掩码
mask1 = theta <= x1;          % 第一角度区间(0 ≤ θ ≤ x1)
mask2 = theta > x1 & theta <= x2; % 第二区间(x1 < θ ≤ x2)
mask3 = theta > x2;           % 第三区间(x2 < θ ≤ alfy)

%% 4. 预计算公共项(提高计算效率)
sin_theta = sin(theta);
cos_theta = sqrt(cos(theta)); % 切趾因子
exp_term = exp(-(sin_theta/sin(alfy)).^2); % 高斯截断项
phase_factor = exp(1i*k*zz.*cos(theta)); % 相位因子(矩阵化)

% 贝塞尔函数项
J1_kr = besselj(1, k*rr.*sin_theta); % J1(kr sinθ)
J1_2beta = besselj(1, 2*sin_theta/sin(alfy)); % J1(2β)
J0_kr = besselj(0, k*rr.*sin_theta); % J0(kr sinθ)

%% 5. 紧聚焦场计算(矢量衍射积分)
% 初始化场矩阵
Er = zeros(grid_size, grid_size); % 径向电场
Ez = zeros(grid_size, grid_size); % 纵向电场
Hphi = zeros(grid_size, grid_size); % 角向磁场

% 分段积分计算
for i = 1:Nx
    % 公共积分核
    integrand_base = sin_theta(i)/sin(alfy)^2 * cos_theta(i) * exp_term(i) * J1_2beta(i);
    
    % 径向分量Er
    term_Er = integrand_base * sin(2*theta(i)) * J1_kr(:,:,i);
    % 纵向分量Ez
    term_Ez = 2 * integrand_base * sin_theta(i)^2 * J0_kr(:,:,i);
    % 角向磁场Hphi
    term_Hphi = integrand_base * sin_theta(i) * J1_kr(:,:,i);
    
    % 相位项(矩阵化)
    phase = phase_factor(:,:,i);
    
    % 分段加权
    if mask1(i)
        weight = dtheta;
    elseif mask2(i)
        weight = -dtheta; % 中间区间符号反转
    else
        weight = dtheta;
    end
    
    % 累加各场分量
    Er = Er + 1i * term_Er .* phase * weight;
    Ez = Ez + 1i * term_Ez .* phase * weight;
    Hphi = Hphi + 1i * term_Hphi .* phase * weight;
end

% 场幅度归一化
eta = E0 * pi * f * n / lambda;
Er = eta * Er;
Ez = eta * Ez;
Hphi = eta * Hphi;

%% 6. 光强计算
Ir = abs(Er).^2;       % 径向偏振光强
Iz = abs(Ez).^2;       % 纵向偏振光强
I_total = Ir + Iz;     % 总光强
I_normalized = I_total / max(I_total(:)); % 归一化光强

%% 7. 金纳米粒子光学力计算
% 极化率
gamma = 4*pi*a^3 * eps1*(eps2 - eps1) / (eps2 + 2*eps1);

% 光学截面
C_abs = k * n * imag(gamma) / eps1; % 吸收截面
C_sca = k^4 * abs(gamma)^2 / (6*pi); % 散射截面

% 坡印廷矢量(轴向和径向)
Sz = real(Er .* conj(Hphi)); % 轴向分量
Sr = -real(Ez .* conj(Hphi)); % 径向分量(注意负号)
S = sqrt(Sz.^2 + Sr.^2);     % 总坡印廷矢量大小

% 光学力
F_abs = n * S * C_abs / c;   % 吸收力
F_sca = n * S * C_sca / c;   % 散射力
F_total = (F_abs + F_sca) / 2; % 总光学力(近似)

% 梯度力
dr = r(2) - r(1);            % 径向步长
dz = z(2) - z(1);            % 轴向步长
[gradIr, gradIz] = gradient(I_total, dr, dz); % 光强梯度
F_grad = real(gamma) * eps0 * gradIr / 4; % 梯度力

%% 8. 可视化
% 光强分布
figure(1);
imagesc(r/lambda, z/lambda, I_normalized');
axis equal tight;
xlabel('径向位置 r/λ');
ylabel('轴向位置 z/λ');
title('归一化光强分布');
colorbar;
colormap hot;

% 光学力分布
figure(2);
subplot(1,2,1);
plot(r/lambda, F_total(grid_size/2,:), 'LineWidth', 2);
xlabel('径向位置 r/λ');
ylabel('光学力大小');
title('径向光学力分布');

subplot(1,2,2);
plot(z/lambda, F_total(:,grid_size/2), 'LineWidth', 2);
xlabel('轴向位置 z/λ');
ylabel('光学力大小');
title('轴向光学力分布');

toc; % 显示总计算时间