function Q = CalculateQ( t,R1,R2,n1,n2,n3,m, lambdaf )
dl=0.00005;%dl需要根据Q值变化做调整，调整方法看figure1
lambda1=lambdaf-dl;
lambda2=lambdaf+dl;
lambda=linspace(lambda1,lambda2,10000);
k=2*pi./lambda;
SS=length(k);
Fm=ones(SS,1);
a1=besselj(m,n1.*R1.*k);
a2=n1*(besselj(m-1,n1.*R1.*k)-besselj(m+1,n1.*R1.*k))/2;
b1=besselj(m,n2.*R1.*k);
b2=n2*(besselj(m-1,n2.*R1.*k)-besselj(m+1,n2.*R1.*k))/2;
b3=besselj(m,n2.*R2.*k);
b4=n2*(besselj(m-1,n2.*R2.*k)-besselj(m+1,n2.*R2.*k))/2;
c1=besselh(m,n2.*R1.*k);
c2=n2*(besselh(m-1,n2.*R1.*k)-besselh(m+1,n2.*R1.*k))/2;
c3=besselh(m,n2.*R2.*k);
c4=n2*(besselh(m-1,n2.*R2.*k)-besselh(m+1,n2.*R2.*k))/2;
fm1=besselh(m,n3.*R2.*k);
fm2=n3*(besselh(m-1,n3.*R2.*k)-besselh(m+1,n3.*R2.*k))/2;
for i=1:SS
AA=[-a1(i),b1(i),c1(i),0;-a2(i),b2(i),c2(i),0;0,b3(i),c3(i),-fm1(i);0,b4(i),c4(i),-fm2(i)];
BB=[0;0;-a1(i);-a2(i)];
% x=AA\BB;
[L,U]=lu(AA);
x=U\(L\BB);
Fm(i)=x(4);
end
d=abs(Fm).^2;
dd=abs(d-max(d)/2);
B= find(diff(sign(diff(dd)))==2)+1;
halfright=lambda(B(2));
halfleft=lambda(B(1));
Q=lambdaf/(halfright-halfleft);
figure(1)
plot(lambda,d)
hold on
%红点和蓝点对应半峰宽位置
plot(halfright,d(B(2)),'O','MarkerSize',8,'MarkerEdgeColor','r','MarkerFaceColor','r');
hold on
plot(halfleft,d(B(1)),'O','MarkerSize',8,'MarkerEdgeColor','b','MarkerFaceColor','b');
hold off
ylabel('散射场强度')
xlabel('wavelength (\mum)')
end
