format long
clear
n1=1.45
n2=1.447
a=4.5
k=(a^2)*(n1^2-n2^2);
u1=(a^2)*(n1^2);
w1=(a^2)*(n2^2);
XiShu1=(n1^2+n2^2)/(2*n1^2);
XiShu2=(n1^2-n2^2)/(2*n1^2);
i=1;
n=n2;
for V=0:0.001:1
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.000001:n1
%n
nn=n+0.000001
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2)
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2)
if(U==0 || W==0)
break;
%disp('sss')
else
z1=besselj(0,U)/(U*besselj(1,U));
z2=1/(U^2)+XiShu1*(1/W^2+besselk(0,W)/(W*besselk(1,W)))-sqrt((XiShu2*(1/W^2+besselk(0,W)/(W*besselk(1,W))))^2+(BeiTa/n1)^2*(V^2/(U2*W2))^2);
if(abs(z1-z2)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
%disp('nnn')
n=BeiTa;
break;
end
end
end
end
for V=1:0.01:6
k02=V^2/k;
k01=sqrt(k02);
for BeiTa=n:0.00001:n1
%n
U2=u1*k02-(a^2)*k02*(BeiTa^2);
U=sqrt(U2);
W2=(a^2)*k02*(BeiTa^2)-w1*k02;
W=sqrt(W2);
if(U==0 || W==0)
break;
else
z1=besselj(0,U)/(U*besselj(1,U));
z2=1/(U^2)+XiShu1*(1/W^2+besselk(0,W)/(W*besselk(1,W)))-sqrt((XiShu2*(1/W^2+besselk(0,W)/(W*besselk(1,W))))^2+(BeiTa/n1)^2*(V^2/(U2*W2))^2);
if(abs(z1-z2)<0.01)
x(i)=V;
y(i)=BeiTa;
i=i+1;
z1;
z2;
n=BeiTa;
break;
end
end
end
end
plot(x,y)
save SL_V.mat x
save SL_N.mat y
评论1