求指教,总是提示内部矩阵维度必须一致,但不知道怎么改?谢谢了!
%%%%国际单位制,研究贝塞尔高斯关联涡旋光束的微粒操控
clear all;clc;
tic
sigma=2*1e-3;f=3.0001*1e-3;lambda=0.6328*1e-6;bs=2*pi/lambda;
m=0;%先设置为零,模拟完全部内容后,更改m=1,2,3,研究不同拓扑荷对横向,轴向梯度力的影响,
%%%研究梯度力对轴向散射力的影响
c=3*10e8;%光速度
a=3*1e-9;%粒子半径
n1=1.8;%粒子折射率
n2=1.33;%周围环境折射率
mm=n1/n2;%粒子的相对折射率
nr=0.0000020001;%相当于2um,可调整
nz=0.0000030001;%相当于3um,可调整
delta=0.5*1e-3;beta=3; %可调整,
k1=20;l1=20;%可调整,可通过增大k1,l1来调整精度,增大到模拟出的结果几乎没变化
%x=0.0000000001;%模拟x轴的情况时的轴设置,
y=0.0000000001;%模拟x轴的情况时的轴设置,
x=-nr:nr/20:nr;%模拟z轴的情况时的轴设置,
z=f+0.0000000001;%模拟z轴的情况时的轴设置
P=1;
A1=P/(2*pi*(sigma*1e2));
rho=sqrt(x.^2+y.^2);
%z=f+0.0000000001:nz/20:(f+nz);
AA=1-z/f; BB=f; CC=-1/f;DD=0;
zeta=1/(4*sigma^2)+1/(2*delta^2)+1i*bs*AA./(2*BB);
eta=1/(4*sigma^2)+1/(2*delta^2)-1i*bs*AA./(2*BB);
It0=0;
for k=-k1:1:k1;
for l=-l1:1:l1;
for r1=0:.00005:0.008;
for r2=0:.00005:0.008;
It0=It0+A1*(bs./BB).^2.*exp(-zeta.*r1.^2-eta.*r2.^2)...
.*besselj(k,beta/delta.*r1).*besselj(k,beta/delta.*r2)...
.*besseli(abs(k+l+m),r1.*r2/delta^2)...
.*(1/2*besselj(l,bs*r1.*rho/BB).*(besselj(-1+l,bs*r2.*rho/BB)...
-besselj(1+l,bs*r2.*rho/BB))*bs*r2./BB*sign(x)...
+1/2*(besselj(-1+l,bs*r1.*rho/BB)-besselj(1+l,bs*r1.*rho/BB))...
.*bs*r1/BB*sign(x).*besselj(l,bs*r2.*rho/BB)).*r1.*r2*0.0000000025;
end
end
end
end
%It=abs(It0);
%bb=max(aa1(1:6));
%plot(z,It/bb);
%aa1=It/bb;
Fg=2*pi*n2*a^3/c*(mm^2-1)/(mm^2+2)*It0*10^12;
%plot(z,Fg)%沿z轴的变化画图
plot(x,real(Fx))%沿x轴的变化画图
mytimer1=toc;
disp(mytimer1)
%%%%国际单位制,研究贝塞尔高斯关联涡旋光束的微粒操控
clear all;clc;
tic
sigma=2*1e-3;f=3.0001*1e-3;lambda=0.6328*1e-6;bs=2*pi/lambda;
m=0;%先设置为零,模拟完全部内容后,更改m=1,2,3,研究不同拓扑荷对横向,轴向梯度力的影响,
%%%研究梯度力对轴向散射力的影响
c=3*10e8;%光速度
a=3*1e-9;%粒子半径
n1=1.8;%粒子折射率
n2=1.33;%周围环境折射率
mm=n1/n2;%粒子的相对折射率
nr=0.0000020001;%相当于2um,可调整
nz=0.0000030001;%相当于3um,可调整
delta=0.5*1e-3;beta=3; %可调整,
k1=20;l1=20;%可调整,可通过增大k1,l1来调整精度,增大到模拟出的结果几乎没变化
%x=0.0000000001;%模拟x轴的情况时的轴设置,
y=0.0000000001;%模拟x轴的情况时的轴设置,
x=-nr:nr/20:nr;%模拟z轴的情况时的轴设置,
z=f+0.0000000001;%模拟z轴的情况时的轴设置
P=1;
A1=P/(2*pi*(sigma*1e2));
rho=sqrt(x.^2+y.^2);
%z=f+0.0000000001:nz/20:(f+nz);
AA=1-z/f; BB=f; CC=-1/f;DD=0;
zeta=1/(4*sigma^2)+1/(2*delta^2)+1i*bs*AA./(2*BB);
eta=1/(4*sigma^2)+1/(2*delta^2)-1i*bs*AA./(2*BB);
It0=0;
for k=-k1:1:k1;
for l=-l1:1:l1;
for r1=0:.00005:0.008;
for r2=0:.00005:0.008;
It0=It0+A1*(bs./BB).^2.*exp(-zeta.*r1.^2-eta.*r2.^2)...
.*besselj(k,beta/delta.*r1).*besselj(k,beta/delta.*r2)...
.*besseli(abs(k+l+m),r1.*r2/delta^2)...
.*(1/2*besselj(l,bs*r1.*rho/BB).*(besselj(-1+l,bs*r2.*rho/BB)...
-besselj(1+l,bs*r2.*rho/BB))*bs*r2./BB*sign(x)...
+1/2*(besselj(-1+l,bs*r1.*rho/BB)-besselj(1+l,bs*r1.*rho/BB))...
.*bs*r1/BB*sign(x).*besselj(l,bs*r2.*rho/BB)).*r1.*r2*0.0000000025;
end
end
end
end
%It=abs(It0);
%bb=max(aa1(1:6));
%plot(z,It/bb);
%aa1=It/bb;
Fg=2*pi*n2*a^3/c*(mm^2-1)/(mm^2+2)*It0*10^12;
%plot(z,Fg)%沿z轴的变化画图
plot(x,real(Fx))%沿x轴的变化画图
mytimer1=toc;
disp(mytimer1)