MATLAB实现聚焦光束仿真

一.基本概念

二. 实验仿真计算

聚焦光束的理论实现

聚焦光束的实现主要是由于基于Debye理论得到，聚焦光束由于其特殊的物理性质在现在科学中有着不可替代的作用，我们现在所熟知的光镊技术，其基础就在于适用聚焦光束对微纳粒子进行捕获操作，其光束的聚焦基本图示如下：

假设现在输入的光束的是一个径向极化的光束，此时在进行紧聚焦之后，根据Debye理论可以得到聚焦之后的电场表达式为

E(r)=−ikf2π∫0θmax∫02πA(θ,ϕ)exp⁡(ik⋅r)sin⁡θdϕdθ\boldsymbol{E}\left( \boldsymbol{r} \right) =-\frac{ikf}{2\pi}\int\limits_0^{\theta _{max}}{\int\limits_0^{2\pi}{\boldsymbol{A}\left( \theta ,\phi \right) \exp \left( i\boldsymbol{k}\cdot \boldsymbol{r} \right) \sin \theta d\phi d\theta}} E(r)=−2πikf0∫θmax0∫2πA(θ,ϕ)exp(ik⋅r)sinθdϕdθ

同样的，对于磁场的表达式，只需要根据麦克斯韦方程中的磁场和电场之间的关系就可以得到，这里我就不在列出(不要问我为什么还有磁场，因为光波是电磁场，包括电场和磁场分量)。

当我们完成上式的计算之后，我们就可以对公式进行仿真计算，这里我采用的是径向极化的高斯光束作为入射光束，当然还存在圆极化，线极化，方位角极化等情况，不同极化方式的光束只是在光束的表达式上存在差异，但是具体的推导过程都是一样的，所以不需要进行过多的考虑。

function [E_rho,E_phi,E_z]=opticalfield(m,phi,r,z)
lambda0=1.064e-6;               %入射光波波长
n_i=1.518;                      %折射率
NA=1.32;                        %
power=0.1;                      % 入射光线的功率
f=1.2e-3;                       % 透镜的聚焦焦距
k0=2*pi/lambda0;                %波数
k=k0*n_i;
A=k*f/2;
theta_min=0;
theta_max=asin(NA/n_i);
beta0=1.5;
E_rho1=zeros(size(phi,1),size(phi,2)); %预分配内存
E_phi1=zeros(size(phi,1),size(phi,2));
E_z1=zeros(size(phi,1),size(phi,2));theta1=34.66*pi/180;
theta2=43.88*pi/180;
theta3=51.94*pi/180;
theta4=57.89*pi/180;   %角度系数% theta1=0;
% theta2=0;
% theta3=0;
% theta4=0;   %角度系数
theta_num=100;
dtheta=(theta_max-theta_min)/theta_num;
theta=theta_min:dtheta:theta_max;
for i=1:length(theta)thetai=ones(size(phi))*theta(i);if (theta(i)>theta1&&theta(i)<=theta2)||(theta(i)>theta3&&theta(i)<=theta4)E_rho1(:,:,i)=-A.*1i^(m).*exp(1i*m*phi).*power.*sqrt(cos(thetai)).*sin(2*thetai).*exp(-beta0^2*(sin(thetai)).^2./(sin(theta_max))^2).*besselj(m,2*beta0*sin(thetai)/...sin(theta_max)).*exp(1i*k*z.*cos(thetai)).*(besselj(m+1,k*r.*sin(thetai))-besselj(m-1,k*r.*sin(thetai)));E_phi1(:,:,i)=A.*1i^(m+1).*exp(1i*m*phi).*power.*sqrt(cos(thetai)).*sin(thetai).*exp(-beta0^2*(sin(thetai)).^2./(sin(theta_max))^2).*besselj(m,2*beta0*sin(thetai)/...sin(theta_max)).*exp(1i*k*z.*cos(thetai)).*(besselj(m+1,k*r.*sin(thetai))+besselj(m-1,k*r.*sin(thetai)));E_z1(:,:,i)=-2*A.*1i^(m+1).*exp(1i*m*phi).*power.*sqrt(cos(thetai)).*(sin(thetai)).^2.*exp(-beta0^2*(sin(thetai)).^2./(sin(theta_max))^2).*besselj(m,2*beta0*sin(thetai)/...sin(theta_max)).*exp(1i*k*z.*cos(thetai)).*besselj(m,k*r.*sin(thetai));elseE_rho1(:,:,i)=A.*1i^(m).*exp(1i*m*phi).*power.*sqrt(cos(thetai)).*sin(2*thetai).*exp(-beta0^2*(sin(thetai)).^2./(sin(theta_max))^2).*besselj(m,2*beta0*sin(thetai)/...sin(theta_max)).*exp(1i*k*z.*cos(thetai)).*(besselj(m+1,k*r.*sin(thetai))-besselj(m-1,k*r.*sin(thetai)));E_phi1(:,:,i)=-A.*1i^(m+1).*exp(1i*m*phi).*power.*sqrt(cos(thetai)).*sin(thetai).*exp(-beta0^2*(sin(thetai)).^2./(sin(theta_max))^2).*besselj(m,2*beta0*sin(thetai)/...sin(theta_max)).*exp(1i*k*z.*cos(thetai)).*(besselj(m+1,k*r.*sin(thetai))+besselj(m-1,k*r.*sin(thetai)));E_z1(:,:,i)=2*A.*1i^(m+1).*exp(1i*m*phi).*power.*sqrt(cos(thetai)).*(sin(thetai)).^2.*exp(-beta0^2*(sin(thetai)).^2./(sin(theta_max))^2).*besselj(m,2*beta0*sin(thetai)/...sin(theta_max)).*exp(1i*k*z.*cos(thetai)).*besselj(m,k*r.*sin(thetai));end
end
content1=floor(theta1/dtheta);
content2=floor(theta2/dtheta);
content3=floor(theta3/dtheta);
content4=floor(theta4/dtheta);
E_rho=sum(E_rho1(:,:,[1:1:content1])*dtheta,3)+sum(E_rho1(:,:,[1:1:content2])*dtheta,3)+...sum(E_rho1(:,:,[1:1:content3])*dtheta,3)+sum(E_rho1(:,:,[1:1:content4])*dtheta,3)+sum(E_rho1*dtheta,3);
E_phi=sum(E_phi1(:,:,[1:1:content1])*dtheta,3)+sum(E_phi1(:,:,[1:1:content2])*dtheta,3)+...sum(E_phi1(:,:,[1:1:content3])*dtheta,3)+sum(E_phi1(:,:,[1:1:content4])*dtheta,3)+sum(E_phi1*dtheta,3);
E_z=sum(E_z1(:,:,[1:1:content1])*dtheta,3)+sum(E_z1(:,:,[1:1:content2])*dtheta,3)+...sum(E_z1(:,:,[1:1:content3])*dtheta,3)+sum(E_z1(:,:,[1:1:content4])*dtheta,3)+sum(E_z1*dtheta,3);%计算各个方向的光场
end