# MATLAB实现聚焦光束仿真_单点聚焦matlab

MATLAB实现聚焦光束仿真

一.基本概念

二. 实验仿真计算

• 聚焦光束的理论实现

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

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

  $$\mathbf{E}\left(\mathbf{r}\right)=-\frac{ikf}{2\pi }\underset{0}{\overset{{\theta }_{max}}{\int }}\underset{0}{\overset{2\pi }{\int }}\mathbf{A}\left(\theta ,\varphi \right)\exp \left(i\mathbf{k}\cdot \mathbf{r}\right)\sin \theta d\varphi d\theta$$

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

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

  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));
      else
          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));
      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
  1
  2
  3
  4
  5
  6
  7
  8
  9
  10
  11
  12
  13
  14
  15
  16
  17
  18
  19
  20
  21
  22
  23
  24
  25
  26
  27
  28
  29
  30
  31
  32
  33
  34
  35
  36
  37
  38
  39
  40
  41
  42
  43
  44
  45
  46
  47
  48
  49
  50
  51
  52
  53
  54
  55
  56
  57