opticalduijies.m

是matlb语言

clear all
n1=1.5075;
n2=1.5;
r1=9.8*10^(-6);
r2=4.3*10^(-6);
for jj=1:51
    r(jj)=r1;
    rh(jj)=r1/100;
end
for jj=52:101
    r(jj)=r2;
    rh(jj)=r2/100;
end
%r2=2.5*10^(-5);
%注意
dn=(n1-n2)/n1;
for ii=1:101
    r1(ii)=(r(ii)/100)*ii;
end
n=n2.*ones(101,1);
for ii=1:60
    n(ii,1)=n1*sqrt((1-2*dn*(r1(ii)/r1(60))));
end
%[x,y]=meshgrid(1:101,1:101);mesh(n);
%for i=1:101 for j=1:101 if r1(i,j)<=r2 n(i,j)=n2; end end end

k0=2*pi/(1.55*10^(-6));
m=1;
zh=20*10^(-8);
n0=1.503;
for ii=1:101
    b(ii,1)=(k0^2)*(rh(1)^2)*(n(ii,1)^2-(n0)^2)-2*(2-m^2/2);
end%系数b
for jj=1:101
for ii=1:101
        f(ii,jj)=(rh(jj)^2)*(k0^2)*(n(ii,1)^2-(n0)^2)-(rh(jj)^2)*(m^2)/r(jj)^2-2;
end
end%系数f
for jj=1:101
a(jj)=4*(k0^2)*(n0^2)*(rh(jj)^2)/(1+i*k0*n0*zh);
end
c=(1+exp(i*m*pi))*(2-m^2/2);
d=(1-k0*n0*i*zh)/(1+k0*n0*i*zh);
e=zeros(101,1);
g=zeros(101,1);
for ii=2:101
    e(ii,1)=1-1/(2*(ii-1));
end
for ii=2:101
      g(ii,1)=1+1/(2*(ii-1));  
end%系数整理结束
 r0=2.9815*10^(-6);
for ii=1:101
    r1(ii)=rh(1)*(ii-1);
end
for jj=1:101
    u(1,jj)=0.25*10^6*r1(jj)*exp(-3.33*10^5*(r1(jj)-r0))^2;
end
 
for m1=1:101%大程序开始
A=zeros(101,101);
A(1,1)=a(m1)+b(1,1);
A(1,2)=c;
for ii=2:101
    A(ii,ii-1)=e(ii,1);
end
for ii=2:100
    A(ii,ii)=f(ii,m1)+a(m1);
end
for ii=2:99
    A(ii,ii+1)=g(ii,1);
end%左面的系数矩阵A结束
B=zeros(101,101);
B(1,1)=a(m1)+b(1,1)*d;
B(1,2)=c*d;
for ii=2:101
    B(ii,ii-1)=e(ii,1)*d;
end
for ii=2:100
    B(ii,ii)=f(ii,m1)*d+a(m1);
end
for ii=2:99
    B(ii,ii+1)=g(ii,1)*d;
end%右面的系数矩阵B结束
d1=ones(101,101);
y=ones(101,101);
%z=ones(101,101);
d1(:,m1)=B*u.';%转置
  d2=ones(101,101);
for ii=1:101
    d2(m1,ii)=d1(ii,m1);%变成一行一行的
end
    %下面用追赶法
  h2(m1,1)=c./(a(m1)+b(1,1));
    for jj=2:100
    h2(m1,jj)=g(jj,1)./(f(1,m1)+a(m1)-e(jj,1)*h2(m1,jj-1));
    end
    y(m1,1)=d2(m1,1)/(a(m1)+b(1,1));
    for jj=2:101
    y(m1,jj)=(d1(jj,m1)-e(jj,1).*y(m1,jj-1))./(f(1,m1)+a(m1)-e(jj,1).*h2(m1,jj-1));%problem
    end
    z(m1,101)=y(m1,101);
    for jj=100:-1:1
    z(m1,jj)=y(m1,jj)-h2(m1,jj).*z(m1,jj+1);
    end
    for jj=1:101
        u(1,jj)=z(m1,jj);
    end
    for jj=1:101
    q(m1,jj)=z(m1,jj);
    end
    end%最大循环的end
    [x,y]=meshgrid(0:1:100,0:zh:99*zh);
    %mesh(abs(q));
    mesh(abs(q));
    surf(abs(q));
    shading interp