yuanzhus.m

是matlab语言

clear all
n1=1.5075;
n2=1.5;
r=9.8*10^(-6);
dn=(n1-n2)/n1;
for ii=1:101
    r1(ii)=(r/100)*ii;
end
for ii=1:60
    n(ii)=n1*sqrt((1-2*dn*(r1(ii)/r1(60))));
end
for ii=61:101
    n(ii)=n2;
end
rh=9.8*10^(-8);
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^2)*(n(ii)^2-(n0)^2)-2*(2-m^2/2);
end%系数b
for ii=1:101
        f(ii,1)=(rh^2)*(k0^2)*(n(ii)^2-(n0)^2)-(rh^2)*(m^2)/r^2-2;
end%系数f
a=4*(k0^2)*(n0^2)*(rh^2)/(1+i*k0*n0*zh);
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%系数整理结束
 %x=0;%光开始输入
 %h1=0.03; 
 %for jj=1:101
  %  u(1,jj)=exp(-((x+jj)*h1)^2);
 %end%光结束
 r0=2.9815*10^(-6);
for ii=1:101
    r(ii)=rh*(ii-1);
end
for jj=1:101
    u(1,jj)=0.25*10^6*r(jj)*exp(-3.33*10^5*(r(jj)-r0))^2;
end
for m1=1:101%大程序开始
A=zeros(101,101);
A(1,1)=a+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,1)+a;
end
for ii=2:99
    A(ii,ii+1)=g(ii,1);
end%左面的系数矩阵A结束
B=zeros(101,101);
B(1,1)=a+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,1)*d+a;
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+b(1,1));
    for jj=2:100
    h2(m1,jj)=g(jj,1)./(f(1,1)+a-e(jj,1)*h2(m1,jj-1));
    end
    y(m1,1)=d2(m1,1)/(a+b(1,1));
    for jj=2:101
    y(m1,jj)=(d1(jj,m1)-e(jj,1).*y(m1,jj-1))./(f(1,1)+a-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));
    surf(abs(q));
    shading interp